#
# IRAF cl script to demonstrate some of the routines in the imaging polarimetry 
# package impol, simply plotting of a model reflection nebula
#
# 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 demo *** "
print "                                     "
#
# (First delete any model reflection nebulae left over from previous demos)
#
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
#
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
#
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"
#
# Overplot the polarization vector map on the image already displayed
#
imdisp_pos image="reflneb.fits"
lpar imdisp_pos
#
# [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 inpol="reflneb1po.fits" inpopa="reflneb1pa.fits" poi_err=no inpoer="" inpaer="" \
xybin=1 poscal=0.80 polow=1.0 pohigh=100.0 pacorr=0.0 xysam=9 polab=no title="" xlabel="" \
ylabel="" device="imdw" left=0.3402 right=0.6598 bottom=0.3402 top=0.6598 erase=yes
#
# Also plot a labelled version on the soft device
#
polimplot inpol="reflneb1po.fits" inpopa="reflneb1pa.fits" poi_err=no inpoer="" inpaer="" \
xybin=1 poscal=0.80 polow=1.0 pohigh=100.0 pacorr=0.0 xysam=7 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
#
# End of demo
#
print "                                     "
print "    *** End of impol demo *** "
print "                                     "
#