;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Extrapolation of minimum r0, L0, noise for DIM struct. functions
; Input:
;    arrays of variances, mult. and add. constants, L0s
; 
FUNCTION MINIMUM_PARAMETERS, sig, Cs, as, L0s

  N_L = N_ELEMENTS(L0s)

  mn = (WHERE(sig eq min(sig)))[0]

  if ((mn eq N_L-1) or (mn eq 0)) then begin
     L0 = L0s[mn]
     C = Cs[mn]
  endif else begin
     x = sig[mn-1] - sig[mn+1]
     x = x / max([sig[mn-1] - sig[mn], sig[mn+1] - sig[mn]])
     x = x / 2.
     L0 = L0s[mn] * (L0s[mn+1]/L0s[mn])^x
     if (x gt 0) then begin
        C = Cs[mn] + (Cs[mn+1] - Cs[mn]) * x
     endif else begin
        C = Cs[mn] + (Cs[mn] - Cs[mn-1]) * x
     endelse
  endelse

  r0 = C^(-3./5.)
  
  sigma = as[mn]/2.		; division by 2 to get noise
  if (sigma gt 0) then begin
     sigma=sqrt(sigma)
  endif else begin
     sigma=0
  endelse

  RETURN, [r0,L0,sigma]
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Calculating of r0, PSD etc. from Zernike decomposition
; Input:
;    phase : 2-dim array of phases
;	(either phase or cent_x and cent_y must be given)
;    disp1 : if set, display phase PSD
;    disp2 : if set, display Zernike PSDs
;    disp3 : if set, display Zernike variances
;    disp4 : if set, display von Karman Zernike residuals
;    disp5 : if set, display r0s from Zernikes (Kolmogorov) with error bars
;    bhw : 1 if Blackmann-Harris window is to be used in PSD calcs
;    
; Output: structure of type
;	res ={	r0:		FLTARR(N_orders), 	$
;			r0's from each Zernike order
;		dr0:		FLTARR(N_orders),	$
;			dr0's from each Zernike order
;		pow_ph:		FLTARR(N_frames/2+1),	$
;			phase PSD
;		sigz:		FLTARR(N_terms),	$
;			Zernike variances sig_i^2
;		pow_Z:		FLTARR(N_frames/2+1,N_terms)}
;			Zernike mode power spectra
; 
FUNCTION CALC_R0_ZERNIKE,phase  
  eff_tel_dia=8.992             ; =16.*0.562
  
  disp1=0
  disp2=0
  disp3=0
  disp4=0
  disp5=0
    
; Restore reconstruction matrices
  RESTORE, filename='/kroot/rel/ao/qfix/data/IdlParms/zern_rec-keck_55.idl'

  Nx_ph = N_ELEMENTS(phasemask[*,0])

  phasemap = WHERE(phasemask)

  N_phases = N_ELEMENTS(phasemap)
  N_modes = N_ELEMENTS(zernrec[*,0])
  N_terms = N_ELEMENTS(zernrec[0,*])
  N_frames = N_ELEMENTS(phase[0,*])
  
  C_z_ind=phase##TRANSPOSE(zernrec##vmodes_ph)
  sig_z=TOTAL(C_z_ind^2,2)/N_frames 
  
  var_kol = FLTARR(N_terms)
; # of radial orders
  N_orders = FIX(SQRT(0.25 + 2 * N_terms) - 0.5 )
  
  r0 = FLTARR(N_orders)
  dr0 = FLTARR(N_orders)
  var_data = FLTARR(N_orders)

