Skip to content
shark_simulation.pro 12.9 KiB
Newer Older
Marco Stangalini's avatar
Marco Stangalini committed
;+
; SHARH@LBT simulation tool 
; 
; This code simulates the coronographic psf of LBT/SHARK in real seeing conditions.
; A param file in the same working directory is used to configure the code.
; The code works with either simulated or measured AO phase residuals
; See param file for more information on the setup of the code and a complete list of input parameters
; 
; HISTORY:
; Sept 2013 added configuration file easy setup (param)
; Sept 2013 telescope jitter simulation added 
; Oct 2013 FLAO mode: measured AO residuals (good up to the R=200 mas)
; Nov 2013 LBT interferometric mode
; Nov 2013 
;          -- LBTI Code by F. Pedichini implemented and integrated 
;          -- New configuration file added
;          -- Sim report file added 
;          -- Bugs fixed 
; Dec 2013 
;          -- A nasty bug affecting the pixel scale in LBTI mode fixed
;          -- pix variable introduced. It contains the pixel scale in mas as obtained from the prop lib
;             after the propagation of the beam using the function prop_get_sampling()   
;          -- brief and incomplete description in the save file added
;          -- LBT Focal plane PSF estimation added 
;           
; Feb 2014
;          -- Added possibility to assess occulter off axis displacement. This is done in both single and interferometric mode. 
;          -- Report filename fixed 
;          -- 
;          
; Sept 2017
;          -- SOUL LBT AO Upgrade added
;          
; Jul 2018
;         -- Long time series phase screens capabilities added 
; 
; Author: Marco Stangalini marco.stangalini@inaf.it
;First release August 2013
; May 2017 Bug fixed - Working version
;-


@param

label='SHARK_'
HWRS=0

PRINT, '--------------------------------'
PRINT, 'SHARK SIMULATION RUN...
PRINT, '--------------------------------'
PRINT, SYSTIME()
IF FLAO eq 1 THEN begin
  PRINT, 'FLAO ON SKY MODE '
  label=label+'FLAO_'

ENDIF
PRINT, 'LAMBDA [MU]', LAMBDA0
PRINT, 'LYOT STOP [%]', RAGGIO_LYOT
PRINT, 'OCCULTER [MU]', RAGGIO_OCCULTER
PRINT, 'OCC. OFF AXIS DISPLACEMENT [mm]', XOFF, YOFF 
IF SWING EQ 1 THEN begin
  PRINT, 'TELESCOPE SWING YES'
ENDIF
PRINT, '# STEP', NSTEP
PRINT, 'PHASE RESIDUALS', PATH0
print, 'PUPIL SAMPLING [m/pixel]', 8.4/(beam_ratio*lsize)
PRINT, 'BEAM RATIO ', BEAM_RATIO

IF interfmode eq 1 then begin 
label=label+'LBTI_'
print, 'LBT INTERFEROMETRIC MODE ACTIVE'
print, 'PUPIL SAMPLING [m/pixel]', Dtot/(beam_ratio*lsize)
print, 'Effective beam ratio = ', (8.4/Dtot)*beam_ratio
s_sampling=Dtot/(beam_ratio*lsize)
npix_pupil=round(8.4d/s_sampling)
if (2*(npix_pupil/2) ne npix_pupil) then begin 
  npix_pupil=npix_pupil+1
endif
print, 'NPIX PUPIL', NPIX_PUPIL, '(it must be close to 150)'
print, 'WIND SPEED ', WIND_SPEED
PRINT, 'AO LOOP TIME ', phase_sampling, '[ms]
delta_step=round((Cent2Cent_SEP/wind_speed)/(PHASE_SAMPLING))
PRINT, 'WFs step difference ', delta_Step
IF tot_screens lt (step0+nstep+delta_step) then begin
  print, '-- WARNING -- number of steps exceeds the number of phase screens available'
  print, 'Increase the number of phase screens or reduce the wind speed'
  stop
ENDIF
endif
PRINT, '--------------------------------'
PRINT, '--------------------------------'


