#!/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