;xy_offsets ; ;PURPOSE: ; This script will return the difference in coordinates between a ; ShackHartman image centroid and the corresponding point on the ; reference grid. I.e. it will give the offsets between the ; Image centroid and the supposed true center point of a ShackHartman ; spot ; ; It will also calculate the magnitude and angle of the deflection ; vector. ; ;OUTPUT: ; The program will write to the file ; [stats_directory]/[raw_file_name].offsets.txt ; ; The following columns will be written ; ; index|xcoord|ycoord|xoffset|yoffset|magnitude|angle pro routine_xy_offsets,Date=Date,stats_directory=stats_directory,$ x_center_file=x_center_file,y_center_file=y_center_file,$ exam_date=exam_date,HalfSpot=HalfSpot If (not keyword_set(exam_date)) then exam_Date='1_31' If (not keyword_set(Date)) then Date='10' If (not keyword_set(stats_directory)) then stats_directory='stats_'+exam_date+'/' If (not keyword_set(x_center_file)) then x_center_file='newest_x_centers.txt' If (not keyword_set(y_center_file)) then y_center_file='newest_y_centers.txt' If (not keyword_set(HalfSpot)) then HalfSpot = 5 ;-----------------------------------------STATS AND RAW .fits IMAGES stats_directory='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/'+stats_directory+'/' RawDir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+Date+$ '/ShackHartman/raw/' Stats_files=file_search(stats_Directory+'*5132*_stats.txt',count=num_of_stats_files) RawFiles=file_search(RawDir+'h*.fits',count=NumRawFiles) ;OBTAIN THE FITTED CENTERS FOR ALL IMAGES readcol,stats_Directory+x_center_file,x_fit_center,format='X,X,f' readcol,stats_Directory+y_center_file,y_fit_center,format='X,X,f' ;-----------------------------------------------------------CONSTANTS num_of_files=N_elements(x_fit_center) num=25 ; subx=2*indgen(num)-25 ;The x-reference numbers suby=indgen(num)-12 ;The y-reference numbers ;GENERATE LARGE ARRAYS TO HOLD THE AVERAGE OFFSETS OF THE SPOTS xoffset_array=fltarr(num*num,num_of_files) yoffset_array=fltarr(num*num,num_of_files) xcoord_array=fltarr(num*num,num_of_files) ycoord_array=fltarr(num*num,num_of_files) mag_array=fltarr(num*num,num_of_files) bolo_mag_array=fltarr(num*num,num_of_files) phi_array=fltarr(num*num,num_of_files) Xerr_array=fltarr(num*num,num_of_files) Yerr_array=fltarr(num*num,num_of_files) Sharp_array=fltarr(num*num,num_of_files) print,num_of_stats_files print,num_of_files ;-------------------------------------------------------------------- ;=================================================================== ;FIRST LOOP ----EXTRACT ORIGINAL OFFSETS for file_number=0,num_of_files-1 do begin stats_path_parts=strsplit(stats_files(file_number),'/',/extract) stats_file_ref=strsplit(stats_path_parts(N_elements(stats_path_parts)-1)$ ,'stats.txt',/extract,/regex) stats_file_ref=stats_file_ref(0) ;Take the centroid values from the stats.txt files stats_file_name=stats_files(file_number) offsets_file_name=stats_directory+stats_file_ref+'offsets.txt' readcol,stats_file_name,grid_index,x_coord,y_coord,$ XErr, YErr, bolo_mags, Sharp, $ format='f,f,f,f,f,f,f' grid_spots=N_elements(x_coord) ;Figure out the new reference grid around x_center and y_center x_ref_col=12.5*subx+x_fit_center(file_number) y_ref_col=25.0*suby+y_fit_center(file_number) 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) not_fitted=where(x_coord eq -1,complement=fitted) ;Obtain the difference between centroid and grid point x_offset=x_coord-gridx y_offset=y_coord-gridy magnitude=sqrt(x_offset^2+y_offset^2) phi=atan(y_offset,x_offset) ;If point was not fitted, set to -1 x_offset(not_fitted)=-10 y_offset(not_fitted)=-10 y_coord(not_fitted)=-10 x_coord(not_fitted)=-10 magnitude(not_fitted)=-10 phi(not_fitted)=-10 bolo_mags(not_fitted)=-10 Xerr(not_fitted)=-10 Yerr(not_fitted)=-10 Sharp(not_fitted)=-10 xoffset_array(*,file_number)=x_offset yoffset_array(*,file_number)=y_offset xcoord_array(*,file_number)=x_coord ycoord_array(*,file_number)=y_coord mag_array(*,file_number)=magnitude bolo_mag_array(*,file_number)=bolo_mags phi_array(*,file_number)=phi Xerr_array(*,file_number)=Xerr Yerr_array(*,file_number)=Yerr Sharp_array(*,file_number)=Sharp get_lun,offsets openw,offsets,offsets_file_name for stat_index=0, num*num-1 do begin printf,offsets,FORMAT='(%"%i\t%f\t%f\t%f\t%f\t%f\t%f\t%f")',$ stat_index,$ x_coord(stat_index),x_offset(stat_index),$ y_coord(stat_index),y_offset(stat_index),$ magnitude(stat_index),phi(stat_index),bolo_mags(stat_index) endfor close,offsets free_lun,offsets endfor ;------------------------------------------------------------------------- ;SECOND LOOP----COMPUTE AVERAGE OFFSETS AND COORDS avg_xoffset_array=fltarr(num*num) avg_yoffset_array=fltarr(num*num) get_lun,offsets openw,offsets,stats_directory+'Average_offsets.txt' for stat_index=0,num*num-1 do begin ;NOW OBTAIN THE AVERAGE OFFSET FOR EACH SPOT xoffset=xoffset_array(stat_index,*) yoffset=yoffset_array(stat_index,*) xcoord=xcoord_array(stat_index,*) ycoord=ycoord_array(stat_index,*) valid=where(xoffset ne -10) if valid(0) eq -1 then begin avg_xoffset=-1 avg_yoffset=-1 avg_xcoord=-1 avg_ycoord=-1 endif else begin avg_xoffset=mean(xoffset(valid)) avg_yoffset=mean(yoffset(valid)) avg_xcoord=mean(xcoord(valid)) avg_ycoord=mean(ycoord(valid)) endelse avg_xoffset_array(stat_index)=avg_xoffset avg_yoffset_array(stat_index)=avg_yoffset printf,offsets,FORMAT='(%"%i\t%f\t%f\t%f\t%f")',$ stat_index,$ avg_xcoord,avg_xoffset,$ avg_ycoord,avg_yoffset endfor close,offsets free_lun,offsets ;---------------------------------------------------------------------- ;THIRD LOOP avg_xoffset100_array=fltarr(num*num) avg_yoffset100_array=fltarr(num*num) rms_xoffset100_array=fltarr(num*num) rms_yoffset100_array=fltarr(num*num) AvgXCooArr=fltarr(num*num) AvgYCooArr=fltarr(num*num) AvgXerrArr=fltarr(num*num) AvgYerrArr=fltarr(num*num) AvgSharpArr=fltarr(num*num) num_stats_files=num_of_files-1 ;***************************************************SHOW AVG DIS. BY 100 tv_off=0 ;SET TV_OFF=0 to skip if tv_off eq 1 then begin set_plot,'z' erase device,set_resolution=[1200,900] !p.charsize=2 !p.charthick=1.8 !x.thick=3 !y.thick=3 !p.thick=1.2 !p.color=0 loadct,39 white='FFFFFF'x black='000000'x red='FF0000'x !p.noerase=0 image_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/' plot,[100,820],[180,880],/xstyle,/ystyle,xrange=[100,820],yrange=[180,880],$ /nodata,color=black,background=white,ytitle='Y coordinate (pixels)',$ xtitle='X coordinate (pixels)' xyouts,200,750,'1 pix' !p.noerase=1 for time_ind=0, num_stats_files/100-1 do begin for stat_index=0,num*num-1 do begin ;NOW OBTAIN THE AVERAGE OFFSET FOR EACH SPOT startInd=time_ind*100 finInd=(time_ind+1)*100 xoffset=xoffset_array(stat_index,startInd:finInd) yoffset=yoffset_array(stat_index,startInd:finInd) xcoord=xcoord_array(stat_index,startInd:finInd) ycoord=ycoord_array(stat_index,startInd:finInd) valid=where(xoffset ne -10) if valid(0) eq -1 then begin avg_xoffset=-1 avg_yoffset=-1 avg_xcoord=-1 avg_ycoord=-1 endif else begin avg_xoffset=mean(xoffset(valid)) avg_yoffset=mean(yoffset(valid)) avg_xcoord=mean(xcoord(valid)) avg_ycoord=mean(ycoord(valid)) endelse AvgXCooArr(stat_index)=avg_xcoord AvgYCooArr(stat_index)=avg_ycoord avg_xoffset100_array(stat_index)=avg_xoffset avg_yoffset100_array(stat_index)=avg_yoffset endfor AvgXCooArr=[AvgXCooArr,220] AvgYCooArr=[AvgYCooArr,750] avg_xoffset100_array=[avg_xoffset100_array,1] avg_yoffset100_array=[avg_yoffset100_array,0] valid=where(AvgXCooArr ne -1) partvelvec,$ avg_xoffset100_array(valid),$ avg_yoffset100_array(valid),$ AvgXCooArr(valid),AvgYCooArr(valid),$ yrange=[100,820],xrange=[180,880],$ background=white,$ color=240-floor(240.*(time_ind*100./num_stats_files)),$ xstyle=5,ystyle=5 endfor jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,image_dir+'Average_Distortions_By100_'+$ strtrim(string(Date),2)+'.png',$ jpgimg,reds,greens,blues erase ;GET RMS OF X AND Y OFFSETS !p.noerase=0 image_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/final_paper_figs/' plot,[100,820],[180,880],/xstyle,/ystyle,xrange=[100,820],yrange=[180,880],$ /nodata,color=black,background=white,ytitle='Y coordinate (pixels)',$ xtitle='X coordinate (pixels)' xyouts,200,750,'1 pix' !p.noerase=1 for time_ind=0, num_stats_files/100-1 do begin for stat_index=0,num*num-1 do begin startInd=time_ind*100 finInd=(time_ind+1)*100 xoffset=xoffset_array(stat_index,startInd:finInd) yoffset=yoffset_array(stat_index,startInd:finInd) xcoord=xcoord_array(stat_index,startInd:finInd) ycoord=ycoord_array(stat_index,startInd:finInd) valid=where(xoffset ne -10) if valid(0) eq -1 then begin rms_xoffset=-1 rms_yoffset=-1 avg_xcoord=-1 avg_ycoord=-1 endif else begin rms_xoffset=sqrt(mean((xoffset(valid))-avg_xoffset_array(stat_index))^2) rms_yoffset=sqrt(mean((yoffset(valid))-avg_yoffset_array(stat_index))^2) avg_xcoord=mean(xcoord(valid)) avg_ycoord=mean(ycoord(valid)) endelse AvgXCooArr(stat_index)=avg_xcoord AvgYCooArr(stat_index)=avg_ycoord rms_xoffset100_array(stat_index)=rms_xoffset rms_yoffset100_array(stat_index)=rms_yoffset endfor AvgXCooArr=[AvgXCooArr,220] AvgYCooArr=[AvgYCooArr,750] rms_xoffset100_array=[rms_xoffset100_array,1] rms_yoffset100_array=[rms_yoffset100_array,0] valid=where(AvgXCooArr ne -1) partvelvec,$ rms_xoffset100_array(valid),$ rms_yoffset100_array(valid),$ AvgXCooArr(valid),AvgYCooArr(valid),$ yrange=[100,820],xrange=[180,880],$ background=white,$ color=240-floor(240.*(time_ind*100./num_stats_files)),$ xstyle=5,ystyle=5 endfor jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,image_dir+'Average_RMS_By100_'+$ strtrim(string(Date),2)+'.png',$ jpgimg,reds,greens,blues ;GET AVG OF X AND Y ERRORS DEALT BY IRAF DAOFIND erase !p.noerase=0 plot,[100,820],[180,880],/xstyle,/ystyle,xrange=[100,820],yrange=[180,880],$ /nodata,color=black,background=white,ytitle='Y coordinate (pixels)',$ xtitle='X coordinate (pixels)' xyouts,200,750,'.05 pix' !p.noerase=1 for time_ind=0, num_stats_files/100-1 do begin for stat_index=0,num*num-1 do begin ;NOW OBTAIN THE AVERAGE OFFSET FOR EACH SPOT startInd=time_ind*100 finInd=(time_ind+1)*100 Xerr=Xerr_array(stat_index,startInd:finInd) Yerr=Yerr_array(stat_index,startInd:finInd) Sharp=Sharp_array(stat_index,startInd:finInd) xcoord=xcoord_array(stat_index,startInd:finInd) ycoord=ycoord_array(stat_index,startInd:finInd) valid=where(xcoord ne -10) if valid(0) eq -1 then begin Avg_Xerr=-1 Avg_Yerr=-1 Avg_Sharp=-1 avg_xcoord=-1 avg_ycoord=-1 endif else begin Avg_Xerr=mean(Xerr(valid)) Avg_Yerr=mean(Yerr(valid)) Avg_sharp=mean(Sharp(valid)) avg_xcoord=mean(xcoord(valid)) avg_ycoord=mean(ycoord(valid)) endelse AvgXCooArr(stat_index)=avg_xcoord AvgYCooArr(stat_index)=avg_ycoord AvgXErrArr(stat_index)=Avg_Xerr AvgYErrArr(stat_index)=Avg_Yerr AvgSharpArr(stat_index)=Avg_Sharp endfor AvgXCooArr=[AvgXCooArr,220] AvgYCooArr=[AvgYCooArr,750] AvgXErrArr=[AvgXErrArr,.05] AvgYErrArr=[AvgYErrArr,0] valid=where(AvgXCooArr ne -1) partvelvec,$ AvgXErrArr(valid),$ AvgYErrArr(valid),$ AvgXCooArr(valid),AvgYCooArr(valid),$ yrange=[100,820],xrange=[180,880],$ background=white,$ color=240-floor(240.*(time_ind*100./num_stats_files)),$ xstyle=5,ystyle=5,title='X and Y error' endfor jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,image_dir+'Average_XErr_YErr_By100_'+$ strtrim(string(Date),2)+'.png',$ jpgimg,reds,greens,blues erase !p.noerase=0 plot,[100,820],[180,880],/xstyle,/ystyle,xrange=[100,820],yrange=[180,880],$ /nodata,color=black,background=white,ytitle='Y coordinate (pixels)',$ xtitle='X coordinate (pixels)',title='Sharpness for each spot' xyouts,200,750,'1 pix' !p.noerase=1 for stat_index=0,num*num-1 do begin ;NOW OBTAIN THE AVERAGE OFFSET FOR EACH SPOT Sharp=Sharp_array(stat_index,0:num_stats_files) xcoord=xcoord_array(stat_index,0:num_stats_files) ycoord=ycoord_array(stat_index,0:num_stats_files) valid=where(xcoord ne -10) if valid(0) eq -1 then begin Avg_Sharp=-1 avg_xcoord=-1 avg_ycoord=-1 endif else begin Avg_sharp=mean(Sharp(valid)) avg_xcoord=mean(xcoord(valid)) avg_ycoord=mean(ycoord(valid)) AvgYCooArr(stat_index)=avg_ycoord AvgSharpArr(stat_index)=Avg_Sharp plot,[AvgXCooArr(stat_index)],[AvgYCooArr(stat_index)],$ psym=6,symsize=4*Avg_Sharp,xrange=[100,820],yrange=[180,880] endelse endfor jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,image_dir+'Average_Sharpness_'+$ strtrim(string(Date),2)+'.png',$ jpgimg,reds,greens,blues stop endif ;*********************************************************************** ;FOURTH LOOP ;Run through again 1) subtract offsets ; 2) Get Spot Characteristics for file_number=0,num_of_files-1 do begin stats_path_parts=strsplit(stats_files(file_number),'/',/extract) stats_file_ref=strsplit(stats_path_parts(N_elements(stats_path_parts)-1)$ ,'stats.txt',/extract,/regex) stats_file_ref=stats_file_ref(0) ;Take the centroid values from the stats.txt files stats_file_name=stats_files(file_number) offsets_sub_file_name=stats_directory+stats_file_ref+'sub_offsets.txt' fits_read,RawFiles(file_number),Im get_lun,offsets openw,offsets,offsets_sub_file_name Xcoord=xcoord_array(*,file_number) Ycoord=ycoord_array(*,file_number) Spots=where(Xcoord ne -10,complement=NonSpots) Xoffset=xoffset_array(*,file_number) Yoffset=yoffset_array(*,file_number) Xoffset(Spots)=Xoffset(Spots)-avg_xoffset_array(Spots) Yoffset(Spots)=Yoffset(Spots)-avg_yoffset_array(Spots) magnitude=sqrt(x_offset^2+y_offset^2) phi=atan(y_offset,x_offset) bolo_mags=bolo_mag_array(*,file_number) magnitude(NonSpots)=-10 phi(NonSpots)=-10 for SpotNo=0, num*num-1 do begin if XCoord(SpotNo) ne -10 then begin ;USE FITS FILE TO GET SPOTS Spot=Im(floor(Xcoord(SpotNo))-HalfSpot:$ floor(Xcoord(SpotNo))+HalfSpot,$ floor(Ycoord(SpotNo))-HalfSpot:$ floor(Ycoord(SpotNo))+HalfSpot) Spot=Spot/max(Spot) Flux=total(Spot) XSpot=rebin(findgen(2*HalfSpot+1)-HalfSpot,$ 2*HalfSpot+1,2*HalfSpot+1) YSpot=rebin(transpose(findgen(2*HalfSpot+1)-HalfSpot),$ 2*HalfSpot+1,2*HalfSpot+1) ;QUADROPOLE MOMENTS--------------- I11=total(XSpot*XSpot*Spot) I22=total(YSpot*YSpot*Spot) I12=total(XSpot*YSpot*Spot) ;ELLIPTICITY COMPONENTS----------- d=sqrt(.5*(I11+I22)/Flux) e1=(I11-I22)/(I11+I22) e2=2*I12/(I11+I22) e=sqrt(e1^2+e2^2) pa=0.5*atan((I12),(I11-I22)) I11=I11/flux I12=I12/flux I22=I22/flux endif else begin I11=-10 I22=-10 I12=-10 d=-10 e1=-10 e2=-10 e=-10 pa=-10 flux=-10 endelse printf,offsets,$ FORMAT=$ '(%"%i\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f")',$ SpotNo,$ Xcoord(SpotNo),Xoffset(SpotNo),$ Ycoord(SpotNo),Yoffset(SpotNo),$ magnitude(SpotNo),phi(SpotNo),bolo_mags(SpotNo),$ I11,I22,I12,d,e1,e2,e,pa,flux endfor close,offsets free_lun,offsets endfor stop end