; speckle simulator, monochromatic and noiseless, exp.time only
; Oct 24-25, 2023 AT. See  shsim2a.pro as presursor
; Lower screen translates in X, upper screen in Y
; All parameters are in the code
; polychromatic speckle
;-----------------------------------------
pro simspec4, seeing

; All parameters are hard-coded at the moment
ngrid = 100 ; half-size of the image, pixels
Nfr = 400   ; number of frames to simulate
;seeing = 0.8  ; seeing at lambda, arcsec

D = 4.1        ; telescope diameter, m
eps = 0.25     ; central obscuration, relative
pixel = 0.01575 ; pixel, arcsec
lambda = 0.80 ; reference imaging wavelength, micron
wind1 = 8.     ; lower wind wpeed, m/s
wind2 = 40.    ; upper wind, m/s
texp  = 0.030   ; exposure time, second
tstep = 8     ; steps per exposure time
chunk = Nfr/20 ; progress bar display
highfrac = 0.5 ; fraction of turbulence in the high layer 

; spectral response of I, filter+CCD-QE
nlam = 5
wlam = 0.75 + 0.05*findgen(5)
resp = [0.724, 0.694, 0.539, 0.347, 0.190] ; 822nm effective 
resp = resp/total(resp) 

; monochromatic speckle
;nlam = 1
;wlam = lambda
;resp = 1.

lameff = total(wlam*resp)

sampl = lameff*1e-6/D*206265./pixel ;  pixels per lambda_eff/d
fc = 2.*ngrid/sampl ; cutoff frequency, pixels
Rpix = ngrid/sampl ; pupil radius in the grid

r0 = 0.98*lambda*1e-6/(seeing/206265.) ; at imaging wavelength
;r0 = (lambda/0.5)^(1.2)*0.98*lambda*1e-6/(seeing/206265.)

Dr0 = D/r0
print, 'Effective wavelength: ', lameff
print, 'Fried parameter at imaging wavelength, D/r0: ', r0, Dr0
print, 'Cutoff frequency [pix]: ', fc
print, 'Ready for the phase screen simulation'

; Seed for random-number generator
  seed = long(systime(/seconds))

print, 'Generating large phase screens...'

 kphase = 2  ; number of fine pixels per phase-screen pixel
 ngrida = 512 ; half-size of atmospheric grid
 size = 206265.*(1e-6*lambda)/pixel
 spixel = size/(2*ngrid)*kphase ; phase-screen pixel [m]
 phsize = spixel*2*ngrida  ; phase-screen size, m

 r01 = r0*(1. - highfrac)^(-3./5.) 
 r02 = r0*(highfrac)^(-3./5.) 
; check: r0 = (r01^(-5./3.) + r02^(-5./3.))^(-3./5.)

; fact = sqrt(0.023)*(phsize/r0)^(5./6.) 
 fact1 = sqrt(0.023)*(phsize/r01)^(5./6.) 
 fact2 = sqrt(0.023)*(phsize/r02)^(5./6.) 

 wshift1 = wind1*texp/spixel ; wind shift per exposure, pix, low layer
 wshift1 = round(wshift1) 
 wspeed1 = wshift1*spixel/texp
 wshift2 = wind2*texp/spixel ; wind shift per exposure, pix, low layer
 wshift2 = round(wshift2) 
 wspeed2 = wshift2*spixel/texp
 print, ' Low wind speed [m/s] & shift [pix]: ', wspeed1, wshift1
 print, 'High wind speed [m/s] & shift [pix]: ', wspeed2, wshift2

 wshift1a = wshift1*kphase/tstep ; shift in fine pixels
 wshift2a = wshift2*kphase/tstep ; during exposure

 tmp = shift(dist(2*ngrida),ngrida,ngrida) ; radius, large grid

; low phase screen
  tmp1 = fact1*tmp^(-11./6.)*dcomplex(randomn(seed,2*ngrida,2*ngrida,/double),randomn(s,2*ngrida,2*ngrida,/double))  
  tmp1[ngrida,ngrida] = dcomplex(0.,0.)
  phasescreen1 = float(shift(fft(shift(tmp1,ngrida,ngrida),/inverse),ngrida,ngrida))
