; Extract binary parameters from the speckle power spectra
; Jul 18, 2008. AT
; dir1 ='/uwa1/atokovin/20080714/'
;-------------------------------------------------
;-------------------------------------------------
FUNCTION binfunc, ix, par
; Levenberg-Marquart Gaussian function: value and derivatives
; ix is the index of the point in fitgood array
; parameters = [A,B,dx,dy]
; function F = ref*( A + B*cos(fx*dx + fy*dy))

COMMON fitpar, refsp, fitgood, fx, fy

index = fitgood[ix]
fx1 = fx[index] & fy1 = fy[index]

  n = n_elements(fx[*,0])
  arg = 2.*!pi*(par[2]*fx1 + par[3]*fy1)/n
  c1 = cos(arg)
  c2 = -2.*!pi*sin(arg)/n

  tmp = [par[0]*(1.+ par[1]*c1), 1.+par[1]*c1, par[0]*c1, par[1]*par[0]*c2*fx1, par[1]*par[0]*c2*fy1]

 return, tmp*refsp[index]
end
;-----------------------------------------------------
;--------------------------------------------------------
function radav, ima, r
; Radially average the image around the  center
; returns 2D array 

n = n_elements(ima(*,0))

id = fix(r)
nrad = fltarr(n) ; number of pixels at given distance
frad =  fltarr(n) ; sum of fluxes

for i=0,n-1 do for j=0,n-1 do begin
  ind = id[i,j]
  nrad[ind] +=1
  frad[ind] += ima[i,j] 
endfor

;sel = indgen(n/2)
;frad = frad[sel]/nrad[sel]
;frad[n/2:n-1] = 0.

frad = frad/nrad

res = fltarr(n,n)
for i=0,n-1 do for j=0,n-1 do res[i,j] = frad[id[i,j]] 

;stop

 return, res

end
;-------------------------------------------------------
;-----------------------------------------------------
pro b1, dir, name, par=par, rad=rad, dm0=dm0, self=self

; optional par=[dx,dy] to indicate the components
; rad is the min. radius for ACF peak search, 3 pix by default
; the reference file must be processed first, saved in ref.idl
; set keyword dm to fix the dm 
; set keyword self to use radially-averaged power sp as reference

COMMON fitpar, refsp, fitgood, fx, fy


   fname =  dir+'spectra/powps-'+name+'.fits'

  tmp = file_search(fname, count=cnt)
  if (cnt eq 0) then begin
    print, 'File '+fname+' is not found!'
    return
  endif else begin 
    powsp = readfits(fname, header) 
;  powsp = powsp[100:399,100:399] ; for 400-pixel frames
   n = n_elements(powsp[*,0])
  endelse

; for Blanco
  D = 3.93 ; telescope diameter, m
  pixel = 0.015 ; pixel scale, arcsec
  lambda = 0.5535E-6 ; wavelength, m

  samp = lambda/D*206265./pixel ;  pixels per lambda/d
  fc = n/samp ; cutoff frequency, pixels
  print, 'Sampling, pixels per lambda/D: ', samp

    xbad = [-80,-40,40,80,0,0] ; list of bad pixels, rel. to origin
    ybad = [80, 40, -40, -80,67,-67]
   for i = 0, n_elements(xbad)-1 do powsp[n/2+xbad[i],n/2+ybad[i]] = $
0.5*(powsp[n/2+xbad[i]-1,n/2+ybad[i]] + powsp[n/2+xbad[i]+1,n/2+ybad[i]])

  r = shift(dist(n,n), n/2, n/2)

  out = where(r gt fc)
  nbias = total(powsp[out])/n_elements(out) ; noise offset in the spectrum  
  powsp -= nbias

; ---  determine the initial binary parameters from ACF peak
 pow1 = powsp ; temporary array
 sel = where( (r lt 0.1*fc) and (r gt 0.05*fc))
 const = total(pow1[sel])/n_elements(sel)
 pow1[where(r lt 0.1*fc)] = const

 acf = float(shift(fft(shift(pow1,n/2,n/2),/inverse),n/2,n/2)) 
 writefits, dir+name+'acf.fits', acf

   if (keyword_set(rad)) then minrad=rad else minrad=5.

  acf0 = acf[n/2,n/2]
  tmp = acf
  tmp[where(r lt minrad)] = 0. ; mask below 3-pix radius
  tmp = tmp[*,n/2:n-1] ; upper part

; give optionally input parameters as par=[dx,dy] 
  if  keyword_set(par) then begin 
    dx = par[0]
    dy = par[1]
    acfmax = acf[n/2+dx, n/2+dy]
  endif else begin   
    acfmax = max(tmp)
    i = (where(tmp eq acfmax))[0] 
    dy = i / n
    dx = (i mod n) - n/2
  endelse
 
  peakrat = acfmax/acf0
  print, 'Binary coordinates and peak ratio: ', dx, dy, peakrat
 
  tmp[n/2+ dx, dy] = min(tmp)
;  tmp[n/2+ dx, dy] = -0.2*acf0
  tvscl, congrid(tmp,400,200)

   print, 'Ready for fitting! Type dx=. & dy=. to change, .c to fit'

  ; m=40   
  ; shade_surf, acf[n/2-m:n/2+m-1,n/2-m:n/2+m-1],zr=[0,acf0] ; az=60
  ;  device, fi='acf032rc.ps', bits=8, /color
  
stop

    ; frequency domain for fitting, only upper half 

  fx = (findgen(n) - n/2) # replicate(1, n)
  fy = transpose(fx)
  fmin = 0.1*fc
  fmax = 0.8*fc

;  fmin = 0.01*fc
;  fmax = 0.05*fc

