# Example command script to compute I, Q, and U images, and then
# generate contour plot with polarization vectors superimposed.
# 
# J. Biretta -- 22 Dec. 97
#
# Uncomment imdel and !rm to delete previous results and re-run.
#
# This script will plot polarization vectors 10 pixels apart. 
# If you want a different spacing, you will need to change
# numbers at points labeled "vector spacing" in script.
# 

stsdas
tools
imgtools
# Compute I, Q, U images.
#imdel i_pol
#imdel q_pol
#imdel u_pol
# These equations are for example of R Mon observed in F555W 
# with POLQ+ WF2, WF3, and WF4.  The input images wf2r, wf3r, wf4r 
# have been wmosaic'ed, aligned, and rotated north-up.  
# Parameters are results from new pol calib tool 18dec97.  
imcalc wf2r,wf3r,wf4r i_pol "1.375*im1 + 0.098*im2 + 1.375*im3"
imcalc wf2r,wf3r,wf4r q_pol "2.157*im1 - 2.695*im2 + 0.459*im3"
imcalc wf2r,wf3r,wf4r u_pol "0.760*im1 + 1.525*im2 - 2.241*im3"

# The rest of the script from here onwards is fairly generic, 
# should work for any polarization image set with only minor 
# adjustments.

# Average images in 10x10 pixel blocks.
# Vector spacing 10 pixels.  For different spacing change 10 and 10.
#imdel i_pol_10
#imdel q_pol_10
#imdel u_pol_10
blkavg i_pol i_pol_10 10 10 option="average"
blkavg q_pol q_pol_10 10 10 option="average"
blkavg u_pol u_pol_10 10 10 option="average"

# Compute fractional pol and direction.
stsdas
toolbox
imgtools
#imdel f_pol
#imdel d_pol
# Compute fractional polarization with clipping.
imcalc i_pol_10,q_pol_10,u_pol_10 f_pol.hhh \
   "if im1.gt.5 then sqrt(im2*im2+im3*im3)/im1 else 0."
# Compute polarization direction.
imcalc i_pol_10,q_pol_10,u_pol_10 d_pol.hhh \
   "if (im2.ge.0.) then 57.3 * atan(im3/im2) / 2.0 + 90.\
   else 57.3 * atan(im3/im2) / 2.0 + 180." 
 
# Generate table needed by fieldplot.
bye
ttools
#imdel pol_frac
#imdel pol_dir
#!rm pol.tab
#!rm pol.dat
# Make table with X, Y, fractional pol.
# wcs logical means use actual pixel numbers from image.
imtab f_pol pol.tab pol_frac pname="x" wcs="logical" 
# Add pol direction to the table.
imtab d_pol pol.tab pol_dir 
# Correct pixel numbers in table for 10x10 blkavg.
# Vector spacing 10 pixels.  If you want different spacing
# change the 10 and 5=10/2 below.
tcalc "pol.tab" "x" "10*(x1-1)+5" datatype=real
tcalc "pol.tab" "y" "10*(x2-1)+5" datatype=real
# Convert iraf table to ascii text table.
tdump pol.tab datafile=pol.dat columns="x y pol_frac pol_dir"

# Generate contour plot.
# Adjust vx1, vx2, vy1, vy2 to get proper aspect ratio.
#!rm pol.meta
contour i_pol floor=0. ceiling=200. zero=0. ncontou=10 \
    interva=0. nhi=-1 dashpat=528 device=stdgraph \
    title="R Mon" preserv+ label- fill- xres=64 yres=64 \
    perimet+ vx1=0.1 vx2=0.9 vy1=0.2 vy2=0.818 subsamp- \
    append- >>G pol.meta

# Overlay vector plot onto contours.
bye
bye
graphics
stplot
fieldplot pol.dat rtheta+ degrees+ magscal=50 crdpsn="center" \
   head- headsiz=0.0075 psnmark=INDEF marksize=0.0075 zeroplo+ \
   axes- margin- title='R Mon' dvpar.device='imdr' \
   dvpar.append+ dvpar.left=0 dvpar.right=0 \
   dvpar.bottom=0 dvpar.top=0 >>G pol.meta

# Convert meta-script to postscript.
#!rm pol.ps
psikern pol.meta devic=psi_port output=pol.ps