; Procedure to compare SLODAR and SODAR profiles
; Input: month and date of the evening
pro slosod, mm, dd

; slodar data selection
restore, 'slodar.idl'
; ---------slodar structure:
; prof = { UTyear: 0, UTmonth: 0, UTday: 0, UT: 0.0, $
;	  r0: 0.0, alt:FLTARR(9), Intcn2:FLTARR(9) }
; wrong (evening) dates in the CN2 files!
eve = where((slodat.UTmonth eq mm) and (slodat.UTday eq dd) )
morn = -1
nmorn = 0

neve = n_elements(eve) 
nmorn =  n_elements(morn)
nslo = neve + nmorn

if (nslo eq 0) then return
prof = slodat[0]
out = replicate(prof,nslo)

out(0:neve-1) = slodat(eve)
;out(neve:nslo-1) = slodat(morn)

utslo = fltarr(nslo)
for i=0,nslo-1 do utslo[i] = out[i].UT

print, 'Total SLODAR profiles found: ', nslo
print, 'UT range: ',utslo[0], utslo[nslo-1]
nalt = 8
tmp = fltarr(nslo,nalt) 
for i=0,nslo-1 do tmp(i,*) = out[i].intcn2(0:nalt-1)

; --------------------  SODAR data selection
 restore, 'sodar1.idl'
; prof = { UTyear: 0, UTmonth: 0, UTday: 0, UThour:0.0, UTmin:0.0, $
;	    nlevels:40, alt:FLTARR(40), cn2:FLTARR(40) }

eve = where((spdat.UTmonth eq mm) and (spdat.UTday eq dd) and (spdat.UThour gt 12) )

dd1 = dd+1   & mm1 = mm
if (dd1 gt 30) then begin dd1=1 & mm1=mm+1 & endif

morn = where((spdat.UTmonth eq mm1) and (spdat.UTday eq dd1) and (spdat.UThour lt 12) )

neve = 0 & nmorn = 0

if (eve[0] ge 0) then neve = n_elements(eve)
if (morn[0] ge 0) then nmorn =  n_elements(morn)
nsod = neve + nmorn

prof = spdat[0]
sod = replicate(prof,nsod)


if (eve[0] ge 0) then sod(0:neve-1) = spdat(eve)
if (morn[0] ge 0) then sod(neve:nsod-1) = spdat(morn)
nalt = 40

utsod = fltarr(nsod)
for i=0,nsod-1 do utsod[i] = sod[i].UThour ; +sod[i].UTmin/60.
; utsod(0:neve-1) =  utsod(0:neve-1)-24.
 utsod = utsod -1.  ; Correct UT error in sdar.txt

print, 'Total SODAR profiles found: ', nsod
print, 'UT range: ',utsod[0], utsod[nsod-1]

;----------------------------------------------------
; Data analysis
;-------------------------------------------------------


; Make running plot
jsod = -1

; outer loop on SODAR profile

sodbin = fltarr(nsod, 8)
slobin = fltarr(nsod, 8)
utbin = fltarr(nsod)

for nn=0,nsod-1 do begin

; select the Sodar profile number nn
  ut = utsod[nn]
  sodprof = sod[nn].cn2*1e13   ; in 1e-13 m^1/3 units
  sodalt = sod[nn].alt -20.
  sel = where( (utslo ge ut - 1./3.) and (utslo le ut) )

  if (sel[0] lt 0) then continue

; Average the slodar profiles (wide binary) over 20-min. SODAR integration time
  sloalt = fltarr(8) & sloprof =  fltarr(8) & nprof=0
  slo = out[sel]
  for i=0,n_elements(slo)-1 do begin
    dh = slo[i].alt[1] - slo[i].alt[0]
    if (dh le 500.) then begin
      nprof = nprof+1
      sloalt = sloalt + slo[i].alt[0:7]
      sloprof = sloprof + slo[i].intcn2[0:7]
    endif
  endfor
  
  if (nprof lt 2) then continue

  jsod = jsod+1
  sloalt = sloalt/nprof
  slobin(jsod,*) = sloprof/nprof/100.
  utbin[jsod] = utsod(nn) - 1./6.
 
  dh = sloalt[1]-sloalt[0]
  for i=0,37 do for j=0,7 do if (abs(sloalt[j] - sodalt[i]) le 0.5*dh) then sodbin[jsod,j] += sodprof[i] 

endfor

utbin = utbin(0:jsod)
slobin = slobin(0:jsod,*)
sodbin = sodbin(0:jsod,*)


!p.charsize = 1.5
!x.style=1
sh = 5

plot, utbin, slobin(*,0), xtitle='UT, hours', ytitle='Integral, 10!U-13!N m!U1/3',yr=[0,4*sh]
oplot, utbin, sodbin(*,0), linestyle=2
  oplot, utbin, slobin(*,0), psym=2


for j=1,3 do begin
  oplot, utbin, slobin(*,j)+sh*j
  oplot, utbin, sodbin(*,j)+sh*j, linestyle=2
  oplot, utbin, slobin(*,j)+sh*j, psym=2
  oplot, utbin, utbin*0.+sh*j, linestyle=1

endfor

  xyouts, utbin[2],sh*3.5, +string(dd,format='(I2)')+'/'+string(mm,format='(I2)')
; stop

end