;spot_drift ; ;PURPOSE: ; This script will read in all of the shackhartman data from a series ; of images with path given by ; ; [Directory]hr5132[letter]_file_number_stats.txt ; ; and record into a data cube the centroids of the spots ; ;CALLING SEQUENCE ; ; spot_drift,[Directory=Directory,[letter=letter],[id_array]]] ; ;INPUTS ; ; Directory The path to the directory where the Shackhartman data ; are located ; letter The letter included in the files. This ranges from ; a-j, corresonding to early-late times ; id_array The id numbers (stored in a one-dimensional array) ; of the points that should be plotted pro shackhartman, Directory=Directory,letter=letter,id_array=id_array ;********************************************************************* ; ;Here are 27 x reference points and 27 y reference points that define where the spot actually ;should be. There are a total of 729 possible points x_ref_col=[225,250,278,302,332,353,380,404,428,452,476,502,528,$ 554,578,604,630,654,680,704,730,754,780,806,832] y_ref_col=[164,188,214,238,264,290,314,340,366,390,414,440,464,$ 490,516,542,568,592,616,642,668,692,718,742,768] spot_count=long(N_elements(x_ref_col)*N_elements(y_ref_col)) ;********************************************************************* ;set directory path delimiter if (!version.os_family eq 'unix') then begin delim='/' endif else begin delim='\' endelse If (not keyword_set(Directory)) then Directory='/a/pippin01/Volumes/u08/lsst/CPanalysis/2005-05-10/ShackHartman/' If (not keyword_set(id_array)) then id_array=[300,301,302] If (not keyword_set(letter)) then letter='a' stack_x=findgen(spot_count,99) stack_y=findgen(spot_count,99) for file_number=1,99 do begin if file_number lt 10 then begin mag_file_name=Directory+'hr5132'+letter+'_'+$ string(file_number,format='(I1)')+'_raw.fits.mag.1.txt' coo_file_name=Directory+'hr5132'+letter+'_'+$ string(file_number,format='(I1)')+'_raw.fits.coo.1' endif if file_number ge 10 and file_number lt 100 then begin mag_file_name=Directory+'hr5132'+letter+'_'+$ string(file_number,format='(I2)')+'_raw.fits.mag.1.txt' coo_file_name=Directory+'hr5132'+letter+'_'+$ string(file_number,format='(I2)')+'_raw.fits.coo.1' endif if file_number ge 100 and file_number lt 1000 then begin mag_file_name=Directory+'hr5132'+letter+'_'+$ string(file_number,format='(I3)')+'_raw.fits.mag.1.txt' coo_file_name=Directory+'hr5132'+letter+'_'+$ string(file_number,format='(I3)')+'_raw.fits.coo.1' endif ;Read in the data columns from the output IRAF mag file. They are ; XCENTER, YCENTER, XERR, YERR, MAG readcol,mag_file_name,x_coord,y_coord,x_err,y_err,mag,format='f,f,f,f,f' ;Read in the data columns from the output IRAF coo file. They are ; SHARP, ROUNDNESS readcol,coo_file_name,sharp,roundness,format='X,X,X,f,f',SKIPLINE=41 ;This loop will go through all of the possible places where a point could be ;and check to see if there is an image point near it for col_num=0,N_elements(y_ref_col)-1 do begin for row_num=0,N_elements(x_ref_col)-1 do begin ;Pick points where x might be near and those where y might ;be near x_index_pos=where(abs(x_coord-x_ref_col(row_num)) lt 5) y_index_pos=where(abs(y_coord-y_ref_col(col_num)) lt 5) for y_sub_index=0,N_elements(y_index_pos)-1 do begin x_sub_index=where(x_index_pos eq $ y_index_pos(y_sub_index)) if ((x_sub_index(0) ne -1) and (y_index_pos(0) ne -1) and (x_index_pos(0) ne -1)) then begin x_index=x_index_pos(x_sub_index) y_index=y_index_pos(y_sub_index) x_ref_dif=x_coord(x_index)-x_ref_col(row_num) y_ref_dif=y_coord(y_index)-y_ref_col(col_num) printf,1,FORMAT='(%"%3.0d\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f\t%0.3f")',col_num*N_elements(x_ref_col)+$ row_num,x_coord(x_index),y_coord(y_index),$ x_err(x_index),y_err(y_index),mag(x_index),$ sharp(x_index),roundness(x_index),$ x_ref_dif,y_ref_dif found_point=1 break endif else begin found_point=0 endelse endfor if found_point eq 0 then printf,1,$ FORMAT='(%"%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d")',$ col_num*N_elements(x_ref_col)+row_num,-1,-1,-1,-1,-1,-1,-1,-1,-1 endfor endfor ;create data cube of (x,y) coordinate pairs stack_x(*,file_number-1)=x_coord stack_y(*,file_number-1)=y_coord endfor weighted_x_arr=fltarr(spot_count) weighted_y_arr=fltarr(spot_count) variance_x_arr=fltarr(spot_count) variance_y_arr=fltarr(spot_count) valid_index=0 set_plot,'ps' device,filename='hr5132'+letter+'_reference'+'.ps' 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=7.0,/inches ;Width of plot in x device,yoffset=1.0,/inches ;Y position of lower left corner !P.CHARSIZE=2 !P.THICK=10. plot,[200,900],[100,900],/NODATA,title='Plot of ShackHartman Focal Plane',xtitle='X coordinate',ytitle='Y coordinate',background=white,color=000,/xstyle,/ystyle,position=[.1,.1,.9,.9] ;For each spot, go through and calculate its average (x,y) coordinate for spot_number=0,spot_count-1 do begin ;Any spots where the x_coord and y_coord are -1 will not be ;included. valid_point_indices are the array indices where this ;is the case valid_point_indices=where((stack_x(spot_number,*) gt 0) and (stack_y(spot_number,*) gt 0)) ;print,N_elements(valid_point_indices) ;For some spots, there will be no valid fits if (N_elements(valid_point_indices) gt 70) then begin for valid_point_index=0, N_elements(valid_point_indices)-1 do begin oplot,[stack_x(spot_number,valid_point_index)],$ [stack_y(spot_number,valid_point_index)],$ psym=3,color=(float(max(valid_point_indices)-valid_point_index)/max(valid_point_indices)*256) endfor weighted_x=mean(stack_x(spot_number,valid_point_indices)) weighted_y=mean(stack_y(spot_number,valid_point_indices)) variance_x=variance(stack_x(spot_number,valid_point_indices)) variance_y=variance(stack_y(spot_number,valid_point_indices)) weighted_x_arr(valid_index)=weighted_x weighted_y_arr(valid_index)=weighted_y variance_x_arr(valid_index)=variance_x variance_y_arr(valid_index)=variance_y valid_index=valid_index+1 print,N_elements(valid_point_indices) endif endfor weighted_x_arr=weighted_x_arr(where(weighted_x_arr gt 0)) weighted_y_arr=weighted_y_arr(where(weighted_y_arr gt 0)) variance_x_arr=variance_x_arr(where(weighted_x_arr gt 0)) variance_y_arr=variance_y_arr(where(weighted_y_arr gt 0)) oplot,weighted_x_arr,weighted_y_arr,psym=3,color=256 device,/close stop end