diff --git a/shark_simulation.pro b/shark_simulation.pro
new file mode 100644
index 0000000000000000000000000000000000000000..42126868fe066abb158b4a512b2ac299eb5890f1
--- /dev/null
+++ b/shark_simulation.pro
@@ -0,0 +1,337 @@
+;+
+; 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