openw, 33, path+'Report.dat' 
PRINTf,33, '--------------------------------'
PRINTf,33, 'SHARK SIMULATION RUN...
PRINTf,33, '--------------------------------'
PRINTf,33, SYSTIME()
PRINTF,33,'FLAO', flao
PRINTf,33, 'LAMBDA [MU]', LAMBDA0
PRINTf,33,'LYOT STOP [%]', RAGGIO_LYOT
PRINTf,33, 'OCCULTER [MU]', RAGGIO_OCCULTER
PRINTF,33, 'OCC. OFF AXIS DISPLACEMENT [mm]', XOFF, YOFF
PRINTF,33, 'LYOT OFF AXIS DISPLACEMENT [mm]', LYOT_XOFF, LYOT_YOFF
printf,33, 'SWING', swing
PRINTf,33, '# STEP', NSTEP
PRINTf,33, 'PHASE RESIDUALS', PATH0
printf,33, 'PUPIL SAMPLING [m/pixel]', 8.4/(beam_ratio*lsize)
PRINTf,33, 'BEAM RATIO ', BEAM_RATIO
  printf, 33,'LBT INTERFEROMETRIC MODE', interfmode
  printf,33, 'PUPIL SAMPLING [m/pixel]', Dtot/(beam_ratio*lsize)
  printf,33, 'Effective beam ratio = ', (8.4/Dtot)*beam_ratio
  s_sampling=Dtot/(beam_ratio*lsize)
  printf, 33,'WIND SPEED ', WIND_SPEED
  PRINTf, 33,'AO LOOP TIME ', phase_sampling, '[ms]
  delta_step=round((Cent2Cent_SEP/wind_speed)/(PHASE_SAMPLING))
  PRINTf,33, 'WFs step difference ', delta_Step
PRINTf, 33,'--------------------------------'
PRINTf,33, '--------------------------------'
close, 33



;SWING FREQUENCY
AMPLITUDE_TILT_CALIBRATION=0.045e-6
t=FINDGEN(TOT_SCREENS)*PHASE_SAMPLING
MODULATION=AMPLITUDE_TILT_CALIBRATION*(SWING_AMPL/4.12)*SIN(2.*!PI*SWING_FREQ*T)
pupil_distance=Cent2Cent_SEP/2.
pix=0.
XOFF=XOFF/1000.   ;CONVERSIONE IN METRI
YOFF=YOFF/1000.   ;CONVERSIONE IN METRI
LYOT_XOFF=LYOT_XOFF/1000.   ;CONVERSIONE IN METRI
LYOT_YOFF=LYOT_YOFF/1000.   ;CONVERSIONE IN METRI


shortcount=0

nframes=fix(inttime/(PHASE_SAMPLING*1000))
nframe_count=0
psf_short=fltarr(lsize,lsize)
psfao_short=fltarr(lsize,lsize)

;prop_run, 'psd_error_map', psf0, lambda0(0), lsize, dx, passvalue=passval_speckle
;print, 'Error map created'
if SOUL eq 1 then begin
  WFS_DATA=readfits(path0)
  maschera=shift(dist(220,220), 110,110)
  maschera=maschera lt 107
endif

if FLAO eq 1 then begin
  restore, path0, /v
endif

if SYNTH_FFT eq 1 then begin
  paths=file_search(path0+'*.fits')
  WFS_DATA=readfits(paths(0))
  maschera=shift(dist(220,220), 110,110)
  maschera=maschera lt 107
endif

  

