;phase_screen_multi_layer ; ;PURPOSE: ; This script will take in phase screens generated by turb2d.c (written ; by Garret Jernigan) and simulate N layers with ; structure constants CN2(n) drifting across a telescope aperture. ; ; Assuming geometrical optics, the ellipticity of a star that would be ; formed by the wavefront is calculated. High and low pass filters are ; applied and the ellipticity of these filtered wavefronts is calculated. ; In addition, if the Covariance keyword is set, the covariance of each ; wavefront will be calculated. ; ;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) ; ApNum: ; The number of pixels to average over in order to get the wavefront ; slopes ; 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 ; PlotStyle: ; 0 = Use an X Window ; 1 = Plot each set of data separately ; 2 = Plot the figures for the paper pro phase_screen_multi_layer_image,$ xWindSpeed=xWindSpeed,yWindSpeed=yWindSpeed,$ PhaseDim=PhaseDim,$ Lo=Lo,lcut=lcut,$ MultiWindSpeeds=MultiWindSpeeds, AvgSub=AvgSub,$ UseFid=UseFid,Covariance=Covariance, TvFlag=TvFlag,$ OldWay=OldWay,MySeed=MySeed,ApDim=ApDim,$ PlotStyle=PlotStyle,ComputeAnyway=ComputeAnyway,$ DifferenceMeth=DifferenceMeth,$ ApNum=ApNum If (not keyword_set(AvgSub)) then AvgSub =1 If (not keyword_set(FittedOnly)) then FittedOnly = 2 If (not keyword_set(UseFid)) then UseFid =2 If (not keyword_set(ApNum)) then ApNum = 6 If (not keyword_set(Covariance)) then Covariance = 1 If (not keyword_set(PhaseDim)) then PhaseDim=long(250) If (not keyword_set(ApDim)) then ApDim=25 If (not keyword_set(Lo)) then Lo=100. If (not keyword_set(DeltaX)) then DeltaX = 0.17/ApNum If (not keyword_set(Ro)) then Ro=0.1 If (not keyword_set(MultiWindSpeeds)) then MultiWindSpeeds=1 If (not keyword_set(ExpTime)) then ExpTime=.03 If (not keyword_set(TvFlag)) then TvFlag = 0 If (not keyword_set(Spacing)) then Spacing = 1 If (not keyword_set(MySeed)) then MySeed= 30 If (not keyword_set(PlotStyle)) then PlotStyle=0 If (not keyword_set(ComputeAnyway)) then ComputeAnyway=0 If (not keyword_set(DifferenceMeth)) then DifferenceMeth=1 ;*************************************************************************** ;FIDUCIAL VALUES NonFiducial=return_fiducial() complement,findgen(626),NonFiducial,Fiducial Border=1 ;Only Need 1 unit for Fi Diff ;************************************************************************* ;FILES ;********************INPUT AtmosphereDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' PostscriptDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/' ScreenDirectory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/lsst_sim/atmosphere/' FinalStatsDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/FinalNumbers/' FileExists=0 ;********************OUTPUT If DifferenceMeth eq 1 then begin DifferenceMethSt="FCD" endif else if DifferenceMeth eq 2 then begin DifferenceMethSt="FFD" endif else if DifferenceMeth eq 3 then begin DifferenceMethSt="FBD" endif else if DifferenceMeth eq 4 then begin DifferenceMethSt="ComPade" endif If UseFid eq 1 then begin FidStr="_Fid_" endif else if UseFid eq 0 then begin FidStr="_NoFid_" endif else begin FidStr='_' endelse ParamStr='MultiLayerPhase'+FidStr+DifferenceMethSt+"_"+'Avg'+strtrim(ApNum,2)+$ 'Aps_'+strtrim(MySeed,2)+'Seed_'+strtrim(Lo,2)+'Low_'+$ strtrim(PhaseDim,2)+'Finite_Exposure_Image' FileName=AtmosphereDir+ParamStr+'DimCovar.txt' If File_Test(FileName) then begin FileExists=1 print,"File: "+FileName+"Already exists...." print,"Taking Data from File" endif If PlotStyle eq 3 then begin PSFileName=PostscriptDir+ParamStr+'.ps' endif else if PlotStyle eq 2 then begin PSFileName=PostscriptDir+ParamStr+'IndPlots.ps' endif Date=['10','11','12'] SOARPowerFileNames=FinalStatsDir+'PowerEllip'+Date+'.txt' SOAREHistFileNames=FinalStatsDir+'EHistogram'+Date+'.txt' ;========================================================================== ;************************************************************************* ;***********************************FILEEXISTS=0 CALCULATE VALUES If ComputeAnyway eq 1 then FileExists=0 If FileExists eq 0 then begin if MultiWindSpeeds eq 1 then begin readcol,'/nfs/slac/g/ki/ki08/lsst/CPanalysis/winds/PhaseScreenWinds.txt',$ NumFiles,Wx1,Wy1,Wgt1,Wx2,Wy2,Wgt2,Wx3,Wy3,Wgt3,Wx4,Wy4,Wgt4,$ format='d,f,f,f,f,f,f,f,f,f,f,f,f' endif ;*******************************************************VALUES FOR LOOP KMax=12 ApDim=25 NumLayers=4 TotalPhaseScreens=Total(NumFiles) PhaseScreens=dblarr(ApDim,ApDim,TotalPhaseScreens) IValArr = dblarr(3, TotalPhaseScreens) ;================================================================================= ;********************************************************************************* ;Loop For Each Segment of the Night NumFile(FNo) Screens -- NumFile(FNo)*NumLayers Seed=MySeed GlobalScreenNum=0 for FNo = 0, N_elements(NumFiles)-1 do begin if NumLayers eq 4 then begin XWinds=[Wx1(FNo),Wx2(FNo),Wx3(FNo),Wx4(FNo)] YWinds=[Wy1(FNo),Wy2(FNo),Wy3(FNo),Wy4(FNo)] CN2s=[Wgt1(FNo),Wgt2(FNo),Wgt3(FNo),Wgt4(FNo)] endif else if NumLayers eq 1 then begin XWinds=[Wx4(FNo)] YWinds=[Wy4(FNo)] CN2s=[Wgt1(FNo)] endif for LocalScreenNum=0, NumFiles(FNo)-1 do begin ;Final Phase at pupil Screens=fltarr(PhaseDim,PhaseDim,NumLayers) ;Generate NumLayers Phase Screens for LayerScreenNum=0, NumLayers-1 do begin ;Ind. Phase Layers Screens(*,*,LayerScreenNum)=return_kolmogorovscreen(PhaseDim,Deltax,$ Lo,Ro,Seed) endfor ;Loop through the time increment and layers and get one value for Ixx, Iyy, Ixy IValArr(*,GlobalScreenNum) = return_mullayphasescreen_image_ellip(Screens, $ XWinds, YWinds, CN2s, $ ExpTime, NumLayers, PhaseDim, ApNum, ApDim, DifferenceMeth, TvFlag) GlobalScreenNum=GlobalScreenNum+1 print, GlobalScreenNum endfor print,FNo endfor ;Form Ellipticity Array to Pass on to the Hist_Ellipticities function EllipticityArr = Fltarr(4, TotalPhaseScreens) EllipticityArr(0,*) = (IValArr(0,*) - IValArr(1,*))/(IValArr(0,*) + IValArr(1,*)) EllipticityArr(1,*) = (IValArr(2,*))/(IValArr(0,*) + IValArr(1,*)) EllipticityArr(2,*) = sqrt(EllipticityArr(0,*)^2+EllipticityArr(1,*)^2) EllipticityArr(3,*) = .5*atan(EllipticityArr(1,*), EllipticityArr(0,*)) ;================================================================================= ;********************WRITE VALUES TO FILE HistEllip = Hist_Ellipticities_Image(EllipticityArr) EHistFileName=AtmosphereDir+ParamStr+'EHistogramImage.txt' get_lun,E openw,E,EHistFileName printf,E,format='(%"%s\t%s\t%s\t%s")',$ 'E Locs','E Hist','Theta Locs', 'ThetaHist' for i=0, N_elements(HistEllip(0,*))-1 do begin printf,E,$ format='(%"%f\t%f\t%f\t%f")',$ HistEllip(0,i),HistEllip(1,i),HistEllip(2,i),HistEllip(3,i) endfor close,E free_lun,E SOARDataFileName='/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt' endif else begin ;=========================================================================== ;*************************************************************************** ;**********************************************READ FROM FILE FILEEXISTS=1 EHistFileName=AtmosphereDir+ParamStr+'EHistogramImage.txt' SOARDataFileName='/nfs/slac/g/ki/ki08/lsst/CPanalysis/Ellip_Wind_Pointing.txt' endelse ;========================================================================== ;**********************************************PLOTTING If PlotStyle eq 1 then begin Set_plot,'z' device,set_resolution=[1400,1200] !p.charsize=3 !p.charthick=5 !x.thick=4 !y.thick=4 !p.thick=3.5 !p.noerase=0 loadct, 39 white='FFFFFF'x black='000000'x red='FF0000'x !P.multi=[0,1,1] endif else if PlotStyle eq 2 or PlotStyle eq 3 then begin Set_plot,'ps' device,filename=PSFileName loadct,39 ;load color table 39 device,/color ;allow color on the postscript device,ysize=8.5,/inches ;Height of plot in y device,xsize=12.0,/inches ;Width of plot in x device,yoffset=0.0,/inches ;Y position of lower left corner device,xoffset=0.0,/inches !p.charsize=1 !p.thick=4.5 !p.charthick=4 !x.thick=4 !y.thick=4 white='FFFFFF'x black='000000'x red='FF0000'x green='00FF00'x blue='0000FF'x endif else begin set_plot,'X' endelse PlotE=plot_elliphistogram_ff(EHistFileName,SOARDataFileName,PlotStyle,ParamStr,$ PostScriptDir) if PlotStyle eq 2 or PlotStyle eq 3 then begin device,/close endif set_plot,'x' stop end