Example showing derivation of Stokes images and polarization vector plot.
J. Biretta 01/28/1997
The three images of R Mon taken through the POLQ filter in its unrotated position using CCDs WF2, WF3, and WF4. We obtained the uncalibrated .d0h files from the HST data archive, and then calibrated them using the most recent version of CALWP2. We used the F555W+POLQ reference flat ga711093u.r4h, which is also available in the archive.
Data files were as follows:
image filters CCD PA_V3 --------- ---------- --- ----- u2m70401t F555W+POLQ WF2 255.25 u2m70405t F555W+POLQ WF3 255.25 u2m70409t F555W+POLQ WF4 255.25
After calibration, cosmic rays were removed. Since there were no CR-SPLITs performed, this was done with a rather complicated algorithm, which I will not explain.
Next, it was necessary to remove geometric distortion from the images. This was done using "wmosaic" in IRAF stsdas.hst.wfpc. Since the three images were also taken at the same PA_V3, this has the added effect of rotating the images to the same orientation. There is a small redisual distortion remaining from the POLQ filter itself, but we will ignore this.
wmosaic u2m70401t.c0h u2m70401t.hhh wmosaic u2m70405t.c0h u2m70405t.hhh wmosaic u2m70409t.c0h u2m70409t.hhh
At this point the images had the same orientation, but were shifted since they were taken on different CCDs. We used "imcopy" and "imlintran" in IRAF "images" to extract aligned 401x401 pixel subimages from the approx. 1600x1600 pixel output from "wmosaic." Below "xin" and "yin" represent the location of some feature in the inout images, and "xout" and "yout" represent the desired position of this feature in the output image.
Extract roughly aligned 401x401 pixel sub-images:
imcop u2m70401t.hhh[221:621,172:572] wf2 imcop u2m70405t.hhh[964:1364,159:559] wf3 imcop u2m70409t.hhh[959:1359,933:1333] wf4
Now for more exact alignment, we apply "imlintran" to each image:
Image Reduction and Analysis Facility
PACKAGE = images
TASK = imlintran
input = wf3 Input data
output = wf3s Output data
xrotatio= 0. X axis rotation in degrees
yrotatio= 0. Y axis rotation in degrees
xmag = 1. X output pixels per input pixel
ymag = 1. Y output pixels per input pixel
(xin = 119.8) X origin of input frame in pixels
(yin = 105.7) Y origin of input frame in pixels
(xout = 119.2) X origin of output frame in pixels
(yout = 104.7) Y origin of output frame in pixels
(ncols = 401.) Number of columns in the output image
(nlines = 401.) Number of lines in the output image
(interpo= linear) Interpolant (nearest,linear,poly3,poly5,spline3)
(boundar= nearest) Boundary extension (nearest,constant,reflect,wra
(constan= 0.) Constant boundary extension
(fluxcon= yes) Preserve image flux?
(nxblock= 256) X dimension of blocking factor
(nyblock= 256) Y dimension of blocking factor
(mode = ql)
Finally, we used "imlintran" again to rotate the aligned images to north-up orientation. A careful look at figure 3.11 in the WFPC2 Handbook (v.4) will convince you that the rotation needed for a "wmosaic" image is approximately (PA_V3) - 225 degrees (accuracy about 0.2 degrees). Hence the images were rotated as below. The output images are 500x500 pixels. This was applied to all three of the images.
Image Reduction and Analysis Facility
PACKAGE = images
TASK = imlintran
input = wf3s Input data
output = wf3r Output data
xrotatio= 30.25 X axis rotation in degrees
yrotatio= 30.25 Y axis rotation in degrees
xmag = 1. X output pixels per input pixel
ymag = 1. Y output pixels per input pixel
(xin = 200.) X origin of input frame in pixels
(yin = 200.) Y origin of input frame in pixels
(xout = 250.) X origin of output frame in pixels
(yout = 250.) Y origin of output frame in pixels
(ncols = 500.) Number of columns in the output image
(nlines = 500.) Number of lines in the output image
(interpo= linear) Interpolant (nearest,linear,poly3,poly5,spline3)
(boundar= const) Boundary extension (nearest,constant,reflect,wra
(constan= 0.) Constant boundary extension
(fluxcon= yes) Preserve image flux?
(nxblock= 256) X dimension of blocking factor
(nyblock= 256) Y dimension of blocking factor
(mode = ql)
At this point we went to the WFPC2 calibration WWW tool, and generated recipes for I, Q, and U. Inputs were as follows:
We then hit the "calculate" button, and these results come back:
To compute I, Q, and U, we use "imcalc" in stsdas.tools.imgtools along with the result from the WWW tool:
stsdas toolbox imgtools imcalc wf2r,wf3r,wf4r i_pol "1.372*im1+0.104*im2+1.372*im3" imcalc wf2r,wf3r,wf4r q_pol "2.430*im1-2.696*im2+0.182*im3" imcalc wf2r,wf3r,wf4r u_pol "1.248*im1+1.525*im2+-2.726*im3" imcalc i_pol,u_pol,q_pol f_pol "sqrt(im2*im2+im3*im3)/im1"
Finally, after alll this work, we would like to make a nice display of the images. Currently this is best done in the AIPS software package from NRAO. To get the data into AIPS, we need FITS files. Hence write the images out to files with "wfits" in IRAF "dataio" using below command. Note output must have capital letters.
wfits i_pol.hhh IPOL.FITS
Then start AIPS in another window. When it is up, type "free" which will give a list of disk area and available space:
>free AIPS 1: Disk Volume name Total Full Free Timd Access AIPS 1: # blocks % blocks days AIPS 1: 1 /ubach/data1/DA01 963342 68 280125 99.0 Alluser AIPS 1: 2 /ubach/data2/DA01 1759749 57 686336 99.0 Alluser AIPS 1: 3 /ubach/data3/DA01 3940910 77 808735 99.0 Alluser
We see that disk 3 has a bunch of space. We go back to the IRAF window and type:
!cp *_POL.FITS /ubach/data3/DA01
Which copies the files where AIPS can get them. Next we need to read the data into AIPS. For this we use IMLOD in AIPS. Note that we specify "DA03" since the data are on the third AIPS disk.
>task 'imlod' >outnam 'rmon' >outcl 'ipol' >outdisk 3 >infile 'da03:i_pol.fits' >input imlod AIPS 1: IMLOD: Task to store an image from a FITS or IBM-CV tape AIPS 1: Adverbs Values Comments AIPS 1: ---------------------------------------------------------------- AIPS 1: INTAPE 1 Input tape drive # (0 => 1) AIPS 1: OUTNAME 'RMON ' Image name (name) AIPS 1: OUTCLASS 'IPOL ' Image name (class) AIPS 1: OUTSEQ 0 Image name (seq. #) AIPS 1: 0 => highest unique number AIPS 1: -1 => FITS tape value AIPS 1: OUTDISK 3 Disk drive # (0 => any) AIPS 1: NCOUNT 0 Number of files to load. AIPS 1: DOTABLE 1 True (1.0) means load tables AIPS 1: NFILES 0 # of files to advance on tape AIPS 1: NMAPS 1 # IBM maps to advance on tape AIPS 1: INFILE 'DA03:I_POL.FITS Disk file name (FITS only) AIPS 1: AIPS 1: >go imlod
And then repeat this for each of the 4 images. AIPS will later become very unhappy if we do not set the image scale to the proper value. We do this with the following incantation. Note that the pixel scale on X is negative; this is the AIPS convention for RA. Assuming the first image is in AIPS disk 3, catalog number 1:
>indisk 3; getn 1 >keyword='crdelt1'; keyval=-0.09965/3600.,0; puth >keyword='crdelt2'; keyval= 0.09965/3600.,0; puth >imhe AIPS 1: Image=I_POL[1/ (MA) Filename=RMON .IMAP . 1 AIPS 1: Telescope= Receiver=WFPC2 AIPS 1: Observer= User #= 1103 AIPS 1: Observ. date=05-FEB-1995 Map date=25-JAN-1997 AIPS 1: Minimum= 0.00000000E+00 Maximum= 1.08949385E+04 AIPS 1: ---------------------------------------------------------------- AIPS 1: Type Pixels Coord value at Pixel Coord incr Rotat AIPS 1: RA---TAN 500 06 39 09.980 214.00 -0.099650 0.00 AIPS 1: DEC--TAN 500 08 44 23.700 115.00 0.099650 0.00 AIPS 1: ---------------------------------------------------------------- AIPS 1: Coordinate equinox 1950.00 AIPS 1: Maximum version number of extension files of type HI is 1
The last command will display the modified header, so you can check for the correct image scale. Then repeat for the other 3 polarization images.
We need to compute an image of the polarization angle direction. For this we use the COMB task in AIPS. The dialog might look something like below.
>task 'comb' >indisk 3; mcat AIPS 1: Catalog on disk 3 AIPS 1: Cat Usid Mapname Class Seq Pt Last access Stat AIPS 1: 1 1103 RMON .IMAP . 1 MA 27-JAN-1997 19:00:29 AIPS 1: 2 1103 RMON .QMAP . 1 MA 27-JAN-1997 19:23:22 AIPS 1: 3 1103 RMON .UPOL . 2 MA 27-JAN-1997 19:23:22 AIPS 1: 4 1103 RMON .FPOL . 1 MA 27-JAN-1997 19:23:03 [note where your images are] >indisk 1; getn 2; get2 3 >outnam 'rmon'; outcl 'pola' >doalign=-2; opcode 'pola' >blc 0; trc 0; aparm 0; bparm 0 >inp AIPS 1: COMB: Task to combine in many ways two overlapping images AIPS 1: Adverbs Values Comments AIPS 1: ---------------------------------------------------------------- AIPS 1: USERID 0 User ID. 0 => current user, AIPS 1: 32000 => any user. AIPS 1: INNAME 'RMON ' First image name (name) AIPS 1: INCLASS 'QMAP ' First image name (class) AIPS 1: INSEQ 1 First image name (seq. #) AIPS 1: INDISK 3 First image disk drive # AIPS 1: IN2NAME 'RMON ' Second image name (name) AIPS 1: IN2CLASS 'UPOL ' Second image name (class) AIPS 1: IN2SEQ 2 Second image name (seq. #) AIPS 1: IN2DISK 3 Second image disk drive # AIPS 1: DOALIGN -2 Should images be coincident? AIPS 1: (See HELP.) AIPS 1: OUTNAME 'RMON ' Output image name (name) AIPS 1: OUTCLASS ' ' Output image name (class) AIPS 1: OUTSEQ 0 Output image name (seq. #) AIPS 1: OUTDISK 3 Output image disk drive # AIPS 1: BLC *all 0 Bottom left corner AIPS 1: TRC *all 0 Top right corner AIPS 1: OPCODE 'POLA' Algorithm type: AIPS 1: 'SUM ','DIV ','SPIX','POLI', AIPS 1: 'POLA','MULT','OPTD','CLIP' AIPS 1: 'REAL','IMAG','MEAN','RM ' AIPS 1: 'POLC' AIPS 1: APARM *all 0 Parameters for algorithm: AIPS 1: (1) - (4) scale and offset AIPS 1: (8) > 0 => blank with 0.0 AIPS 1: (9) Map1 clip level AIPS 1: (10) Map2 clip level AIPS 1: see HELP COMB AIPS 1: BPARM *all 0 Noise/control parameters: AIPS 1: (1) Map1 noise level AIPS 1: (2) Map2 noise level AIPS 1: (3) > 0 => output noise AIPS 1: (4) < 0.5 => clip w inputs AIPS 1: > 1.5 => clip w S/N AIPS 1: else => clip w noise AIPS 1: (5) minimum ok S/N or AIPS 1: maximum ok noise AIPS 1: (6) max output noise AIPS 1: see HELP COMB > >go comb
At long last, we generate a contour image of the target, with the polarization vectors displayed. We use AIPS task PCNTR. You can get a help file by typing HELP PCNTR. The three input images and other parameters are filled in is below. Stuff in [] are just comments.
>task pcntr [set names of I, F, and PA images] >indisk 3; getn 1; get2 4; get3 6 [plot both contours and vectors] >docont 1; dovect 1 [set desired image region in pixels] >blc 100, 83; trc 301, 239 [select style of plot] >ltype 8 [set contour levels; here we use percent] >plev 1; clev 0; levs 0.5,1,2,3,4,6,8,10,15,20,30,40,60,80 [set length and spacing of vectors] factor 20; xinc 7; yinc 7 [set cutoffs where no vectors are plotted] >pcut 0.003; icut 20 >inp AIPS 1: PCNTR: Task to generate plot file for contour plus pol. vectors AIPS 1: Adverbs Values Comments AIPS 1: ---------------------------------------------------------------- AIPS 1: DOCONT 1 Draw contours? > 0 => yes AIPS 1: DOVECT 1 Draw pol. vectors? > 0 => yes AIPS 1: USERID 0 Image owner ID number AIPS 1: Total intensity image: AIPS 1: INNAME 'RMON ' Image name (name) AIPS 1: INCLASS 'IMAP ' Image name (class) AIPS 1: INSEQ 1 Image name (seq. #) AIPS 1: INDISK 3 Disk unit # AIPS 1: Polarization intensity image: AIPS 1: IN2NAME 'RMON ' (name) blank => INNAME AIPS 1: IN2CLASS 'FPOL ' (class) blank => 'PPOL' AIPS 1: IN2SEQ 2 (seq. #) 0 => high AIPS 1: IN2DISK 3 Disk drive #, 0 => any AIPS 1: Polarization angle image: AIPS 1: IN3NAME 'RMON ' (name) blank => INNAME AIPS 1: IN3CLASS 'POLA ' (class) blank => 'PANG' AIPS 1: IN3SEQ 5 (seq. #) 0 => high AIPS 1: IN3DISK 3 Disk drive #, 0 => any AIPS 1: BLC 100 83 Bottom left corner of images AIPS 1: 1 1 1 1 AIPS 1: 1 AIPS 1: TRC 301 239 Top right corner of images AIPS 1: 1 1 1 1 AIPS 1: 1 AIPS 1: XYRATIO 0 X to Y axis plot ratio. 0=> AIPS 1: header inc or window ratio. AIPS 1: LTYPE 8 Type of labeling: 1 border, AIPS 1: 2 no ticks, 3 standard, 4 rel AIPS 1: to center, 5 rel to subim cen AIPS 1: 6 pixels, 7-10 as 3-6 with AIPS 1: only tick labels AIPS 1: <0 -> no date/time AIPS 1: PLEV 1 Percent of peak for levs. AIPS 1: CLEV 0 Absolute value for levs AIPS 1: (used only if PLEV = 0). AIPS 1: CLEV=PLEV=0 => PLEV=10 AIPS 1: LEVS 0.5 1 Contour levels (up to 30). AIPS 1: 2 3 4 6 AIPS 1: 8 10 15 20 AIPS 1: 30 40 60 80 AIPS 1: *rest 0 AIPS 1: FACTOR 20 Mult. factor for Pol vector AIPS 1: (see HELP) AIPS 1: XINC 7 X-inc. of Pol vectors. 0=>1 AIPS 1: YINC 7 Y-inc. of Pol vectors. 0=>1 AIPS 1: PCUT 0.003 Pol. vector cutoff. P units. AIPS 1: ICUT 30 Int. vector cutoff. I units. AIPS 1: DOALIGN -2 Maps must align? > 0 => yes AIPS 1: See HELP DOALIGN! AIPS 1: DOCIRCLE -1 > 0 => extend ticks to form AIPS 1: coordinate grid AIPS 1: INVERS 0 STar file version number. AIPS 1: STFACTOR 0 Scale star sizes: 0 => none. AIPS 1: > 0 crosses with no labels AIPS 1: < 0 crosses with labels AIPS 1: CBPLOT -1 Position for beam plot: AIPS 1: -1: don't plot beam AIPS 1: 1: lower left (default) AIPS 1: 2: lower right AIPS 1: 3: upper right AIPS 1: 4: upper left AIPS 1: 6-9 : fill in a little AIPS 1: 11-14: more filled AIPS 1: 16-19: even more AIPS 1: DOTV -1 > 0 Do plot on the TV, else AIPS 1: make a plot file AIPS 1: GRCHAN 0 Graphics channel 0 => 1. AIPS 1: TVCORN 0 0 TV pixel location of bottom AIPS 1: left corner of image 0=> self AIPS 1: scale, non 0 => pixel scale. >go pcntr
This will generate either a TV display (input DOTV=1 above) or a disk file which can be used to generate hardcopy (DOTV=-1). To actually get a hardcopy, you also need to run AIPS task LWPLA; this will send the plot to a printer or make a postscript file on disk.
We later used "xv" and "xfig" to generate the two-panel figure with the gray scale image and contour/vector plot.
Let me know if you have problems with this. Good luck. - John B. (biretta@stsci.edu)