; high phase screen
  tmp1 = fact2*tmp^(-11./6.)*dcomplex(randomn(seed,2*ngrida,2*ngrida,/double),randomn(s,2*ngrida,2*ngrida,/double))  
  tmp1[ngrida,ngrida] = dcomplex(0.,0.)
  phasescreen2 = float(shift(fft(shift(tmp1,ngrida,ngrida),/inverse),ngrida,ngrida))


  phx1 = 0 ; initialize shift of low screen
  phy1 = 0 
  phx2 = 0 ; initialize shift of high screen
  phy2 = 0 
  slide = fix(512/(ngrid/kphase)) ; number of patches in the large screen

print, 'Screens are ready, starting the loop'
;stop

; prepare the variables
cube = fltarr(2*ngrid,2*ngrid,Nfr)
r = shift(dist(2*ngrid),ngrid,ngrid) ; radial distance from grid center, pixels
r[ngrid,ngrid] = 1.e-3

ur = fltarr(2*ngrid,2*ngrid)    ; zero array
inside = where((r le Rpix) and (r ge eps*Rpix))
ur(inside)=1.    ; input amplitude: real part
u0 = complex(ur, ur*0.) ; unit amplitude inside the pupil

; prepare coordinates for stretching images
xx = findgen(2*ngrid) # replicate(1., 2*ngrid)
yy = transpose(xx)


; ---- Main loop --------
; for shifts during exposure shift the phase-screen fragment

for i=0,Nfr-1 do begin

; carve the pieces of the phase screens and compute atmospheric phase
  phx1++ ; left corner in the large phase screen
  if (phx1 mod slide) eq 0 then phy1++ ; slide in y
;  print, 'i, phase screen shifts: ',i, phx,phy
  tmp1 = (shift(phasescreen1,wshift1*phx1,wshift1*phy1))[0:2*ngrid/kphase-1, 0:2*ngrid/kphase-1] 

  phx2++ ; left corner in the large phase screen
  if (phx2 mod slide) eq 0 then phy2++ ; slide in y
  tmp2 = (shift(phasescreen2,wshift2*phy2,wshift2*phx2))[0:2*ngrid/kphase-1, 0:2*ngrid/kphase-1] 

 aphase1 = rebin(tmp1, 2*ngrid, 2*ngrid) ; re-binned on fine pixels
 aphase2 = rebin(tmp2, 2*ngrid, 2*ngrid) ; re-binned on fine pixels

; utmp = u0 ; un-blurred image
; utmp[inside] = complex(cos(aphase[inside]), sin(aphase[inside])) 
; imh =  abs(shift(fft(shift(uampl,ngrid,ngrid),/inverse),ngrid,ngrid))^2
; stop

 img = fltarr(2*ngrid, 2*ngrid) ; accumulate average image
 for j=0,tstep-1 do begin ; inter-exposure  shifts
   atmp = shift(aphase1, j*wshift1a,0) +  shift(aphase2, 0, j*wshift2a) 
   atmp = atmp[inside] ; phase at lambda
   for k=0,nlam-1 do begin  ; spectral-response loop
     f = lambda/wlam[k]
     uampl = u0 
     uampl[inside] = complex(cos(atmp*f), sin(atmp*f)) 
     imh = abs(shift(fft(shift(uampl,ngrid,ngrid),/inverse),ngrid,ngrid))^2
     imh1 = bilinear(imh, ngrid + (xx - ngrid)*f, ngrid +(yy - ngrid)*f) ; stretch

     img += imh1*(resp[k]/total(imh1))
;     stop
   endfor ; k
 endfor ; j
; print, 'Exposure averaging done'
 tvscl, img
; stop  
 cube[*,*,i] = img
 if ((i mod chunk) eq 0) then print, '*', format='(A1,$)' ; progress bar

endfor ; i

writefits, 'simspec4.fits', cube
print, 'Simulated speckle data saved in simspec4.fits!'
;stop

end