FOR OFF=0,N_ELEMENTS(XOFF)-1 DO BEGIN
      for j=0,n_elements(raggio_lyot)-1 do begin
          for k=0,N_elements(raggio_occulter)-1 do begin
            FOR LAM=0,N_ELEMENTS(LAMBDA0)-1 DO BEGIN   
              lyot_r=raggio_lyot(j)
              occrad=raggio_occulter(k)
              print, lyot_r, ' %', occrad, ' lambda/D', lambda0, ' nm'
              WFLambda=750. ;nm
              WFLambda=wflambda *1e-9 ; in meters
              occrad=occrad*1.e-6
              
                WF_res=fltarr(220,200)
                WF_turb=fltarr(220,200)

              psf_shark=fltarr(lsize,lsize)
              psf_AO=fltarr(lsize,lsize)
              psf_fp=fltarr(lsize,lsize)
              
              
              count=0
              counter=0
              device, decomposed=0
              loadct, 5
              psf_shark_8th=fltarr(lsize,lsize)
              psf_shark_gauss=fltarr(lsize,lsize)
              
              
              cnt_synth_blocknum=0
              cnt_synth=0
              
              
              ;Time LOOP =======================================================
              for index = step0, step0+nstep-1 do begin
                  
                    lambda=lambda0(LAM)*1000*1e-9 ;in meters

                    if SOUL eq 1 then begin
                      wf_res = maschera*WFS_data[index,*,*]*1e-9  ;in Dlambda [nm]
                    endif else begin
                      
                      cnt_synth_blocknum++
                      if cnt_synth_blocknum eq 5000 then begin
                        cnt_synth++
                        WFS_DATA=readfits(paths(cnt_synth))
                        cnt_synth_blocknum=1
                      endif
                      
                      wf_res = maschera*WFS_data[*,*,cnt_synth_blocknum-1]*1e-9  ;in Dlambda [nm]
                      
                    endelse
                    
                    

                    phase_error_res=reform(wf_Res);/(2.*!pi)
                    temp=dblarr(lsize, lsize)
                    phase_error_res=congrid(phase_error_res, lsize*beam_ratio, lsize*beam_ratio, cubic=-0.5)
                    size_temp=lsize*beam_ratio
                    temp[lsize/2-size_temp/2:lsize/2+size_temp/2-1, lsize/2-size_temp/2:lsize/2+size_temp/2-1]=phase_error_res
                    phase_error_res=temp  
                  

                  	print,index-step0,  lyot_r, ' %', occrad, ' lambda/D', lambda0(LAM), ' nm'
                  	passval={beam_ratio:beam_ratio, fase:phase_error_res,   occulter:occrad, lyotstop:lyot_R, OSCILLATION:SWING, SWING_a:MODULATION(INDEX),$
                  	  pupil_sep:pupil_distance, D0:Dtot, sampl:pix, sampl_fp:pix, SUPPRESSOR:raggio_suppressor, xoff:xoff(OFF), yoff:yoff(OFF),$
                  	   LYOT_XOFF:LYOT_XOFF, LYOT_YOFF:LYOT_YOFF, Zer_val:Zer_val, occtype:occtype,$
                  	   HWRS:HWRS, zer_num:zer_num}
              
                    ;prop_run,'psf_lbt_single', psffp, lambda0(LAM), lsize, dx_fp, passvalue=passval,/quiet
                  	prop_run,'pisces_on_axe_phased', psf0, lambda0(LAM), lsize, dx, passvalue=passval,/quiet
                    prop_run,'pisces_on_axe_mask_2', psf_cor, lambda0(LAM), lsize, dx, passvalue=passval,/quiet
                    
                    
                    tvframe, psf0^0.3, /asp                    
 

              
                  	M=MAX(PSF0, MA)
                  	psf_shark=psf_shark+psf_cor
                  	psf_AO=psf_AO+psf0
                  	count=count+1
                  	counter=counter+1.
                
                
                  	mas=indgen(lsize)*passval.sampl
                  	nframe_count=1+nframe_count
                  	psf_short=psf_short+psf_cor
                  	psfao_short=psfao_short+psf0

                  	
                  	if nframe_count eq nframes then begin
                  	  IF SHORT_EXP EQ 1 THEN BEGIN

                  	    ;tvframe, psfao_short^0.5
                  	    

                  	    if occtype eq 'gauss' or occtype eq 'Gauss' then begin
                  	      save, filename=path+SUBDIR+label+strtrim(fix(lambda0(0)*1000),1)+'nm_'+'_NSTEP_'+strtrim(shortcount,1)+$
                  	        '_Gauss_Occulter_'+strtrim(occrad,1)+'_OFF_AXIS_DISPL_'+strtrim(XOFF(OFF),1)+'_STOP_'+strtrim(lyot_r,1)+'_STOP_DISPL_'+strtrim(lyot_XOFF,1)+'.sav', $
                  	        psf_short, psfAO_short, LSIZE, DX, mas,$
                  	        description='mas = array [marcsec] -- dx=pixel scal in m -- pix=pixel scale in mas, occrad=occulter R in mu'
                  	    endif else begin
                  	      save, filename=path+SUBDIR+label+strtrim(fix(lambda0(0)*1000),1)+'nm_'+'_NSTEP_'+strtrim(shortcount,1)+$
                  	        '_Lyot_Occulter_'+strtrim(occrad,1)+'_OFF_AXIS_DISPL_'+strtrim(XOFF(OFF),1)+'_STOP_'+strtrim(lyot_r,1)+'_STOP_DISPL_'+strtrim(lyot_XOFF,1)+'.sav', $
                  	        psf_short, psfAO_short, LSIZE, DX, mas,$
                  	        description='mas = array [marcsec] -- dx=pixel scal in m -- pix=pixel scale in mas, occrad=occulter R in mu'
                  	    endelse
                        shortcount++
                        psf_short=psf_short*0.
                        psfao_short=psfao_short*0.


                  	    nframe_count=0
                  	  ENDIF
                  	endif                
                
                
                endfor
                


                print,'pixel=',passval.sampl,' marcsec'
                mas=indgen(lsize)*passval.sampl
                mas_fp=indgen(lsize)*passval.sampl_fp
                
              device, decomposed=0
              psf_AO=psf_AO/(counter)
              psf_shark_lyot=psf_shark/counter
              
              cartop, psf_AO, lsize/2, lsize/2,lsize/2,findgen(lsize/2), 2.*!pi*findgen(360)/360,frt
              siz=size(frt)    
              psf_AO_prof=total(frt,2)/siz(2)
              
              cartop, psf_SHARK_lyot, lsize/2, lsize/2,lsize/2,findgen(lsize/2), 2.*!pi*findgen(360)/360,frt
              siz=size(frt)
              psf_SHARK_LYOT_prof=total(frt,2)/siz(2)
          
           pix=passval.sampl
           pix_fp=passval.sampl_fp
           
           
          PSF_AO=psf_AO
          PSF_CORO=psf_SHARK_LYOT
          PSF_CORO_PROF=PSF_SHARK_LYOT_prof
          PSF_AO_PROF=psf_ao_prof
          
          plot, psf_ao_prof/max(psf_ao_prof), /ylog, yr=[1e-7,1]
          oplot, psf_coro_prof/max(psf_ao_prof), thick=4
          
          
         
          
          if occtype eq 'gauss' or occtype eq 'Gauss' then begin
            save, filename=path+label+strtrim(fix(lambda0(LAM)*1000),1)+'nm_'+'_NSTEP_'+strtrim(nstep,1)+$
              '_Gauss_Occulter_'+strtrim(occrad,1)+'_OFF_AXIS_DISPL_'+strtrim(XOFF(OFF),1)+'_STOP_'+strtrim(lyot_r,1)+'_STOP_DISPL_'+strtrim(lyot_XOFF,1)+'.sav', $
              PSF_AO, nstep, mas, PSF_CORO, dx, PIX,  lambda0,  PSF_CORO_PROF, occrad, LYOT_R, PSF_AO_PROF,$
            description='mas = array [marcsec] -- dx=pixel scal in m -- pix=pixel scale in mas, occrad=occulter R in mu'
          endif else begin
            save, filename=path+label+strtrim(fix(lambda0(LAM)*1000),1)+'nm_'+'_NSTEP_'+strtrim(nstep,1)+$
              '_Lyot_Occulter_'+strtrim(occrad,1)+'_OFF_AXIS_DISPL_'+strtrim(XOFF(OFF),1)+'_STOP_'+strtrim(lyot_r,1)+'_STOP_DISPL_'+strtrim(lyot_XOFF,1)+'.sav', $
              PSF_AO, nstep, mas, PSF_CORO, dx, PIX,  lambda0,  PSF_CORO_PROF, occrad, LYOT_R, PSF_AO_PROF,$
            description='mas = array [marcsec] -- dx=pixel scal in m -- pix=pixel scale in mas, occrad=occulter R in mu'
          endelse
          
          print, occrad, Lyot_r, lambda
          
          
           
           
            ENDFOR
          endfor
      endfor
ENDFOR

end