#
# IRAF cl script to demonstrate the routines in the imaging polarimetry package
# impol for WFPC2
#
# Author; J. R. Walsh, ST-ECF. jwalsh@eso.org
#
# REQUIREMENTS:
# An image display tool should be active.
# The package stecf.impol must be loaded.
# The package stsdas.graphics needs to be loaded (to run imdisp_pos for getting
# the NDC's of the display image)
#
set imtype="fits"
#
print "                                     "
print "    *** Start of impol WFPC2 demo *** "
print "    [stecf.impol and stsdas.graphics need to be loaded] "
print "                                     "
#
# (First delete any model reflection nebulae left over from previous demos)
#
print "    Deleting any old demo files "
imdel reflneb1*.fits verify=yes
#
# An input (intensity) image is required. The image reflneb.fits is a 512x512
# real image created by artdata.mkimage with a point source and two elliptical
# nebulae
#
print "    Displaying reflneb.fits on the image display "
display image="reflneb.fits" frame=1 z1=0.0 z2=500.
# 
# Using polimodel a linear polarization (%) and a position angle image will
# be created from this intensity image. The parameters controlling the 
# variation of polarization and position angle are given by the STSDAS 
# table file reflneb.tab. Running trpint on this table produces the 
# following:
#
# row CEN_POL   CEN_X   CEN_Y     EXP   SCALE
#           %  pixels  pixels          pixels
#
#    1   1.000 310.000 290.000   1.000   10.00
#
# The polarization map is simple: it is a centred at pixel 310,290 (i.e. 
# the point source) and increases outward as (radial distance/10). The
# polarization will not be calculated for pixel values less than
# 3.0 and no random errors are applied. The polarization image is
# reflneb1po.fits (in percent) and the position angle image (in degrees) 
# is reflneb1pa.fits
#
print "    Running polimodel to create model images "
polimodel inint="reflneb.fits" simlis="reflneb.tab" coincid=1 intclip=3.0 \
pacorr=0. random=no seed=0 polran=0. paran=0. modpol="reflneb1po.fits" \
modpa="reflneb1pa.fits"
#
# Also needed is an image of the circular polarization. Simply take
# this as a fractional value of the linear polarization image
#
imarith operand1="reflneb1po.fits" op=* operand2=1.0E-4 result="reflneb1cir.fits"
#
# Form the set of input images for the positions of the polarizer. Let
# the HST instrument be WFPC2 with the F502N filter. The position angle of 
# the telescope is 0degrees. An unpolarized background is added to simulate 
# dark current. The electron to ADU conversion factor is assumed to be 7.1
# and the readout noise 5.4 electrons, appropriate for Gain 7. The STSDAS 
# table file wfpc2-pol-filt-001.tab lists the parameters of the WFPC2 
# polarizer filters for all colour filters.
#  
print "    Running hstpolsim to create observed WFPC2 images "
hstpolsim inint="reflneb.fits" inpol="reflneb1po.fits" inpopa="reflneb1pa.fits" \
incir="reflneb1cir.fits" instr="WFPC2" polname="" filtname="F502N" PA_V3=0.0 \
instpolt="wfpc2-pol-filt-001.tab" unback=5.0 e_adu=7.1 read_noi=5.4 seed=3 outroot="reflneb1f"
#
# Make list files of the signal (postfix int) and error (postfix err)
# files.
# First delete any previous list files
#
! rm reflint*.lis
! rm reflerr*.lis
ls -1 reflneb1fint*.fits > reflintlist.lis
ls -1 reflneb1ferr*.fits > reflerrlist.lis
#
# Combine the seven detected signal and their error files and calculate 
# the total signal, linear polarization and position angle. Bin the output
# data 2x2 input pixels and only output pixel values with polarization 
# errors less than 3%
# 
print "    Running hstpolima to calculate polarization images.... "
hstpolima imalis="@reflintlist.lis" errlis="@reflerrlist.lis" \
instpolt="wfpc2-pol-filt-001.tab" xybin=2 binrej=10.0 errlim=3.0 outroot="reflneb1f"
#
# Overplot the derived polarization vector map on the image already displayed
#
imdisp_pos image="reflneb.fits"
#
# [The values of the parameters left right bottom top may have to be adjusted
# to match the polarization vector map on the intensity image.]
#
polimplot.left=imdisp_pos.left
polimplot.right=imdisp_pos.right
polimplot.top=imdisp_pos.top
polimplot.bottom=imdisp_pos.bottom
#
print "    Overplot the polarization vector map on the image "
polimplot inpol="reflneb1fpol.fits" inpopa="reflneb1fpa.fits" poi_err=no inpoer="" inpaer="" \
xybin=1 poscal=0.80 polow=1.0 pohigh=100.0 pacorr=0.0 xysam=5 polab=no title="" xlabel="" \
ylabel="" device="imdw" erase=yes
#
# Also plot a labelled version on the soft device
#
print "    Plot a labelled polarization vector map on stdgraph "
polimplot inpol="reflneb1fpol.fits" inpopa="reflneb1fpa.fits" poi_err=no inpoer="" inpaer="" \
xybin=1 poscal=0.80 polow=1.0 pohigh=100.0 pacorr=0.0 xysam=3 polab=yes title="Model-Nebula" \
xlabel="X" ylabel="Y" device="stdgraph" left=0.1 right=0.9 bottom=0.1 top=0.9 erase=yes
#
# Finally determine the polarization data in three circular apertures centred
# on the point source and the two nebular features. The apertures are listed in 
# STSDAS table file reflnebaper3.tab and are as follows:
#
# row OBJECT_NAME        POS_X   POS_Y  RADIUS     GAP ANNULUS
#                                                             
#
#   1 STAR              310.00  290.00    5.00    0.00    0.00
#   2 REFLNEB_POS_1     200.00  200.00   30.00    0.00    0.00
#   3 REFLNEB_POS_2     300.00  400.00   20.00    0.00    0.00
#
! rm reflneb1aper.tab
print "    Finally determine the polarization in some apertures "
hstpolpoints imalis="@reflintlist.lis" errlis="@reflerrlist.lis" \
instpolt="wfpc2-pol-filt-001.tab" pointab="reflnebaper3.tab" iterej=0 \
sigsky=5.0 outtab="reflneb1aper.tab"
#
# Print out the table of aperture polarization results
#
tprint table="reflneb1aper.tab" prparam=yes prdata=yes
#
# End of demo
#
print "                                     "
print "    *** End of impol WFPC2 demo *** "
print "                                     "
#