;  fitgood = where( (r gt fmin) and (r lt fmax) and (fy ge 0) )
  fitgood = where( (r gt fmin) and (r lt fmax) and (fy ge 0) and (fx ne 0) )

;  if (keyword_set(self)) then begin
    refsp = radav(powsp, r)
    refname = name
    print, 'Self-reference'
;  endif else  restore, 'ref.idl' ; refsp, refname 

  y = powsp[fitgood]
  m = n_elements(fitgood)
  ix = findgen(m)
 
  y_err = (abs(y) + nbias)/sqrt(200.) ; for 200-frame sequence

  A = total(y)/total(refsp[fitgood])
  B = 2.*peakrat*A
  par = [A,B/A,dx,dy]

  if (keyword_set(dm0)) then begin
    alp = 10.^(-0.4*dm0)
    par[1] = 2.*alp/(1. + alp^2) ;
    fixpar=[1,0,1,1]
  endif else fixpar = [1,1,1,1]


  par0 = par

    p1 = A + B*cos(2.*!pi*(dx*fx + dy*fy)/n)
    p1[where(r lt fmin)] = 0.
  plot, pow1[n/2,*]
  oplot, (p1*refsp)[n/2,*], li=2
  print, 'Fitting the model...'

;tmp = lmfit(ix,y,par,function_name='binfunc',ITER=itnum, sigma=sigma, convergence=convflag, measure_errors=y_err, chisq=chisq)  

  ; patch for dm=0
tmp = lmfit(ix,y,par,function_name='binfunc',ITER=itnum, sigma=sigma, convergence=convflag, measure_errors=y_err, chisq=chisq,fita=fixpar)  

 print, itnum,convflag, ' iterations and convergence'
  print, 'CHISQ/nu= ', chisq/m

  res = fltarr(n,n)
  res[fitgood] = tmp

  print, 'Initial parameters: ', par0
  print, 'Final parameters:   ', par
  print, 'Parameter errors:   ', sigma

  sigx = sqrt(sigma[2]^2 + sigma[3]^2) ; average x,y error
;  p = par[1]/par[0] < 1. ; B/A ratio
  p = par[1] < 1. ; B/A ratio
  rat = (1. - sqrt(1. - p^2))/p ; intensity ratio
  dm = -2.5*alog10(rat) ; magnitude difference


; at Blanco, X points to the West, Y to South if the Luca "tip" points North
  rho = sqrt(par[2]^2 +par[3]^2)*pixel & pa = 180./!pi*atan(par[2], par[3])

; for SOAR, mirror image, swap xy 
;  rho = sqrt(par[2]^2 +par[3]^2)*pixel & pa = 180./!pi*atan(par[3], par[2])
 
 if (pa lt 0.) then pa = pa+360.

  print, 'Rho, theta, dm: ', rho, pa, dm
 

  resid =  (pow1-res)
  resid[where(r lt 0.1*fc)] = 0.
  resid = resid[0:n-1,n/2:n-1]

;  tvscl, res[*,n/2:n-1]
;  tvscl, pow1[*,n/2:n-1], n+10,0

   plot, pow1[n/2+1,*],xr=[n/2,n]
  oplot, res[n/2+1,*], li=2

  tvscl, resid 

; stop
; newref =  powsp/(1. + par[1]*cos(2.*!pi*(dx*fx + dy*fy)/n)) ; synthetic reference
  
;; --------- Save parameters in a log-file bin6.dat
   get_lun, u
   openw, u, 'bin6.dat', /append
   fmt='(3F8.2, F8.4, F8.1, F8.3, F8.1, I3, 2A8 )'
   printf, u, par[2], par[3],sigx, rho, pa, dm, chisq/m, convflag, name, refname, format=fmt
   print, par[2], par[3],sigx, rho, pa, dm, chisq/m, convflag, name, refname, format=fmt
   close, u
   free_lun, u

; stop
;  p1 = par[0] + par[1]*cos(2.*!pi*(par[2]*fx + par[3]*fy)/n)
;  ref2 = powsp/p1 ; artificial reference
;  acf1 = float(shift(fft(shift(ref2,n/2,n/2),/inverse),n/2,n/2)) 
;stop

end
;-------------------------------------------------
pro ref, dir, name

; the reference power spectrum is  processed  and saved
; dir ='/uwa0/atokovin/soar-20071024/'
 
  fname =  dir+'spectra/powps-'+name+'.fits'

  tmp = file_search(fname, count=cnt)
  if (cnt eq 0) then begin
    print, 'File '+fname+' is not found!'
    return
  endif else begin 
    powsp = readfits(fname, header) 
  endelse

  n = n_elements(powsp[*,0])

    xbad = [-80,-40,40,80,0,0] ; list of bad pixels, rel. to origin
    ybad = [80, 40, -40, -80,67,-67]
   for i = 0, n_elements(xbad)-1 do powsp[n/2+xbad[i],n/2+ybad[i]] = $
0.5*(powsp[n/2+xbad[i]-1,n/2+ybad[i]] + powsp[n/2+xbad[i]+1,n/2+ybad[i]])

  D = 4.1 ; telescope diameter, m
  pixel = 0.015 ; pixel scale, arcsec
  lambda = 0.6563E-6 ; wavelength, m

  samp = lambda/D*206265./pixel ;  pixels per lambda/d
  fc = n/samp ; cutoff frequency, pixels
;  print, 'Sampling, pixels per lambda/D: ', samp

  r = shift(dist(n,n), n/2, n/2)

  out = where(r gt fc)
  nbias = total(powsp[out])/n_elements(out) ; noise offset in the spectrum  
  powsp -= nbias

  refsp = powsp
  refname = name
  save, filename='ref.idl', refsp, refname
    print, 'The reference power spectrum is saved!' 

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

;-------------------------------------------------
;-------------------------------------------------
;-------------------------------------------------