; simulate noisy data from noiseless speckle cube
;------------------------------------------------------------------
function otf0, arg
 return, 2./!pi*(acos(arg) - arg*(1.- arg^2)^(0.5))
end
;--------------------------------------------------------
function radav, ima, r
; Radially average the image around the  center
; returns 1D array of the size n/2

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]

  return, frad
end
;------------------------------------------------------------------
;--------------------------------------------------------
; General simulator of detector noise in speckle
; INPUT PARAMETERS:
; fname = fits file with simulated noiseless data cube
; fc - cutoff frequency in pixels
; Nph - number of photoms per frame
; ron - readout noise in ADU
; ampl - amplitude of single-photon event in ADU (set to 1 if CMOS)
; CIC - charge-injection count, =0 for CMOS
; threshold - threshold in ADU
; binary = [dx, dy, alpha] optional binary with secondary ratio alpha
; OUTPUT parameters outpar:
; [Nph, Nph2, pbias, pspar[0], pspar[1], dr0]
; ---
pro noisesimul, fname, fc, nph, ron, ampl, cic, threshold, outpar, binary=binary, $
    powsp=powsp, pow1d=pow1d, ncube=ncube


;fname ='simspec.fits' ; data cube
;Nph = 200     ; photons per frame
;CIC = 0.02 ; iXon rate
;fc = 86.  ; cutoff frequency in pixels
;ampl = 10.  ; amplitude of single-el. events in ADU, set zero for CMOS
;threshold = 10.  ; threshold. ADU
;maskrad  = 50. ; mask radius, pixels

if (ampl le 1.) then ampl = 1. 

cube = readfits(fname)

s = size(cube)
nx = s[1]
nz = s[3]
r = shift(dist(nx,nx), nx/2, nx/2)
; out = where(r gt maskrad) ; for masking

ncube = fltarr(nx,nx,nz) ; noise cube
print, 'Adding photons...'

for i=0,nz-1 do begin 
  tmp = cube[*,*,i]
  if keyword_set(binary) then tmp += binary[2]*shift(tmp, binary[0], binary[1])
  tmp = tmp/total(tmp)*Nph ; normalize by average photon number

; Poisson distribution of events in each pixel
  for k=0,nx-1 do for l=0,nx-1 do begin 
    a =  randomn(seed, poisson=tmp[k,l]+CIC)  ; number of events in each pixel
    if (ampl gt 1) then begin ; EM cCD
      s = 0
      if (a gt 0) then for j=0,a-1 do s += -alog(randomu(seed))
      ncube[k,l,i] = ampl*s + ron*randomn(seed)
    endif else ncube[k,l,i] = a + ron*randomn(seed)  ; CMOS
  endfor  ; l
endfor  ; data cube loop

; clipped cube
ccube = fltarr(nx,nx,nz)
flux2 = 0.
avimg = fltarr(nx,nx) ; accumulate clipped image for background calculation
; powsp1 = fltarr(nx,nx) ; uncorrected power spectrum

for i=0,nz-1 do begin 
   tmp = ncube[*,*,i]  ; one frame

   low = where(tmp lt threshold, nlow) ; thresholding as in getpower.pro
   if nlow gt 0 then tmp[low] = 0. 
   ccube[*,*,i] = tmp
   flux2 += total(tmp)
   avimg += tmp  ; average image
;   powsp1 += abs(fft(tmp)^2)
endfor

back = median(avimg[where(r gt (nx/2 - 2))])/nz ; average background in the corners of the image, in each frame

; Compute power spectrum and modified flux
powsp = fltarr(nx,nx)


; Power-spectrum calculation with subtracted background
for i=0,nz-1 do begin 
  tmp = ccube[*,*,i] - back
  powsp += abs(fft(tmp)^2) ; 
endfor ; powsp loop

; for a constant image=C, powsp(0,0) = C^2*Nz

npheff = total(ncube)/nz/ampl ; before thresholding
npheff2 = flux2/nz/ampl       ; after thresholding

print, 'DONE, eff. Nph, Nph2: ',npheff, npheff2
; writefits, 'ncube.fits', ncube

powsp = shift(powsp, nx/2, nx/2)

;powsp1 = shift(powsp1, nx/2, nx/2)
; power at f=0 extrapolated to remove bias
;powz = 0.5*(powsp[nx/2-1,nx/2] + powsp[nx/2,nx/2+1]) ; remove peak at zero
; linear
;powz = (powsp[nx/2-1,nx/2] + powsp[nx/2,nx/2+1]) - 0.5*(powsp[nx/2-2,nx/2] + powsp[nx/2,nx/2+2])

;powz1 = (powsp1[nx/2-1,nx/2] + powsp1[nx/2,nx/2+1]) - 0.5*(powsp1[nx/2-2,nx/2] + powsp[nx/2,nx/2+2])

; quadratic - not good
;powz1 = (2./3.)*(powsp1[nx/2-1,nx/2] + powsp1[nx/2,nx/2+1]) - (1./6.)*(powsp1[nx/2-2,nx/2] + powsp[nx/2,nx/2+2])

;print, 'Ready to normalize power spectrum'
;stop

;powsp[nx/2, nx/2] -=  back^2*nz ; correct for background spike at f=0

powz = powsp[nx/2,nx/2] ; value at zero frequency without background
powsp = powsp/powz  

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

pow1d = radav(powsp1,r) ; radial average

if (min(pow1d) le 0) then pow1d[where(pow1d lt 1e-7)] = 1e-7 
fx = findgen(nx/2)/fc ; 1-D frequency normalized by cutoff

; fitting linear model
fmax = 0.8 & fmin = 0.2 ; fitting interval
sel = where( (fx lt fmax) and (fx gt fmin))
x = fx[sel] &  y = alog10(pow1d[sel])
pspar = linfit(x,y)

sel1 = where((x ge 0.4) and (x lt 0.6))
s = mean(y[sel1])

powlev = 10^pspar[0] ; at zero
dr0 = (0.435/powlev)^(0.5) ; power at half-cutoff, otf0=0.391

;powlev = 10^(pspar[0] + 0.5*pspar[1]) ; at half-cutoff
;dr0 = (0.435*0.391/powlev)^(0.5) ; power at half-cutoff, otf0=0.391

print, 'Estimated D/r0=', dr0

powmod1 = 10^(pspar[0] + (pspar[1])*fx)

s1 = pspar[0] + 0.5*pspar[1]

print, 'S, p0, p1 ', s, pspar, format='(A,3F8.3)'

!p.charsize = 1.5

plot, fx, pow1d, /ylog, xr=[0,1], yr=[1e-6,1], xtitle='f/f!Dc!N', ytitle='Power'
oplot, fx, powmod1, li=2
oplot, fx, otf0(fx)*0.435*dr0^(-2), li=1

; stop
;tvscl, -alog10(powsp > 1e-7)

; Output results
outpar = [npheff, npheff2, pbias, pspar[0], pspar[1], dr0] 

;stop
;m = 50
;tvscl, congrid(sqrt((ncube[nx/2-m:nx/2+m,nx/2-m:nx/2+m,0]+RON) > 0), 2*(2*m+1),  2*(2*m+1))


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