; Procedure to compare SLODAR and SODAR profiles
; Input: month and date of the evening, hour and min. to start binning, duration of the interval  
;
; Output: plot of the profiles
; Jan 15, 2005: new sodar data, use selsod.pro
;-------------------------------------------------

pro slosod, mm, dd, uth, utmin, dutmin

; 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) )
neve = n_elements(eve) 
nslo = neve 

if (nslo eq 0) then return
prof = slodat[0]
out = replicate(prof,nslo)
out(0:neve-1) = slodat(eve)


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

sod = selsod(mm, dd)
nsod = n_elements(sod)

utsod = sod.uthour - 1. 

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

nalt = 40

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


; select and bin the Sodar profiles 
sodprof =  sod[0].cn2*0.
sodalt = sod[0].alt
n = 0
utstart = uth + utmin/60.
utend = utstart  + dutmin/60.

for i=0,nsod-1 do  if ( ( utsod[i] ge utstart)  and ( utsod[i] le utend) ) then begin
   sodprof += sod[i].cn2*1e13   ; in 1e-13 m^1/3 units
   n += 1
;   print, 'UT ', utsod[i]
endif
sodprof = sodprof/n
   print, n, ' SODAR pofiles in the selected UT interval'

sel = where( (utslo gt utstart - 1./3.) and (utslo le utend) )

  if (sel[0] lt 0) then begin
     print, 'No SLODAR matching data! Exiting...'
     return
  endif


; Average the slodar profiles (wide binary) over selected sampling time
sloalt = fltarr(8) & sloprof =  fltarr(8) & nprof=0
slo = out[sel]
nslo = n_elements(slo)
for i=0,nslo-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 1) then begin
     print, 'No SLODAR high-res. data! Exiting...'
     return
endif
  print, nprof, ' SLODAR profiles averaged'

sloalt = sloalt/nprof
sloprof = sloprof/nprof/100.
 
sodbin = fltarr(8)
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[j] += sodprof[i] 


sc1 = max(sloprof)/dh*1000.
sc2 = max(sodprof)/20.*1000.
sc = max([sc1,sc2])*1.2
  xlab = max(sloalt)*0.6


!p.charsize = 1.5
!x.style=1
!y.style=1
;  plot, sloalt, sloprof, xtitle='Altitude, m', ytitle='Integral, 10!U-13!N m!U1/3',yr=[-0.1*sc,sc]
;  oplot, sloalt, sloprof, psym=2  
; oplot, sloalt, sodprof/20.*dh, linestyle=2

  plot, sloalt, sloprof/dh*1000., xtitle='Altitude, m', ytitle='C!Dn!U2!N,  10!U-16!N m!U-2/3',yr=[-0.1*sc,sc]
  oplot, sloalt, sloprof/dh*1000, psym=2  
 oplot, sodalt, sodprof/20.*1000., linestyle=2

  xyouts, xlab,sc*0.9, string(dd,format='(I2)')+'/'+  string(mm,format='(I2)')+' 2004'
  xyouts, xlab,sc*0.7, 'UT '+string(uth,format='(I2)')+'h'+  string(utmin,format='(I2)')+'m, '+string(dutmin,format='(I2)')+'min'


;stop

end