pro phase_screen_ellip_theta_histogram,$ Date=Date,stats_directory=stats_directory,$ depth=depth,x_off=x_off,y_off=y_off,exam_date=exam_date If (not keyword_set(Date)) then Date='10' If (not keyword_set(exam_date)) then exam_date='1_31' If (not keyword_set(stats_directory)) then stats_directory='stats_'+exam_date+'/' If (not keyword_set(depth)) then depth=0 If (not keyword_set(postscript_dir)) then postscript_dir='distorted_postscripts' If (not keyword_set(x_off)) then x_off=0 If (not keyword_set(y_off)) then y_off=0 postscript_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/' nonfiducial=[1,2,3,4,5,6,7,8,9,13,14,17,18,19,20,21,22,23,24,25,$ 26,27,28,29,30,31,32,38,39,45,46,47,48,49,50,$ 51,52,53,54,55,63,64,71,72,73,74,75,$ 76,77,78,79,88,89,98,98,99,100,$ 101,102,103,113,114,123,124,125,$ 126,127,128,138,139,149,150,$ 151,152,163,164,174,175,$ 176,177,188,189,200,$ 201,210,213,214,225,$ 226,238,239,249,$ 251,261,262,263,264,265,266,$ 276,286,287,288,289,290,291,$ 301,302,303,304,305,306,307,308,309,$ 310,311,312,313,314,315,316,$ 317,318,319,320,321,322,323,324,325,$ 326,327,336,337,338,339,340,341,$ 351,361,362,363,364,365,366,$ 376,388,389,$ 401,413,414,425,$ 426,427,438,439,449,450,$ 451,452,454,463,464,473,474,475,$ 476,477,478,488,489,498,499,500,$ 501,502,503,504,513,514,522,523,524,525,$ 526,527,528,529,530,538,539,546,547,548,549,550,$ 551,552,553,554,555,556,563,564,570,571,572,573,574,575,$ 576,577,578,579,580,581,582,588,589,594,595,596,597,598,599,600,$ 601,602,603,604,605,606,607,608,609,610,613,614,617,618,619,620,621,$ 622,623,624,625]-1 complement,findgen(625),nonfiducial,fiducial fitted=fiducial ;ARRAYS OF FILES all_nights=1 if all_nights eq 1 then begin stats_directory1='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-10'+$ '/ShackHartman/'+stats_directory+'/' offset_files1=file_search(stats_directory1+'*5132*sub_offsets.txt',$ count=num_stats_files1) stats_directory2='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-11'+$ '/ShackHartman/'+stats_directory+'/' offset_files2=file_search(stats_directory2+'*5132*sub_offsets.txt',$ count=num_stats_files2) stats_directory3='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-12'+$ '/ShackHartman/'+stats_directory+'/' offset_files3=file_search(stats_directory3+'*5132*sub_offsets.txt',$ count=num_stats_files3) num_stats_files=num_stats_files1+num_stats_files2+num_stats_files3 offset_files=[offset_files1,offset_files2,offset_files3] Date='All_Nights' endif else begin stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/'+stats_directory+'/' offset_files=file_search(stats_directory+'*5132*sub_offsets.txt',$ count=num_stats_files) endelse ;CONSTANTS num=25 grid_spots=625 index=findgen(625) ;ARRAY REFERENCE 25x25=625 FOR GRID subx=2*indgen(num)-25 ;The x-reference numbers suby=indgen(num)-12 ;The y-reference numbers ;Grid for partvelvec x_ref_col=12.5*subx+542 y_ref_col=25.0*suby+460 index=findgen(num*num) gridx=fltarr(num*num) gridy=fltarr(num*num) gridx(index)=x_ref_col(index mod num) gridy(index)=y_ref_col(index/num) ;Indices and such no_bins=8 a=findgen(13)+12 b=findgen(12) subn=[a,b] subm=subn ;HUGE DATA CUBE x_offset_hor=fltarr(num_stats_files,num,num) y_offset_hor=fltarr(num_stats_files,num,num) x_offset_ver=fltarr(num_stats_files,num,num) y_offset_ver=fltarr(num_stats_files,num,num) ;Correlations timee1=fltarr(num_stats_files-depth,num,num) timee2=fltarr(num_stats_files-depth,num,num) timee1num=fltarr(num_stats_files-depth,num,num) timee2num=fltarr(num_stats_files-depth,num,num) timeeden=fltarr(num_stats_files-depth,num,num) timee1_arr=fltarr(num,num) timee2_arr=fltarr(num,num) ;for high pass filter frame_e1=fltarr(num_stats_files) frame_e2=fltarr(num_stats_files) frame_e=fltarr(num_stats_files) frame_thetap=fltarr(num_stats_files) frame_thetam=fltarr(num_stats_files) frame_e_den=fltarr(num_stats_files) kx0_array=fltarr(num_stats_files) ky0_array=fltarr(num_stats_files) phi_array=fltarr(num_stats_files) r_array=fltarr(num_stats_files) ;FILL THE DATA CUBES for file_number=0,num_stats_files-1 do begin ;Take the results from the offset.txt files and fill the cube readcol,offset_files(file_number),x_offset,y_offset,format='X,X,f,X,f' ;Obtain the difference ordered in COLUMN, ROW between centroid and grid point x_offset(where(x_offset eq -10))=0 y_offset(where(y_offset eq -10))=0 avg_sub=1 if avg_sub eq 1 then begin x_offset=x_offset-mean(x_offset(where(x_offset ne 0))) y_offset=y_offset-mean(y_offset(where(y_offset ne 0))) endif else if avg_sub eq 2 then begin x_fft=complexarr(num,num) y_fft=complexarr(num,num) x_fft(*,*)=fft(x_offset) y_fft(*,*)=fft(y_offset) x_fft_re=complexarr(num,num) y_fft_re=complexarr(num,num) for n=0,num-1 do begin for m=0,num-1 do begin x_fft_re(subn(n),subm(m))=x_fft(n,m) y_fft_re(subn(n),subm(m))=y_fft(n,m) endfor endfor kx0_array(file_number)=x_fft_re(12,12) ky0_array(file_number)=y_fft_re(12,12) phi_array(file_number)=atan(y_fft_re(12,12),x_fft_re(12,12)) r_array(file_number)=sqrt(x_fft_re(12,12)^2+y_fft_re(12,12)^2) x_fft_re(12,12)=0 y_fft_re(12,12)=0 for n=0,num-1 do begin for m=0,num-1 do begin x_fft(n,m)=x_fft_re(subn(n),subm(m)) y_fft(n,m)=y_fft_re(subn(n),subm(m)) endfor endfor x_offset=fft(x_fft,/inverse) y_offset=fft(y_fft,/inverse) endif x_offset_hor(file_number,*,*)=x_offset y_offset_hor(file_number,*,*)=y_offset ;Obtain the difference ordered in ROW, COLUMN x_offset_ver(file_number,*,*)=x_offset((index mod 25)*25+index/25) y_offset_ver(file_number,*,*)=y_offset((index mod 25)*25+index/25) endfor readcol,offset_files(0),xcoord,ycoord,format='x,f,x,f,x' for file_number=0,num_stats_files-depth-1 do begin timee1num(file_number,*,*)=x_offset_hor(file_number,*,*)*$ x_offset_hor(file_number+depth,*,*)-$ y_offset_hor(file_number,*,*)*$ y_offset_hor(file_number+depth,*,*) timee2num(file_number,*,*)=(x_offset_hor(file_number,*,*)*$ y_offset_hor(file_number+depth,*,*)$ +x_offset_hor(file_number+depth,*,*)$ *y_offset_hor(file_number,*,*)) timeeden(file_number,*,*)=(x_offset_hor(file_number,*,*)^2+$ y_offset_hor(file_number,*,*)^2) timee1temp=timee1num(file_number,*,*) timee2temp=timee2num(file_number,*,*) timeedtemp=timeeden(file_number,*,*) frame_e_den=mean(timeedtemp(where(timee1temp ne 0))) frame_e1(file_number)=mean(timee1temp(where(timee1temp ne 0)))/$ mean(timeedtemp(where(timee1temp ne 0))) frame_e2(file_number)=mean(timee2temp(where(timee2temp ne 0)))/$ mean(timeedtemp(where(timee2temp ne 0))) frame_e(file_number)=sqrt(frame_e1(file_number)^2+frame_e2(file_number)^2) signe1=abs(frame_e1(file_number))/frame_e1(file_number) signe2=abs(frame_e2(file_number))/frame_e2(file_number) ;Get the position angle associated with the elongation of the object ;Should be between pi/4 and pi/2 if e1 is negative and between 0 and pi/2 if ;e1 is positive. Should be in second quadrant if e2 is negative frame_thetap(file_number)=$ signe2*asin(sqrt(.5*(1+signe1*sqrt(1-frame_e2(file_number)^2/$ frame_e(file_number)^2)))) ;Recalculate e1 and e2 with new cosine angle frame_e1(file_number)=frame_e(file_number)*$ cos(2*frame_thetap(file_number)) frame_e2(file_number)=frame_e(file_number)*$ sin(2*frame_thetap(file_number)) endfor thetamin=-!pi/2 no_bins=20. thetabinwidth=!pi/no_bins thetamax=thetamin+(no_bins-1)*thetabinwidth thetahist=histogram(frame_thetap,Nbins=no_bins,locations=loc_theta,$ reverse_indices=theta_ind,$ max=thetamax,min=thetamin) thetabinsize=(thetamax-thetamin)/(no_bins-1) set_plot,'ps' device,filename=Postscript_dir+'Theta_Histogram_5_'+Date+'.ps' loadct,0 ;load color table 39 device,/color ;allow color on the postscript device,ysize=11.0,/inches ;Height of plot in y device,xsize=8.5,/inches ;Width of plot in x device,yoffset=1.0,/inches ;Y position of lower left corner phimin=-!pi no_bins=20. phibinwidth=2*!pi/no_bins phimax=phimin+(no_bins-1)*phibinwidth phihist=histogram(phi_array,Nbins=no_bins,locations=loc_phi,$ reverse_indices=phi_ind,$ max=phimax,min=phimin) phibinsize=(phimax-phimin)/(no_bins-1) png=1 if png eq 1 then begin set_plot,'z' erase device, set_font='Courier' device,set_resolution=[800,600] !p.charsize=.8 !p.charthick=1.2 !x.thick=2 !y.thick=2 !p.thick=.8 !p.noerase=0 endif white='FFFFFF'x black='000000'x !P.CHARSIZE=1.0 !P.THICK=4. !p.multi=[0,1,1] plot,loc_theta+thetabinsize/2,thetahist,psym=10,xtitle='Theta',$ ytitle='Number of Events',$ title='Histogram of Ellipticity Angle for night of 5-'+Date,$ background=white,color=black xyouts,.4,.3,'Mean='+Strtrim(mean(frame_thetap),2),/normal,color=black if png eq 1 then begin jpgimg=tvrd() tvlct,reds,greens,blues,/get if avg_sub eq 0 then begin write_png,postscript_dir+'Theta_Histogram_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else if avg_sub eq 1 then begin write_png,postscript_dir+'Theta_Histogram_Avg_sub_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else begin write_png,postscript_dir+'Theta_Histogram_kx0_ky0_sub_5-'+$ Date+'-2005.png',jpgimg,reds,greens,blues endelse endif set_plot,'z' erase white='FFFFFF'x black='000000'x !P.CHARSIZE=1.0 !P.THICK=4. !p.multi=[0,1,1] plot,loc_phi+phibinsize/2,phihist,psym=10,xtitle='Theta',$ ytitle='Number of Events',$ title='Histogram of Centroid Angle for night of 5-'+Date,$ background=white,color=black xyouts,.4,.3,'Mean='+Strtrim(mean(phi_array),2),/normal,color=black if png eq 1 then begin jpgimg=tvrd() tvlct,reds,greens,blues,/get if avg_sub eq 0 then begin write_png,postscript_dir+'Phi_Histogram_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else if avg_sub eq 1 then begin write_png,postscript_dir+'Phi_Histogram_Avg_sub_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else begin write_png,postscript_dir+'Phi_Histogram_kx0_ky0_sub_5-'+$ Date+'-2005.png',jpgimg,reds,greens,blues endelse endif erase white='FFFFFF'x black='000000'x !P.CHARSIZE=1.0 !P.THICK=4. !p.multi=[0,1,1] plot,abs(phi_array),abs(frame_thetap),psym=2,$ xtitle='Azimuthal angle of centroid',$ ytitle='Ellipticity angle that maximizes e1 (Theta)',$ title='Scatter Plot of Theta vs. Phi of kx=0,ky=0 for night of 5-'+Date,$ background=white,color=black if png eq 1 then begin jpgimg=tvrd() tvlct,reds,greens,blues,/get if avg_sub eq 0 then begin write_png,postscript_dir+'Theta_Phi_Abs_Scatter_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else if avg_sub eq 1 then begin write_png,postscript_dir+'Theta_Phi_Abs_Scatter_Avg_sub_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else begin write_png,postscript_dir+'Theta_Phi_Abs_Scatter_kx0_ky0_sub_5-'+$ Date+'-2005.png',jpgimg,reds,greens,blues endelse endif erase white='FFFFFF'x black='000000'x !P.CHARSIZE=1.0 !P.THICK=4. !p.multi=[0,1,1] plot,phi_array,frame_thetap,psym=2,$ xtitle='Azimuthal angle of centroid',$ ytitle='Ellipticity angle that maximizes e1 (Theta)',$ title='Scatter Plot of Theta vs. Phi of kx=0,ky=0 for night of 5-'+Date,$ background=white,color=black if png eq 1 then begin jpgimg=tvrd() tvlct,reds,greens,blues,/get if avg_sub eq 0 then begin write_png,postscript_dir+'Theta_Phi_Scatter_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else if avg_sub eq 1 then begin write_png,postscript_dir+'Theta_Phi_Scatter_Avg_sub_5-'+Date+'-2005.png',$ jpgimg,reds,greens,blues endif else begin write_png,postscript_dir+'Theta_Phi_Scatter_kx0_ky0_sub_5-'+$ Date+'-2005.png',jpgimg,reds,greens,blues endelse endif device,/close set_plot,'x' stop end