; New way of calculating Kolmogorov variances
  j = 1
  for n = 1, n_orders-1 do begin
     var_kol_order = gamma(14./3.) * gamma(n - 5./6.) /  $
                     ( 2.^(14./3.) * (gamma(17./6.))^2. * gamma(n + 23./6.) )
     var_kol_order = var_kol_order * 0.046/!pi * (n+1.) / 2.^(5./3.)
     var_kol_order = var_kol_order * (2. * !pi) ^ (11./3.)
                                ; this is directly from Noll 76

     for m = -n, n, 2 do begin
        var_data[n] = var_data[n] + sig_z[j]
        var_kol[j] = var_kol_order
        j = j + 1
     endfor
     var_data[n] = var_data[n] / (n+1.)

     dvar = STDDEV(sig_z[j-n-1:j-1]) ; variance of sigma^2 in
                                ; this radial order
     r0[n] = eff_tel_dia * (var_kol_order/var_data[n])^(3./5.)

     dr0[n] = 3./5. * eff_tel_dia * var_kol_order^(3./5.) * var_data[n]^(-8./5.) * dvar

  endfor
  
  print
  print,"Kolmogorov: r0 (different radial orders)"
  print,r0[2:N_orders-3]
;  print,dr0[2:N_orders-3]

; Fit using von Karman Zernike variances
  RESTORE, filename='/kroot/rel/ao/qfix/data/IdlParms/var_vk.idl' 
  N_L0 = N_ELEMENTS(L0s)

  C2_vk = FLTARR(N_L0)
  res_vk = FLTARR(N_L0)
  r0_vk = FLTARR(N_L0)

  for i = 0,N_L0-1 do begin
     var_th = var_vk[1:N_orders-4,i] ; Noll and Winker mode
     var_dt = var_data[2:N_orders-3] ; numbering are different!!
     C2_vk[i] = TOTAL(var_dt*var_th)/TOTAL(var_th*var_th)
     r0_vk[i] = eff_tel_dia / C2_vk[i]^(3./5.)
     res_vk[i] = TOTAL((C2_vk[i] * var_th - var_dt)^2)
  endfor

; Use the minimum_parameter routine that was originally written for
; DIM analysis, therefore somewhat strange parameters
  pars = MINIMUM_PARAMETERS(res_vk, C2_vk, FLTARR(N_L0), L0s)
  r0_vk = pars[0]*eff_tel_dia
  L0_vk = pars[1]

  print
  print,"von Karman r0, L0:",r0_vk,L0_vk
  print

; Set up array for plot van Karman Zernike vars
  best_vk = WHERE(res_vk eq min(res_vk))
  var_vk_modes = FLTARR(N_terms)
  j = 1
  for n = 1, n_orders-1 do begin
     for m = -n, n, 2 do begin
        var_vk_modes[j] = C2_vk[best_vk] * var_vk[n-1,best_vk]
        j = j + 1
     endfor
  endfor	

; Display Zernike mode variances
  if disp3 then begin
     window,2
     plot,sig_z, psym=2, xtitle="Mode Number", $
          ytitle="Variance [radians^2]", $
          title="Variances of Zernike Mode Coefficients", $
          yrange=[0,max(sig_z[3:5])*1.2]
     C_kol = var_kol*(eff_tel_dia/r0[3])^(5./3.)
     oplot,var_kol*(eff_tel_dia/r0[3])^(5./3.)
     oplot,var_vk_modes,linestyle=2
     print, "Variances of Zernike modes [radians^2]"
     print, "Solid line: Kolmogorov spectrum (fit at 4rd radial order)"
     print, "Dashed line: Van Karman spectrum"
     print,"< Press any key to continue >"
     temp = get_kbrd(1)
     print
     wdelete,2
  endif	


; Display van Karman Zernike residuals
  if disp4 then begin
     window,2
     plot,L0s,res_vk,/xlog,/ylog,psym=-2, $
          title = "Residuals of Best Fit to van Karman Zernike variances", $
          xtitle = "L0 [m]", $
          ytitle = "Best Fit Residuals [arbitrary units]"
     print
     print,"Residuals of Best Fit to van Karman Zernike variances"
     print,"< Press any key to continue >"
     temp = get_kbrd(1)
     print
     wdelete,2
  endif	


