#!/bin/csh
#
# CGS4 version:
# CJD - June 2000
# Tweeked for UIST, file suffix may still be wrong!!
# CJD - August 2002
# Tweeked further for UIST, file suffix should be right now...
# CJD - November 2002
# Changed to process _wce frames from pipeline
# CJD - July 2005
# Converted to use Polpack
# CJD - Aug 2005
echo " "
echo " * * * * * * * * * * * * * * * * * * * * * * * * * * * * * "
echo " * * "
echo " * QUICK-LOOK FOR UIST SPEC-POL * "
echo " * * "
echo " * Assumes frame naming convention is * "
echo " * * "
echo " * u<DATE>_<NNNNN>_wce * "
echo " * * "
echo " * (i.e. uses flat-fielded wce output files from * "
echo " * the pipeline; can edit to reduce QUICK_LOOK raw * "
echo " * data!) * "
echo " * * "
echo " * Note: Processes only sequential list of 8 spec-pol * "
echo " * frames. Assumes one obj-sky pair at each WP * "
echo " * angle. Order of frames MUST be * "
echo " * obj-sky-sky-obj-obj-sky-sky-obj * "
echo " * * "
echo " * Images displayed in GAIA / Spectra in KAPVIEW * "
echo " * * "
echo " * * * * * * * * * * * * * * * * * * * * * * * * * * * * * "
echo " "
###
# Start the FIGARO, KAPPA & TSP routines (quietly)
#
alias echo 'echo > /dev/null'
kappa
figaro
polpack
gdset xw
unalias echo
loop0:
###
# Enter UT Date and number of first frame
#
echo " "
echo -n "UT data of the observations: "
set date="$<"
echo -n "Number of first frame: "
set first="$<"
@ last = $first + 7
echo " "
###
# Read in data frames, setting number to 5 digits
# n - increment to first frame number
# m - counter (1 to 8)
#
set n = $first
set m = 1
while ($m <= 8)
set num=000${n}
if ($n < 10) set num=0000${n}
if ($n > 99) set num=00${n}
if ($n > 999) set num=0${n}
set polframe${m}=u${date}_${num}_wce
echo " Copied " u${date}_${num}_wce "to " polframe${m} "..."
@ n = $n + 1
@ m = $m + 1
end
###
# Echo comment
#
echo " "
echo "--->>> First subtract OBJECT frames from SKY frames"
echo " Will display the four sky-subtracted spectral images in GAIA"
echo " and read/display the WP angle from the fits header "
echo " "
echo " *** MAKE SURE THESE LOOK SENSIBLE *** "
echo " (with 2 +ve and 2 -ve spectra in each)"
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
###
# Subtract skies from obj frames first; display each result
#
sub $polframe1 $polframe2 polframe_0d
echo " "
echo "SUM: polframe1 - polframe2 = polframe_0d "
gaiadisp polframe_0d
fitslist polframe_0d | grep WPLANGLE
sub $polframe4 $polframe3 polframe_45d
echo " "
echo "SUM: polframe4 - polframe3 = polframe_45d "
gaiadisp polframe_45d
fitslist polframe_45d | grep WPLANGLE
sub $polframe5 $polframe6 polframe_22d
echo " "
echo "SUM: polframe5 - polframe6 = polframe_22d "
gaiadisp polframe_22d
fitslist polframe_22d | grep WPLANGLE
sub $polframe8 $polframe7 polframe_67d
echo " "
echo "SUM: polframe8 - polframe7 = polframe_67d "
gaiadisp polframe_67d
fitslist polframe_67d | grep WPLANGLE
echo " "
echo -n " Correct set of frames (y/n)? "
set answer="$<"
if ($answer == "N" || $answer == "n") then
echo " ... Ok, try again from the top"
goto loop0
endif
loop1:
###
# Optimally extract the spectra
echo " "
echo " ... Will assume that the E- and O-beam spectra are roungly in the "
echo " same location on the array for each waveplate angle "
echo " "
echo " --->>> Use the last displayed image to estimate (to within 1-2 pixels) "
echo " the y-axis positions of the spectra "
echo " "
echo -n " Y-axis Centre of TOP (e) +ve beam: "
set tp="$<"
@ pa = $tp + 10
@ pb = $tp - 10
echo -n " Y-axis Centre of TOP (e) -ve beam: "
set tn="$<"
@ pc = $tn + 10
@ pd = $tn - 10
echo -n " Y-axis Centre of BOTTOM (o) +ve beam: "
set bp="$<"
@ pe = $bp + 10
@ pf = $bp - 10
echo -n " Y-axis Centre of BOTTOM (o) -ve beam: "
set bn="$<"
@ pg = $bn + 10
@ ph = $bn - 10
echo " "
echo " --- >>> Will optimally extract the +ve and -ve spectra; coadd these at each "
echo " WP angle for each e- and o-beam. "
echo " Will then display the EIGHT resulting spectra..."
echo " "
echo " NOTE -- Each spectrum should look like a NORMAL raw stellar spectrum "
echo " If it doesn't, maybe the y-axis posn of the spectrum (or "
echo " the parameters used in the extraction) are wrong. "
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
echo " "
echo " *** E-beam (top) spectra at 0/45/22/67 degs *** "
echo " Plotted in red/blue/green/yellow "
echo " "
profile image=polframe_0d degree=3 nreject=5 ystart=$pa yend=$pb profile=pos residual=junk weight=f
optextract image=polframe_0d profile=pos spectrum=bpos weight=false
profile image=polframe_0d degree=3 nreject=5 ystart=$pc yend=$pd profile=neg residual=junk weight=f
optextract image=polframe_0d profile=neg spectrum=bneg weight=false
# Choosing decent plotting range for Y-axis...
isub image=bpos image1=bneg output=spec0e
stats spec0e > /dev/null
set max = `parget Maximum stats`
set hi = `calc exp="$max+($max/10)"`
set lo = `calc exp="($hi)/(-10)"`
splot spectrum=spec0e whole=t autoscale=f high=$hi low=$lo label=E-beam hardcopy=f erase=t colour=red accept
profile image=polframe_45d degree=3 nreject=5 ystart=$pa yend=$pb profile=pos residual=junk weight=f
optextract image=polframe_45d profile=pos spectrum=bpos weight=false
profile image=polframe_45d degree=3 nreject=5 ystart=$pc yend=$pd profile=neg residual=junk weight=f
optextract image=polframe_45d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec45e
splot spectrum=spec45e erase=f colour=blue accept
profile image=polframe_22d degree=3 nreject=5 ystart=$pa yend=$pb profile=pos residual=junk weight=f
optextract image=polframe_22d profile=pos spectrum=bpos weight=false
profile image=polframe_22d degree=3 nreject=5 ystart=$pc yend=$pd profile=neg residual=junk weight=f
optextract image=polframe_22d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec22e
splot spectrum=spec22e erase=f colour=green accept
profile image=polframe_67d degree=3 nreject=5 ystart=$pa yend=$pb profile=pos residual=junk weight=f
optextract image=polframe_67d profile=pos spectrum=bpos weight=false
profile image=polframe_67d degree=3 nreject=5 ystart=$pc yend=$pd profile=neg residual=junk weight=f
optextract image=polframe_67d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec67e
splot spectrum=spec67e erase=f colour=yellow accept
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
#########################
echo " "
echo " *** O-beam (bottom) spectra at 0/45/22/67 degs *** "
echo " Plotted in red/blue/green/yellow "
echo " "
sleep 1
profile image=polframe_0d degree=3 nreject=5 ystart=$pe yend=$pf profile=pos residual=junk weight=f
optextract image=polframe_0d profile=pos spectrum=bpos weight=false
profile image=polframe_0d degree=3 nreject=5 ystart=$pg yend=$ph profile=neg residual=junk weight=f
optextract image=polframe_0d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec0o
splot spectrum=spec0o label=O-beam hardcopy=f erase=t colour=red accept
profile image=polframe_45d degree=3 nreject=5 ystart=$pe yend=$pf profile=pos residual=junk weight=f
optextract image=polframe_45d profile=pos spectrum=bpos weight=false
profile image=polframe_45d degree=3 nreject=5 ystart=$pg yend=$ph profile=neg residual=junk weight=f
optextract image=polframe_45d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec45o
splot spectrum=spec45o erase=f colour=blue accept
profile image=polframe_22d degree=3 nreject=5 ystart=$pe yend=$pf profile=pos residual=junk weight=f
optextract image=polframe_22d profile=pos spectrum=bpos weight=false
profile image=polframe_22d degree=3 nreject=5 ystart=$pg yend=$ph profile=neg residual=junk weight=f
optextract image=polframe_22d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec22o
splot spectrum=spec22o erase=f colour=green accept
profile image=polframe_67d degree=3 nreject=5 ystart=$pe yend=$pf profile=pos residual=junk weight=f
optextract image=polframe_67d profile=pos spectrum=bpos weight=false
profile image=polframe_67d degree=3 nreject=5 ystart=$pg yend=$ph profile=neg residual=junk weight=f
optextract image=polframe_67d profile=neg spectrum=bneg weight=false
isub image=bpos image1=bneg output=spec67o
splot spectrum=spec67o erase=f colour=yellow accept
echo " "
echo " --->>> Listing extracted spectra (spec*)... "
echo " "
ls spec*o.sdf spec*e.sdf
echo " "
echo -n " Redo extraction of spectra? (y/n)? "
set answr="$<"
if ($answr == "Y" || $answr == "y") then
goto loop1
endif
###
# Convert data to 3D for polpack!
echo " "
echo "--->>> Convert data to 3D for polpack"
sleep 2
permaxes in="spec0e(,1,1)" out=spec3d0e perm="[2,3,1]"
permaxes in="spec45e(,1,1)" out=spec3d45e perm="[2,3,1]"
permaxes in="spec22e(,1,1)" out=spec3d22e perm="[2,3,1]"
permaxes in="spec67e(,1,1)" out=spec3d67e perm="[2,3,1]"
permaxes in="spec0o(,1,1)" out=spec3d0o perm="[2,3,1]"
permaxes in="spec45o(,1,1)" out=spec3d45o perm="[2,3,1]"
permaxes in="spec22o(,1,1)" out=spec3d22o perm="[2,3,1]"
permaxes in="spec67o(,1,1)" out=spec3d67o perm="[2,3,1]"
echo " "
echo " --->>> Listing reformatted spectra (spec3d*)... "
echo " "
ls spec3d*o.sdf spec3d*e.sdf
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
echo " "
###
# Preparation for polpack...
#
echo " "
echo " ********************* POLPACK STUFF ************************ "
echo " "
echo " HEADERS "
echo " ------- "
echo " POLPACK needs to know certain items of header information about each"
echo " of the supplied data files. These include the half-wave plate position, "
echo " the orientation of the reference direction, etc. Different instruments may "
echo " store this information in different ways (in different headers). Consequently, "
echo " POLPACK includes an application which finds the required info and stores them "
echo " in a standard place within the data file, known as the POLPACK extension. "
echo " "
echo " The process of copying this information from its original, instrument specific "
echo " location, into the POLPACK extension is known as importing the data, and is "
echo " performed by the POLIMP application. (POLEXT can also be used to explicitly write "
echo " data to the extension; its used below to flag spectra as e- or o-beams, using "
echo " the RAY header.) "
echo " "
echo " However, to do this, POLIMP requires a control table. This is a text file "
echo " containing a line for each item of information to be imported into the POLPACK "
echo " extension. "
echo " "
echo " *** Below, this script should create the appropriate table! *** "
echo " "
echo " An example table is also available from the UKIRT/IRPOL web pages: "
echo " - http://www.jach.hawaii.edu/UKIRT/instruments/irpol/irpol.html "
echo " "
echo " Further info. is also given in the Starlink/Polpack web pages: "
echo " - http://www.starlink.rl.ac.uk/star/docs/sun223.htx/sun223.html"
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
touch polimp_e.table
touch polimp_o.table
echo "WPLATE WPLANGLE 0=0.0,45=45.0,22.5=22.5,67.5=67.5 " >> polimp_e.table
echo "_CHAR IDATE " >> polimp_e.table
echo "_CHAR RUN " >> polimp_e.table
echo "IMGID? RUN//IDATE " >> polimp_e.table
echo "WPLATE WPLANGLE 0=0.0,45=45.0,22.5=22.5,67.5=67.5 " >> polimp_o.table
echo "_CHAR IDATE " >> polimp_o.table
echo "_CHAR RUN " >> polimp_o.table
echo "IMGID? RUN//IDATE " >> polimp_o.table
echo " "
echo "--->>> POLIMP processing of polpack extension for E-beam spectra..."
sleep 3
polimp in='spec3d*e' table=polimp_e.table quiet
polext in='spec3d*e' ray=E
echo " "
echo "--->>> POLIMP processing of polpack extension for O-beam spectra..."
sleep 3
polimp in='spec3d*o' table=polimp_o.table quiet
polext in='spec3d*o' ray=O
echo " "
sleep 2
echo "--->>> Now will use POLCAL to create the polarisation cube "
echo " "
echo " NOTES: "
echo " ------ "
echo " POLCAL converts the eight spectra (4 e-beam and 4 o-beam; one at each WP angle) "
echo " into a cube holding Stokes vectors. Because the input NDFs have to be 3-dimensional "
echo " cubes (converted earlier in this script to 1x1x1024), the output NDF will be "
echo " 4-dimensional."
echo " The axes in the output array have the dimensions 1x1x1024x3, where the final axis "
echo " corresponds to the Stokes parameters I, Q, and U."
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
###
# Now using POLCAL
echo " "
echo " Running POLCAL with params : pmode=linear nodualbeam out=specpol ilevel=3 "
#polcal "spec3d*e,spec3d*o" pmode=linear nodualbeam out=specpol ilevel=3
polcal "spec3d*" pmode=linear nodualbeam out=specpol ilevel=3
echo " ***************************************************** "
echo " OUTPUT DATA CUBE is called specpol.sdf "
echo " * * "
echo " Reminder: frames processed were $first to $last "
echo " ***************************************************** "
echo " "
echo " The bounds of the cube are given below..."
echo " Note that the 1:3 bounds are for I, Q and U. "
echo " "
ndftrace specpol | grep bounds
echo " "
echo -n " (Hit any key to continue...) "
set pause="$<"
# POLVEC to extract data from cube...
#
echo " "
echo " --->>> Finally, extract and plot spectra of polarisation, angle, etc with POLVEC"
polvec in=specpol cat=! P=specpol_sp-P IP=specpol_sp-PI ANG=specpol_sp-TH variance nodebias
polvec in=specpol cat=! I=specpol_sp-I Q=specpol_sp-Q U=specpol_sp-U variance nodebias
echo " "
echo " Extracting the following files from the specpol.sdf datacube: "
echo " "
ls specpol_*
echo " "
echo " WHERE: "
echo " Intensity spectrum - specpol_sp-I "
echo " Q stokes-parameters - specpol_sp-Q "
echo " U stokes-parameters - specpol_sp-U "
echo " Polarisation percentage - specpol_sp-P "
echo " Position Angle - specpol_sp-TH "
echo " Polarised intensity - specpol_sp-PI "
loop2:
echo " "
echo -n " --->>> Plot a spectrum (I, Q, U, P, PI, TH - upper case) : "
set toplot="$<"
ndfcopy specpol_sp-${toplot} temp trim
linplot ndf=temp device=xw
rm temp*
echo " "
echo -n " Another plot (y/n)? "
set answr="$<"
if ($answr == "Y" || $answr == "y") then
goto loop2
endif
###
# Tidy up a bit...
echo " "
echo " Remove temporary frames (polframe*, spec*e, spec*o, polimp*) "
rm polframe*.sdf
rm spec*e.sdf
rm spec*o.sdf
rm polimp*.table
echo " ... DONE!"
echo " "
exit