;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 tv_spot_ellipticity,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*_sub_offsets.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 print,num_of_stats_files print,num_of_files !p.charsize=3 !p.charthick=2 ;*********************************************************************** ;FOURTH LOOP ;Run through again 1) subtract offsets ; 2) Get Spot Characteristics for file_number=2,num_of_files-1 do begin readcol,Stats_files(file_number),grid_index,XCoord,YCoord,$ XErr, YErr, bolo_mags, Sharp, $ format='i,f,X,f,X,f,f,f,f' 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 Spots=where(XCoord ne -10,complement=NonSpots) 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) window,xsize=512,ysize=600 !p.noerase=1 whiteback=fltarr(88,512) whiteback(*)=1000 tv,congrid(whiteback,512,88),0,513 tvscl,congrid(Spot,512,512,/center) 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)) plot,[e,e],[pa-!pi,pa],/polar,color='00FF00'x,$ xrange=[-.2,.2],yrange=[-.2,.2],$ xstyle=5,ystyle=5,$ position=[0,0,1,.8] xyouts,.03,.95,'e='+strtrim(string(FORMAT='(%"%8.5f")',e),2),/norm,$ color='000000'x xyouts,.03,.875,'Angle='+$ strtrim(string(FORMAT='(%"%10.2f")',180*pa/!pi),2),/norm,$ color='000000'x xyouts,.5,.95,'e1='+$ strtrim(string(FORMAT='(%"%8.5f")',e1),2),/norm,$ color='000000'x xyouts,.5,.875,'e2='+$ strtrim(string(FORMAT='(%"%8.5f")',e2),2),/norm,$ color='000000'x stop 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 endelse endfor endfor stop end