; Display Kolmogorv Zernike r0s with error bars
  if disp5 then begin
     window,2
     x = findgen(N_orders-4)+2
     maxy = max(r0[2:N_orders-3]+dr0[2:N_orders-3])
     plot,x,r0[2:N_orders-3],psym=-2, $
          xrange=[0,N_orders-3], yrange=[0,maxy], $
          title = "r0s from Kolmogorov Zernike Variances", $
          xtitle = "Zernike radial order", $
          ytitle = "r0 [m]"
     for i = 2, N_orders-3 do begin
        oplot,[i,i],[r0[i]-dr0[i],r0[i]+dr0[i]]
     endfor
     print
     print,"r0s from Kolmogorov Zernikes with error bars"
     print,"< Press any key to continue >"
     temp = get_kbrd(1)
     print
     wdelete,2
  endif	

; Set up and return results structure
  res ={	r0:		FLTARR(N_orders), 	$
                dr0:		FLTARR(N_orders),	$
                vk:		FLTARR(2),		$
                sigz:		FLTARR(N_terms)}

  res.r0 = r0
  res.dr0 = dr0
  res.vk = [r0_vk,L0_vk]
  res.sigz = sig_z

  RETURN, res
END

PRO ESTIMATESEEING,timestampvec,length,condition=condition,r0vk=r0vk,L0vk=L0vk,fwhm=fwhm,dmloop=dmloop
  
  IF NOT KEYWORD_SET(condition) THEN condition='<='  
  
  nQueries=N_ELEMENTS(timestampvec)
  
  FOR c1=0,nQueries-1 DO BEGIN
     timestamp=timestampvec[c1]
     
     ; error handling on the trsquery call
     CATCH, errorVar
     IF errorVar NE 0 THEN BEGIN
        CATCH, /cancel
        MESSAGE,/INFO,'Could not get data from TRS'
        MESSAGE,/INFO,!error_state.msg
        errormsg={errormsg:!error_state.msg}
        dmloop = 0 ;; added by DLM on 07/28/08 nightlog ticket #k2-15753
        RETURN
     ENDIF
     
; figure out how many samples we need  
     trsdata=trsquery(timestamp,2,'ffb','timestamp',condition=condition)
     timediff=FLOAT(trsdata.timestamp[1]-trsdata.timestamp[0])/1e7
     samples=FIX(length/timediff<32767)
          
     data=trsquery(timestamp,samples,'ffb','dmcommand','timestamp',condition=condition)
     
     CATCH, /cancel
     
     dm0=readdm('flatfile.dm')
     data.dmcommand=data.dmcommand-dm0#replicate(1.,samples)
     
; subtract the flat file, since if the DM is equal to the flatfile, it's flat  
     IF MIN(TOTAL(data.dmcommand^2,1)) EQ 0 THEN BEGIN
        dmloop=0
        MESSAGE,/INFO,'DM loop was open'
     ENDIF ELSE BEGIN
        dmloop=1
        RESTORE, filename='/kroot/rel/ao/qfix/data/IdlParms/zern_rec-keck_55.idl'
        phasemaskext=FLTARR(21,21) 
        phasemaskext[2:18,2:18]=phasemask
        wpme=WHERE(phasemaskext)
        z1=ZERNIKE(21,21,1)
        wz=WHERE(z1)
        z1[wz]=INDGEN(349)+1
        tmp=z1*phasemaskext
        wtmp=WHERE(tmp)
        foo=z1[wtmp-1]
        
        phv=data.dmcommand[foo,*]
        
        phase=phv*0.6*(2.*!pi)/0.5 ; 0.5 microns, 0.6V/micron
        phase=phase-TOTAL(phase,2)/samples#replicate(1.,samples)
        
        r0_zernike = CALC_R0_ZERNIKE(phase)
        
        r0vk=r0_zernike.vk[0]
        L0vk=r0_zernike.vk[1]
        
        fwhm=1.03*500e-9/r0vk*206265.
        fwhm=fwhm*SQRT(1-2.813*(r0vk/L0vk)^0.356); Tokovinin, From DIM to seeing
        
        PRINT, 'Time '+timestamp
        PRINT, 'Seeing at 500 nm ='+STRING(fwhm,format='(f5.2)')+' arcsec'
     ENDELSE
  ENDFOR
END