; Procedure to compare SLODAR and MASS-DIMM profiles
; Input: month and date of the evening
; Nov 30, 2004
;---------------------------------------------------------
pro mdseeing, 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) )
neve = n_elements(eve) 
morn = -1
nmorn = 0

; Correct UT dates
;eve = where((slodat.UTmonth eq mm) and (slodat.UTday eq dd) and (slodat.UThour gt 12) )
;dd1 = dd+1   & mm1 = mm
;if (dd1 gt 30) then begin dd1=1 & mm1=mm+1 & endif
;morn = where((slodat.UTmonth eq mm1) and (slodat.UTday eq dd1) and (slodat.UThour lt 12) )

;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)

;--------MASS-DIMM data --------------------------------
 restore, 'massdimm.idl'
; mddat structure:
; prof = { UTyear: 0, UTmonth: 0, UTday: 0, UTmass:0.0, $
;	  fsee:0.0, Intcn2:FLTARR(6), UTdimm:0.0, see:0.0 }

;eve = where((mddat.UTmonth eq mm) and (mddat.UTday eq dd) and (mddat.UTmass gt 12) )
dd1 = dd+1   & mm1 = mm
if (dd1 gt 31) then begin dd1=1 & mm1=mm+1 & endif

morn = where((mddat.UTmonth eq mm1) and (mddat.UTday eq dd1) and (mddat.UTmass lt 12) )

;neve = n_elements(eve)
neve = 0
nmorn =  n_elements(morn)
nmd = neve + nmorn

prof = mddat[0]
md = replicate(prof,nmd)
; md(0:neve-1) = mddat(eve)
md(neve:nmd-1) = mddat(morn)

utmd = fltarr(nmd)
for i=0,nmd-1 do utmd[i] = md[i].UTmass


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

hmass = [0,0.5,1,2,4,8,16]*1000.   ; MASS layers incl. ground
nmass = 7

nalt = 128      ; rebin altitude grid
altmax = 32000.
altmin = 0.  
altstep = (altmax- altmin)/nalt 
realt = findgen(nalt)*altstep + altmin

; MASS triangular weights
wt = fltarr(nmass,nalt)
fact = alog10(2.)
for i=1,nmass-1 do   wt(i,*) = (1. - abs(alog10(realt/hmass[i])/fact) ) > 0.
for i=1,nalt-1 do if (realt[i] gt 16000) then wt(nMASS-1,i) = 1.
wt(0,*) = 1. - total(wt(1:nmass-1,*),1)
;stop

utslo = out.ut


;------- spatial binning of SLODAR data
sloprof = fltarr(nslo,nmass)

  for i=0,nslo-1 do begin
   dh = out[i].alt[1] -  out[i].alt[0]
   if (dh gt 500.) then begin
  ; 
     nrebin = fix(8*dh/altstep)          
     tmpprof = congrid((out[i].intcn2(0:7)*0.01 > 0.), nrebin)
     if (nrebin gt nalt) then print, 'SLODAR profile too high!!!'
       for k=0,nmass-1 do  sloprof(i,k) = total(tmpprof*wt(k,0:nrebin-1))*(altstep/dh) 
;     print, 'Profile processed at UT ',out[i].UT
;     print, 'Intergals raw, rebinned: ', total(out[i].intcn2(0:7)*0.01),total( sloprof(i,*))
;       plot, hmass, sloprof(i,*)*2., psym=2,yr=[0,3]
;     dt = abs(utmd - out[i].ut)
;     sel = where (dt eq min(dt))
;      mdprof = fltarr(nmass)
;      mdprof(1:6) = md[sel].intcn2
;      mdprof[0] =  ( (md[i].see^(5./3.) - md[i].fsee^(5./3.) ) > 0) *6.8
;       oplot, hmass, mdprof, linestyle=0
;       print, 'MASS-DIMM integral: ', total(mdprof)
;     stop
   endif  
endfor

print, 'Spatial binning done!'
;stop

; temporal binning
utstart = 1.
utend = 9.
utstep = 0.2
nut = fix(utend - utstart)/utstep
rebinut = findgen(nut)*utstep + utstart

reslo = fltarr(nut,nmass)
remd = fltarr(nut,nmass)
numslo = intarr(nut) & nummd = intarr(nut)

for k=0,nut-1 do begin
  for i=0,nslo-1 do begin
  if ((abs(out[i].ut - rebinut[k]) le utstep*0.5) and (sloprof[i,0] gt 0.)) then begin
    numslo[k] += 1
    reslo(k,*) += sloprof(i,*)
  endif
  endfor

  for i=0,nmd-1 do begin
  if (abs(md[i].UTmass - rebinut[k]) le utstep*0.5)  then begin
    nummd[k] += 1
    remd(k,1:6) += md[k].intcn2
    remd[k,0] += ( (md[i].see^(5./3.) - md[i].fsee^(5./3.) ) > 0) *6.8
  endif
  endfor

endfor

sel = where ((nummd gt 1) and (numslo gt 1))
if (sel[0] lt 0) then begin
  print, 'No useful data found!'
  return
endif

rebinut = rebinut(sel)
nummd = nummd(sel) & numslo = numslo(sel)
nut = n_elements(rebinut)
for k=0,nut-1 do begin
  remd(k,*) = remd(sel[k],*)/nummd[k]  
  reslo(k,*) = reslo(sel[k],*)/numslo[k]
endfor

remd = remd(0:nut-1,*)
reslo = reslo(0:nut-1,*)

  ; write data to the file
  name = '04'+string(mm, format='(I2)')+string(dd, format='(I2)')+'mdseeing.dat'
openw, 1, name
print, 'Writing data to the file ',name
for k=0,nut-1 do begin
    printf, 1, rebinut[k], numslo[k], nummd[k],reslo(k,*), remd(k,*),$
    format='(F5.2,2I3,14F6.2)'
endfor
close, 1

!p.charsize = 1.5
!x.style=1
for k=0,nut-1 do begin
  plot, 1e-3*hmass, remd(k,*), linestyle=0, xr=[-0.5,18], xtitle='Altitude, km'
  oplot, 1e-3*hmass,reslo(k,*), psym=2
  print, 'UT, nslo, nmd;: ', rebinut[k], numslo[k], nummd[k]
  stop
endfor  

  print, 'No more profiles for this date1'
stop


;stop
end