#
# IRAF cl script to demonstrate the routines in the imaging polarimetry package
# impol for ACS
#
# 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 ACS 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.
#
# Create table of reflection nebula parameters from ascii file
#
! rm reflneb.tab
tcreate table="reflneb.tab" cdfile="model.col" datafile="reflneb.dat" nskip=0
# 
# 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 ACS with the Wide Field Channel (WFC) Visible 
# Polarizer POL0V, POL60V and POL120V. The selected filter is F55W. The 
# position angle of the telescope is 30degrees. A small unpolarized 
# background is added to simulate dark current. The electron to ADU 
# conversion factor is 4.0 and the readout noise 4.0 electrons. The 
# STSDAS table file acs-pol-wfc-v-filt-001.tab lists the parameters of 
# the ten ACS polarizer-colour filter combinations.
#  
#  
# Create table of polarizing filter parameters from ascii file
#
! rm acs-pol-wfc-v-filt-001.tab
tcreate table="acs-pol-wfc-v-filt-001.tab" cdfile="define.col" \
datafile="acs-pol-wfc-v-filt-001.dat" nskip=3
#
print "    Running hstpolsim to create observed ACS images "
hstpolsim inint="reflneb.fits" inpol="reflneb1po.fits" inpopa="reflneb1pa.fits" \
incir="reflneb1cir.fits" instr="ACS" polname="POL0V" filtname="F555W" PA_V3=30.0 \
instpolt="acs-pol-wfc-v-filt-001.tab" unback=1.0 e_adu=4.0 read_noi=4.0 seed=7 \
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 three 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="acs-pol-wfc-v-filt-001.tab" xybin=3 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
#
# Create table of photometric aperture parameters from ascii file
#
! rm reflnebaper3.tab
tcreate table="reflnebaper3.tab" cdfile="aperture.col" \
datafile="reflnebaper3.dat" nskip=0
#
# 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="acs-pol-wfc-v-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 ACS demo *** "
print "                                     "
#