;PRO TESTCALIMAGE ; ;PURPOSE: ; This script deals with the ShackHartmann test-calibration image. It ; is called TestCalImage.gs6 and opens with DS9. ; ; The optimal coordinates of the grid and their offsets with respect to ; a 25-pixel-spaced uniform grid centered at the same center are printed ; to a file ; ;CALLING SEQUENCE: ; TestCalImage,[[Image_Dir],Stats_Dir] ; ;KEYWORDS: ; Image_Dir: The directory in which TestCalImage.gs6 is stored ; Stats_Dir: The directory where the coordinates and offset will be ; written ; DATE: The string '10','11',or '12' corresponding to ; 2005-05-10,2005-05-11, and 2005-05-12 ; ;OUTPUS: ; [stats_dir]/Calibrated_grid.txt ; pro test_cal_image_png, Date=Date, Image_Dir=Image_Dir, Stats_Dir=Stats_Dir If (not keyword_set(Date)) then Date='11' If (not keyword_set(x_center)) then x_center=[542] If (not keyword_set(y_center)) then y_center=[465] If (not keyword_set(stats_dir)) then stats_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/2005-05-'+date+'/ShackHartman/stats/' If (not keyword_set(image_dir)) then image_dir='/nfs/slac/g/ki/ki08/lsst/CPanalysis/ShackHartmann_Cal/' raw_file=image_dir+'TestCalImage.fits' coo_file=image_dir+'TestCalImage.fits.coo.1' mag_file=image_dir+'TestCalImage.mag.1.txt' center_file=image_dir+'TestCalImage_centers.txt' stat_file=image_dir+'TestCalImage_stats.txt' fits_read,raw_file,Im get_lun,stat openw,stat,stat_file num=25 subx=2*indgen(num)-25 ;The x-reference numbers suby=indgen(num)-12 ;The y-reference numbers box=11 ;The distance used for finding spot grid_spots=625 fiducial=[1,2,3,4,5,6,7,8,18,19,20,21,22,23,24,25,$ 26,27,28,29,30,31,45,46,47,48,49,50,$ 51,52,53,54,55,72,73,74,75,$ 76,77,78,79,98,99,100,$ 101,102,103,123,124,125,$ 126,127,149,150,$ 151,152,174,175,$ 176,200,$ 201,225,$ 238,239,$ 261,262,263,264,265,266,$ 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,$ 336,337,338,339,340,341,$ 361,362,363,364,365,366,$ 388,389,$ 401,425,$ 426,427,449,450,$ 451,452,452,473,474,475,$ 476,477,478,498,499,500,$ 501,502,503,504,522,523,524,525,$ 526,527,528,529,530,546,547,548,549,550,$ 551,552,553,554,555,556,570,571,572,573,574,575,$ 576,577,578,579,580,581,582,594,595,596,597,598,599,600,$ 601,602,603,604,605,606,607,608,609,617,618,619,620,621,$ 622,623,624,625]-1 ;Figure out the new reference grid, applied around x_center and y_center x_ref_col=12.5*subx+x_center[0] y_ref_col=25.0*suby+y_center[0] ;prepare arrays to be used in the linear fit x_coords=fltarr(num*num) y_coords=fltarr(num*num) Mags=fltarr(num*num) Sharps=fltarr(num*num) Xerrs=fltarr(num*num) Yerrs=fltarr(num*num) x_coords(*)=-1 y_coords(*)=-1 Mags(*)=-1 Sharps(*)=-1 Yerrs(*)=-1 Xerrs(*)=-1 index=findgen(num*num) fullx=fltarr(num*num) fully=fltarr(num*num) fullx(index)=subx(index mod num) fully=rebin(suby,num*num) ; XCENTER, YCENTER, XERR, YERR, MAG readcol,mag_file,x_coord,y_coord,x_err,y_err,$ format='f,f,f,f' ;Read in the data columns from the output IRAF coo file. They are ; SHARP, ROUNDNESS readcol,coo_file,sharp,roundness,$ format='X,X,X,X,f',SKIPLINE=41 ;Each of these vectors will have as many elements as the number of points num_of_iraf_points=N_elements(x_coord) for col_num=0,N_elements(y_ref_col)-1 do begin for row_num=0,N_elements(x_ref_col)-1 do begin ;This index will provide the location on the grid ID_index=col_num*N_elements(x_ref_col)+row_num ;Check to see if this point is in the fiducial region check=where(fiducial eq ID_index) If check(0) ne -1 then begin found_point=0 endif else 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 box) y_index_pos=where(abs(y_coord-y_ref_col(col_num)) lt box) ;Do some indexing gymnastics to determine if they match 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)) ;The grid coordinates have been matched to a spot 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) ;Check to see if it is ;Print this to the stats file printf,stat,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")',$ ID_index,x_coord(x_index),y_coord(y_index),$ x_ref_col(row_num),y_ref_col(col_num),$ roundness(x_index),$ x_ref_dif,y_ref_dif ;found_point=1 ;Add it to the full grid array X_coords(ID_index)=x_coord(x_index) Y_coords(ID_index)=y_coord(y_index) Sharps(ID_index)=sharp(x_index) Xerrs(ID_index)=x_err(x_index) Yerrs(ID_index)=y_err(x_index) found_point=1 break endif else begin found_point=0 endelse endfor endelse ;The grid coordinates had no spot within a box if found_point eq 0 then begin printf,stat,$ FORMAT='(%"%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d")',$ ID_index,-1,-1,-1,-1,-1,-1,-1,-1,-1 X_coords(ID_index)=-1 Y_coords(ID_index)=-1 Mags(ID_index)=-1 Sharps(ID_index)=-1 Xerrs(ID_index)=-1 Yerrs(ID_index)=-1 endif endfor endfor close,stat free_lun,stat ;Now we have an array that mimics the contents of the _stats.txt files ;Fill the Arrays with the processed grids, removing the -1s point_indices=where(x_coords ne -1) valid_x_coords=x_coords(point_indices) valid_y_coords=y_coords(point_indices) xerrs=xerrs(point_indices) yerrs=yerrs(point_indices) fullx=fullx(point_indices) fully=fully(point_indices) ;For some reason, some xerrors and yerrors are zero xzeroes=where(xerrs eq 0,xcount) yzeroes=where(yerrs eq 0,ycount) If xcount ne 0 then xerrs(xzeroes)=0.001 If ycount ne 0 then yerrs(yzeroes)=0.001 x_fit=curvefit(fullx,valid_x_coords,xerrs,x_center,$ /noderivative,function_name='linear_function_x',yerror=xerror) y_fit=curvefit(fully,valid_y_coords,yerrs,y_center,$ /noderivative,function_name='linear_function_y',yerror=yerror) print,FORMAT='(%"%f\t%f")',x_center,y_center ;THE CENTER HAS BEEN FOUND, REESTABLISH THE GRID BASED UPON IT x_coord_cent=fltarr(num*num) y_coord_cent=fltarr(num*num) I11Arr=fltarr(num*num) I22Arr=fltarr(num*num) I12Arr=fltarr(num*num) dArr=fltarr(num*num) eArr=fltarr(num*num) e1Arr=fltarr(num*num) e2Arr=fltarr(num*num) pAArr=fltarr(num*num) x_off_grid=fltarr(num*num) y_off_grid=fltarr(num*num) new_x_ref_col=12.5*subx+x_center new_y_ref_col=25.0*suby+y_center get_lun,offsets openw,offsets,image_dir+'CalEllipticities.txt' HalfSpot=5 for col_num=0,N_elements(new_y_ref_col)-1 do begin for row_num=0,N_elements(new_x_ref_col)-1 do begin ;This index will provide the location on the grid ID_index=col_num*N_elements(x_ref_col)+row_num ;Check to see if this point is in the fiducial region check=where(fiducial eq ID_index) If check(0) eq -1 then begin x_coord_cent(ID_index)=new_x_ref_col(row_num) y_coord_cent(ID_index)=new_y_ref_col(col_num) x_off_grid(ID_index)=x_coords(id_index)-new_x_ref_col(row_num) y_off_grid(ID_index)=y_coords(id_index)-new_y_ref_col(col_num) ;USE FITS FILE TO GET SPOTS Spot=Im(floor(x_coord_cent(ID_index))-HalfSpot:$ floor(x_coord_cent(ID_index))+HalfSpot,$ floor(y_coord_cent(ID_index))-HalfSpot:$ floor(y_coord_cent(ID_index))+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 I11Arr(ID_index)=I11 I22Arr(ID_index)=I22 I12Arr(ID_index)=I12 eArr(ID_index)=e e1Arr(ID_index)=e1 e2Arr(ID_index)=e2 endif else begin x_coord_cent(ID_index)=-1 y_coord_cent(ID_index)=-1 I11=-10 I22=-10 I12=-10 d=-10 e1=-10 e2=-10 e=-10 pa=-10 I22Arr(ID_index)=I22 I12Arr(ID_index)=I12 eArr(ID_index)=e e1Arr(ID_index)=e1 e2Arr(ID_index)=e2 endelse printf,offsets,FORMAT='(%"%f\t%f\t%f\t%f\t%f\t%f\t%f")',$ I11,I22,I12,d,e1,e2,e,pa endfor endfor close,offsets free_lun,offsets set_plot,'ps' PlSpots=where(I11Arr ne -10) device,filename=image_dir+'SpotsCal.ps' loadct,39 ;load color table 39 device,/color ;allow color on the postscript device,ysize=10.5,/inches ;Height of plot in y device,xsize=12.0,/inches ;Width of plot in x device,yoffset=1.0,/inches ;Y position of lower left corner white='FFFFFF'x black='000000'x red='FF0000'x green='00FF00'x blue='FF0000'x !P.CHARSIZE=1.2 !P.THICK=4. !P.multi=[0,1,2] !P.multi=[0,1,3] !p.noerase=1 I11Mean100=[I11Arr,10] I22Mean100=[I22Arr,10] EMean100=[EArr,0.4] PaMean100=[PaArr,!pi/4] PlSpots=[PlSpots,625] X_Coord_cent=[X_Coord_cent,760] Y_Coord_cent=[Y_Coord_cent,760] Partvelvec,I11Arr(PlSpots),fltarr(N_elements(PlSpots)),$ x_Coord_cent(PlSpots),Y_Coord_cent(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=80 !p.noerase=1 Partvelvec,fltarr(N_elements(PlSpots)),I22Arr(PlSpots)/5,$ X_Coord_Cent(PlSpots),Y_Coord_Cent(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=160 Partvelvec,EArr(PlSpots)*cos(PaArr(PlSpots)),$ EArr(PlSpots)*sin(PaArr(PlSpots)),$ X_Coord_cent(PlSpots),Y_Coord_cent(PlSpots),$ yrange=[100,820],xrange=[180,880],$ xstyle=5,ystyle=5,position=[0.025,0.4,0.625,1.0],$ color=240 xyouts,.7,.9,'Total I(x,y)x^2 10 pixels^2',color=80,/norm xyouts,.7,.85,'Total I(x,y)y^2 10 pixels^2',color=160,/norm xyouts,.7,.8,'Ellipticity = 0.4',color=240,/norm device,/close get_lun,offsets openw,offsets,image_dir+'average_offsets.txt' for i=0, grid_spots-1 do begin printf,offsets,FORMAT='(%"%d\t%f\t%f")',i,x_off_grid(i),y_off_grid(i) endfor close,offsets free_lun,offsets set_plot,'z' device, set_font='Courier' ;device, set_resolution=[8000,6000] ;!P.CHARSIZE=10 ;!P.CHARTHICK=15. ;!X.THICK=20. ;!Y.THICK=20. ;!P.THICK=50. ;!P.color=0 device,set_resolution=[800,600] !p.charsize=1 !p.charthick=1.2 !x.thick=2 !y.thick=2 !p.thick=.8 !p.color=0 white='FFFFFF'x black='000000'x red='FF0000'x loadct,39 set_plot,'x' partvelvec,x_off_grid(point_indices),y_off_grid(point_indices),$ x_coord_cent(point_indices),y_coord_cent(point_indices),$ yrange=[150,800],xrange=[200,850] jpgimg=tvrd() tvlct,reds,greens,blues,/get write_png,image_dir+'Test_Cal_Image.png',$ jpgimg,reds,greens,blues stop end