; Analyze the speckle series
; Jul 15, 2008 AT
; dir ='/uwa1/atokovin/20080714/'
; fname='tst_20071019.009'
;------------------------------------------
function quadratic, x, a
; z = a[0] + a[1]*x^2 
t1 = x^2
f = a[0]+ a[1]*t1
return, [f, 1.+x*0., t1] ; function and derivatives
end
;------------------------------------------
; input parameters: 
; data directory
; file name (without .fits)
; /BIAS if bias frame is available
;------------------------------------------
pro sp1, dir, fname, BIAS=bias

maskrad = 00 ; radius in pixels
threshold = 0  ; threshold, AU
fixbias = 505L   ; fixed bias, ADU


BiasFlag = 0
if keyword_set(BIAS) then begin 
  restore, bias+'.idl' ; bias array
  BiasFlag = 1
  print, 'Bias frame is ', biasfile
    endif else print, 'Subtracting fixed bias 500ADU'

  img = readfits(dir+fname+'.fits', header) 

  n = n_elements(img[*,0,0])
  nf = n_elements(img[0,0,*])
  print, nf, ' frames in this series'

; ---- Get some parameters from the header --------------------------
   emgain = -1
   exptime = 0.
   emgain = sxpar(header, 'EMGAIN')
   exptime  = sxpar(header, 'EXPTIME')

  if (emgain eq 0) then threshold=0

  imav = float(total(img, 3))/nf
  if (BiasFlag) then imav -= bias else imav -= fixbias

;  tfit  = gauss2dfit(imav,par, /tilt) ; fit Gaussian to the average image
;  angle = par[6]/!pi*180. ; rotation angle in degrees
;  ell =  abs( (par[2]-par[3])/(par[2]+par[3]))
;  fwhm = 2.35*(par[2] + par[3])/2. ; FWHM in pixels
;  fwhm = 2.35*sig/samp ; PSF Gaussian FWHM in lambda/D units
;  offset = par[0] ; constant offset

;  print, 'FWHM, Ell, PA: ', fwhm, ell, angle
;  print, 'Offset [ADU]: ',offset

;  tvscl, imav-tfit
  tvscl, imav
;  stop

  impar = fltarr(3,nf-1) ; xc,yc, flux
  powsp = fltarr(n,n)

  r = shift(dist(n,n), n/2, n/2)
  xx = (findgen(n) - n/2) # replicate(1.,n)
  yy = transpose(xx)
  nempty = 0

  ; prepare for masking the images
 if (maskrad gt 0) then begin
   mask = fltarr(n,n)
   mask[where(r lt maskrad)] = 1.
endif

  print, 'Computing the power spectra...'
  for i=1,nf-1 do begin 
    if (img[0,0,i] gt 20000L) then begin ; reject empty frames
       nf -=1
       print, 'Image number ',i,' empty'
       nempty +=1  
       continue
    endif
    tmp = float(shift(img[*,*,i],1,0)) ; select frame, remove the shift

   if (BiasFlag) then tmp -=bias else tmp -= fixbias ; remove background

; calculate flux and centroids
    flux = total(tmp)
    xc = total(tmp*xx)/flux
    yc = total(tmp*yy)/flux
    impar[*,i-1] = [xc,yc,flux]

if maskrad gt 0. then   tmp = shift(tmp, -fix(xc), -fix(yc))*mask ; re-center and mask
; if threshold gt 0. then   tmp = (tmp-threshold) > 0 ; subtract threshold
if threshold gt 0. then   tmp[where(tmp lt threshold)]=0. ; threshold

    powsp += (abs(fft(tmp)))^2

  endfor
  avflux = total(impar[2,*])/nf-1

  powsp = shift(powsp, n/2,n/2)
  powsp = powsp/powsp[n/2,n/2]   ; normalize

; shade_surf, alog10(powsp)

  tvscl, congrid(-alog10(powsp),400,400)
  print, 'Power spectrum is computed!'

  sxaddpar, header, 'NAXIS', 2 ; correct NAXIS keyword
  sxaddpar, header, 'BITPIX', -32 ; correct BITPIX keyword
  sxdelpar, header, 'NAXIS3' ; remove NAXIS keyword
 
   writefits, dir+'ps-'+fname+'.fits', powsp, header
   writefits, dir+'av-'+fname+'.fits', imav, header
  
stop


  acf = float(shift(fft(shift(powsp,n/2,n/2),/inverse),n/2,n/2))

  peakrad = 3. ; expected radius of speckle peak, pixels

  out = where((r gt peakrad) and (r lt 3.*peakrad)) ; ring around center
  in =  where(r le peakrad)

    arg = r[out] ; argument
    y = acf[out] ; function
    a = [max(y),0.] ; initial parameters
    yfit = lmfit(arg,y,a, function_name='quadratic')     

      acf0 = fltarr(n,n)
      acf0[in] = acf[in] - a[0] - a[1]*r[in]^2 ; subtract the base

      ctr0 = acf0[n/2,n/2]/a[0] ; raw contrast
      print, 'Raw speckle contrast: ', ctr0

     tmp = acf0[n/2-peakrad:n/2+peakrad,n/2-peakrad:n/2+peakrad] 

     tfit  = gauss2dfit(tmp,par) ; fit Gaussian to the average image
     fwhm0 = 2.35*(par[2] + par[3])/2. ; FWHM in pixels
     print, 'Peak FWHM [pix]: ', fwhm0

;  x = findgen(10)
;  plot, acf[n/2:n-1,n/2], xr=[0,10]
;  oplot, a[0]+a[1]*x^2, li=2
       shade_surf, acf

; writefits, 'tmp.fits', acf  
  stop

end
;--------------------------------------------------------
pro getbias, dir, biasfile, biasname

img = readfits(dir+biasfile+'.fits', biasheader) 

n = n_elements(img[0,0,*])-1
nx = n_elements(img[*,0,0])

img = float(img[*,*,0:n-1]) ; throw away the 1st image

;tmp = total(img,1) ; average in x 
;tmp = total(tmp,2) ; average in frame
;tmp = tmp/(n*nx)
;bias = tmp ## replicate(1., nx)
;
bias = total(img,3)/n
print, biasname+' bias frame is generated'

tvscl, bias

save, filename=dir+biasname+'.idl', bias, biasfile

; stop

end
;--------------------------------------------------------
;--------------------------------------------------------
;--------------------------------------------------------
;--------------------------------------------------------
;--------------------------------------------------------