;phase_screen_power_spectrum ; ;PURPOSE: ; This script will take in phase screens generated by turb2d.c (written ; by Garret Jernigan) and compute their power spectrum ;CALLING SEQUENCE: ; phase_screen_ellipticity_filter,[date,[FrameDrift,[depth]]]... ; ;KEYWORDS: ; Windspeed: ; The speed at which the frozen screen will be dragged across ; the screen. The units ; PhasDim: ; The dimension of the phase screen (default is 1024) ; AvgSub: ; 0 = don't subtract any offset ; 1 = subtract the average offset of a given frame ; 2 = set the kx=0,ky=0 componenet of FFT to zero ; 3 = subtract a scaled offset from the distorted grid ; FittedOnly: ; 1 = use only the fitted points on each frame ; 2 = use only the non-fitted points one each frame ; 3 = use only points that have non-zero value ; UseFid: ; 1 = apply the fiducial cut ; 0 = don't apply the fiducial cut ; Covariance: ; 0 = don't do covariance analysis ; 1 = do covariance analaysis ; WindSpeed: ; 1 = Use multiple wind speeds ; 2 = Use only one wind speed pro phase_screen_power_spectrum,$ xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$ PhaseDim=PhaseDim,$ low=low,lcut=lcut,$ MultiWindSpeeds=MultiWindSpeeds, AvgSub=AvgSub,$ UseFid=UseFid,Covariance=Covariance, TvFlag=TvFlag If (not keyword_set(AvgSub)) then AvgSub =1 If (not keyword_set(FittedOnly)) then FittedOnly =1 If (not keyword_set(UseFid)) then UseFid =0 If (not keyword_set(Covariance)) then Covariance = 1 If (not keyword_set(PhaseDim)) then PhaseDim=long(1024) If (not keyword_set(low)) then low=1000.0 If (not keyword_set(lcut)) then lcut=100000.0 If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=2 If (not keyword_set(x_off)) then x_off=0 If (not keyword_set(y_off)) then y_off=0 If (not keyword_set(ExpTime)) then ExpTime=.03 If (not keyword_set(TvFlag)) then TvFlag = 0 ;*************************************************************************** ;Calculate The Inner and Outer scales of the turbulence lowstr=strtrim(string(low,format='(%"%10.1f")'),2) lcutstr=strtrim(string(lcut,format='(%"%10.1f")'),2) LowInMeters=PhaseDim*.17/low HighInMeters=PhaseDim*.17/lcut ;FIDUCIAL VALUES nonfiducial=return_fiducial() complement,findgen(626),nonfiducial,fiducial ;************************************************************************* ;FILES AtmosphereDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' PostscriptDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/phase_screen_figs/' ScreenDirectory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' PhaseFiles=file_search(ScreenDirectory+'phase_'+strtrim(lowstr,2)+'_*.out',$ count=NumPhaseFiles) oldway=0 If oldway eq 1 then begin low=1.7 lcut=100000.0 lowstr=strtrim(string(low,format='(%"%10.1f")'),2) lcutstr=strtrim(string(lcut,format='(%"%10.1f")'),2) PhaseFiles=file_search(AtmosphereDir+'phase_'+strtrim(lowstr,2)+'_'+$ strtrim(lcutstr,2)+'*.out',$ count=NumPhaseFiles) endif KMax=PhaseDim/2 Num=PhaseDim xwinds=[0.] ywinds=[0.] CN2s=[1] NumLayers=N_elements(CN2s) TotalPhaseScreens=NumPhaseFiles/NumLayers PhiPowerArr=Fltarr(PhaseDim,PhaseDim,TotalPhaseScreens) ;*************************************************************************** for GlobalScreenNum=0, TotalPhaseScreens-1 do begin ;Final Phase at pupil Screens=fltarr(PhaseDim,PhaseDim,NumLayers) ;Array to Hold Screens for LayerScreenNum=0, NumLayers -1 do begin ;Ind. Phase Layers ScreenNum=GlobalScreenNum*NumLayers+LayerScreenNum ;Open the phase screen binary file and read floating point (data_type=4) openr,Phase,PhaseFiles(ScreenNum),/get_lun PhaseScreenSoar=read_binary(phase,$ data_dims=[PhaseDim,PhaseDim],data_type=4) close,phase free_lun,phase endfor PhiPower=return_entirepowerspectrum(PhaseScreenSoar,FittedOnly,Num) PhiPowerArr(*,*,GlobalScreenNum)=PhiPower print,"Screen Number", GlobalScreenNum endfor TotalPhasePower=avg(PhiPowerArr,2) stop end