#!/bin/csh
alias echo "echo > /dev/null"
kappa
figaro
unalias echo
gdset xw
gdclear
lutheat
#
echo " "
echo "************************************************************** "
echo " "
echo "  CREATE/DIVIDE STANDARD CUBE INTO SOURCE CUBE AND FLUX CALIB  "
echo " "
echo "                  CJD - 9 Dec 2003                             "
echo " "
echo " 1. Extracted and coadded (wavelength calibrated) -> stdspec     "
echo "    spectra from standard scrunched spec image.                  "
echo " 2. Blank off any lines in standard with isedit   -> stdspec_ed  "
echo " 3. Divide stdspec_ed by a normalised BB.         -> stdspec_bb  "
echo "    Flux calibrate so std will convert to Jy      -> stdspec_cal "
echo " 4. Expand stdspec into a cube                    -> stdcube     " 
echo " 5. Divide this cube into the cube of source                     " 
echo "    (reduced with EXTENDED_SOURCE_NOSTD/MAP_EXTENDED_SOURCE_NOSTD) "
echo " 6. Scale flux in source cube so in Jy         --->> name_cube   "
echo " "
echo "   *** NEED ONLY 2 FRAMES:                                      "
echo "   ***          SCRUNCHED SPECTRAL IMAGE FROM STANDARD and      "
echo "   ***          DATA CUBE FROM SOURCE                           "
echo " "
echo "************************************************************** "
echo " "

### History
#   -------
# January - automated row extraction so now only enter y-value of
#           brightest row.
#           Prompt for waveband (rather than just assume its HK) 
#           Option to coadd 3, 5 or 7 of extracted standard star spectra.
# June 04 - Correct flux calibration.  Should have been dividing the 
#           standard star spectrum by the BB function, not multiplying!
# Aug 04  - Now handles MAP_EXTENDED_SOURCE frames
#
##
#
###################################################
#
# 0.  Input files and scale so in counts/second
#
###################################################
#

ls gu*
                                                                                
echo " "
echo "         0. INPUT TWO FRAMES AND DIV. BY EXP_TIME  "
echo " "
                                                                                
echo -n "Name of scrunched standard star SPECTRAL IMAGE (gu<utdate>): "
  set f2="$<"
echo -n "Name of target DATA CUBE                    (gu<date>_cube): "
  set f1="$<"  
                                 
## Firstly, make sure both data sets are in counts/second.

  set expos1 = `fitsval $f1 EXP_TIME`
  set expos2 = `fitsval $f2 EXP_TIME`

echo ""
echo " Exposure time on Target   data cube  is " $expos1 "seconds"
echo " Exposure time on Standard spec image is " $expos2 "seconds"
 
echo ""
echo "Dividing by these values to give objcube.sdf and stdspecim.sdf "

cdiv in=$f1 scalar=$expos1 out=objcube
cdiv in=$f2 scalar=$expos2 out=stdspecim
  
gaiadisp stdspecim

##
#
###################################################
#
# 1.  Optimally extract standard spectra and cooadd
#
###################################################
#
                                                                                
echo " "
echo "         1. EXTRACT STANDARD STAR SPECTRA  "
echo " "
                              
echo -n " >>> Extract standard-star spectra and coadd (y/n)? "
set answr1="$<"
  if ($answr1 == "Y" || $answr1 == "y") then

  echo ""
  echo " >>> Extracting only +ve standard spectra       "
  echo " >>> Note: Prompts for centre of BRIGHTEST beam "
  echo "     Will then extract 7 adjacent spectra  " 
  echo "     (3 above and 3 below the brightest...)     "
  echo "     Optimal extraction over 20-pixel range in y. "

  echo "" 
   echo -n " >>> Y-axis value (row) of centre of brightest beam: "
   set beam="$<"

   set a2 = `calc exp="$beam+143+10"`
   set a1 = `calc exp="$beam+143-10"`
   set b2 = `calc exp="$beam+96+10"`
   set b1 = `calc exp="$beam+96-10"`
   set c2 = `calc exp="$beam+47+10"`
   set c1 = `calc exp="$beam+47-10"`
   set d2 = `calc exp="$beam+10"`
   set d1 = `calc exp="$beam-10"`
   set e2 = `calc exp="$beam-46+10"`
   set e1 = `calc exp="$beam-46-10"`
   set f2 = `calc exp="$beam-94+10"`
   set f1 = `calc exp="$beam-94-10"`
   set g2 = `calc exp="$beam-137+10"`
   set g1 = `calc exp="$beam-137-10"`

  echo " "
  echo "Extracting spectra..."
  echo " "

  profile image=stdspecim ystart=$d1 yend=$d2 degree=3 nreject=5 profile=pos4 residual=junk \\
  optextract image=stdspecim profile=pos4 spectrum=stdpos4 \\
    stats stdpos4 > /dev/null
    set max = `parget Maximum stats`
    set hi = `calc exp="($max)*(1.3)"`
    set lo = `calc exp="($max)/(-10)"`
  linplot stdpos4 style="color(curves)=1" ybot=$lo ytop=$hi > /dev/null

  profile image=stdspecim ystart=$c1 yend=$c2 degree=3 nreject=5 profile=pos3 residual=junk \\
  optextract image=stdspecim profile=pos3 spectrum=stdpos3 \\
  linplot stdpos3  style="color(curves)=2" noclear ybot=$lo ytop=$hi > /dev/null

  profile image=stdspecim ystart=$e1 yend=$e2 degree=3 nreject=5 profile=pos5 residual=junk \\
  optextract image=stdspecim profile=pos5 spectrum=stdpos5 \\
  linplot stdpos5  style="color(curves)=2" noclear ybot=$lo ytop=$hi > /dev/null

  profile image=stdspecim ystart=$b1 yend=$b2 degree=3 nreject=5 profile=pos2 residual=junk \\
  optextract image=stdspecim profile=pos2 spectrum=stdpos2 \\
  linplot stdpos2  style="color(curves)=3" noclear ybot=$lo ytop=$hi > /dev/null

  profile image=stdspecim ystart=$f1 yend=$f2 degree=3 nreject=5 profile=pos6 residual=junk \\
  optextract image=stdspecim profile=pos6 spectrum=stdpos6 \\
  linplot stdpos6  style="color(curves)=3" noclear ybot=$lo ytop=$hi > /dev/null

  profile image=stdspecim ystart=$a1 yend=$a2 degree=3 nreject=5 profile=pos1 residual=junk \\
  optextract image=stdspecim profile=pos1 spectrum=stdpos1 \\
  linplot stdpos1  style="color(curves)=4" noclear ybot=$lo ytop=$hi > /dev/null

  profile image=stdspecim ystart=$g1 yend=$g2 degree=3 nreject=5 profile=pos7 residual=junk \\
  optextract image=stdspecim profile=pos7 spectrum=stdpos7 \\
  linplot stdpos7  style="color(curves)=4" noclear ybot=$lo ytop=$hi > /dev/null

  echo ""
  echo " >>> Displaying seven extracted spectra: "
  echo "            Beam 1 (center) in black/white,   rows: " $d1 "to " $d2
  echo "            Beam 2 (up-1)   in red,           rows: " $c1 "to " $c2
  echo "            Beam 4 (dn-1)   in red,           rows: " $e1 "to " $e2
  echo "            Beam 2 (up-2)   in green,         rows: " $b1 "to " $b2
  echo "            Beam 4 (dn-2)   in green,         rows: " $f1 "to " $f2
  echo "            Beam 2 (up-3)   in blue,          rows: " $a1 "to " $a2
  echo "            Beam 4 (dn-3)   in blue,          rows: " $g1 "to " $g2
  echo " "
  echo    " >>> From the displayed Standard Star spectra decide how many to use. "
  echo -n " >>> Coadd the best 3, best 5 or best 7 spectra (3, 5 or 7)?   "
  set howmany="$<"

  if ($howmany == "3") then

   echo " "                                       
   echo " >>> Ok, coadding only the central 3 spectra and displaying..."
   add in1=stdpos3 in2=stdpos4 out=temp1
   add in1=temp1   in2=stdpos5 out=stdspec

  else if ($howmany == "5") then

   echo " "                                       
   echo " >>> Ok, coadding only the central 5 spectra and displaying..."
   add in1=stdpos2 in2=stdpos3 out=temp1
   add in1=stdpos4 in2=stdpos5 out=temp2
   add in1=temp1   in2=temp2   out=add1
   add in1=add1    in2=stdpos6 out=add2
   add in1=add1    in2=add2    out=stdspec
   
  else

   echo " "                                       
   echo " >>> Ok, coadding all 7 spectra and displaying..."
   add in1=stdpos1 in2=stdpos2 out=temp1
   add in1=stdpos3 in2=stdpos4 out=temp2
   add in1=stdpos5 in2=stdpos6 out=temp3
   add in1=temp1   in2=temp2   out=add1
   add in1=stdpos7 in2=temp3   out=add2
   add in1=add1    in2=add2    out=stdspec
 
  endif

  splot spectrum=stdspec whole=true autoscale=true label=stdspec.sdf \\

##
# Tidy up...

  rm stdpos*
  rm add*
  rm temp*
  rm pos*

##
##
  echo ""
  echo " The extracted standard spectrum is called stdspec.sdf"
  echo ""

  else

  echo " ... OK, assuming standard spectrum already exists "
  echo " and is called stdspec.sdf "
  echo " Plotting stdspec.sdf "
  splot spectrum=stdspec whole=true autoscale=true label=stdspec.sdf \\

  endif 

##
#
###################################################
#
# 2.  ISEDIT any lines in the standard
#
###################################################
#
                                                                                
echo " "
echo "         2. PATCH OVER (ISEDIT) LINES IN THE STANDARD  "
echo " "
 
echo -n " >>> Patch over lines in standard; Br lines, say - (y/n)? "
set answr2="$<"
  if ($answr2 == "Y" || $answr2 == "y") then

  echo " ... Plotting stdspec in GWM display... "
  splot spectrum=stdspec whole=true autoscale=true label=stdspec.sdf \\
  cp stdspec.sdf edit.sdf

  loop1:
  echo ""
  echo " Plot section of spectrum to be patched "
  echo " Position cursor on either side of feature and hit J "
  echo " Mark two or three lines, say, then hit R "
  echo " To do K band after H band can repeat (see below) "
  echo ""
  isedit image=edit output=stdspec_ed whole=false 
  
  echo ""
  echo " Plotting whole spectrum after patchwork..."
  splot spectrum=stdspec_ed whole=true autoscale=true label=stdspec_ed.sdf \\
  
  echo -n " >>> Patch over other lines/regions (y/n)? " 
      set answr3="$<"
      if ($answr3 == "Y" || $answr3 == "y") then
        rm edit.sdf
        mv stdspec_ed.sdf edit.sdf
        goto loop1
      endif

##
##
  echo ""
  echo " The patched standard spectrum is called stdspec_ed.sdf"
  echo ""

  rm edit.sdf

  else

    echo -n " >>> Does PATCHED spectrum already exist - (y/n)? "
    set answr4="$<"
    if ($answr4 == "Y" || $answr4 == "y") then
       echo " ... Ok, assuming its called stdspec_ed.sdf "
       echo " Plotting stdspec_ed..."
       splot spectrum=stdspec_ed whole=true autoscale=true label=stdspec_ed.sdf \\
     else    
       echo " ... OK, no lines patched "
       echo " Will copy stdspec.sdf to stdspec_ed.sdf for next step "
       echo " Plotting stdspec_ed..."
       cp stdspec.sdf stdspec_ed.sdf
       splot spectrum=stdspec_ed whole=true autoscale=true label=stdspec_ed.sdf \\
 
    endif
 
  endif

##
#
###################################################
#
# 3. Divide standard star spectrum by normalised BB
#    and then scale so in Jy/pixel.
#
###################################################
#
                                                                                
echo " "
echo "         3. FLUX CALIBRATE STANDARD SPEC. "
echo " "
 
  echo " First divide standard spec by BB normalised at reference wavelength"
  echo " Will then scale spectrum so that division by standard cube "
  echo " will give a source cube in Jy. " 
  echo ""
  echo " NOTE: Script may need tweeking for short-J, short-H and long-K grism! "
  echo " (unless SED of standard is very flat;  A or F-type) "
  echo ""

  echo -n " >>> Waveband of spectra (J,H,K,L,Lp or M) : "
  set wband = "$<" 
  echo -n " >>> Enter magnitude of standard star at $wband : " 
  set mag = "$<" 
  echo -n " >>> Enter Temp of standard star           : "
  set temp = "$<" 

# Create a BB spectrum with same axis data as stdspec_ed

  bbody in=stdspec_ed out=bbody TEMP=$temp
   
  echo ""
  echo " ... Plotting black-body function"
  splot spectrum=bbody whole=true autoscale=true label=BB-function \\

## Note - change wavelength ranges if short-J, short-H and long-K grism

  echo ""
  if      ($wband == "J" || $wband == "j") then 
     echo " Calculating signal at 1.25 microns (note the MEAN)"
     istat image=bbody xstart=1.24 xend=1.26
  else if ($wband == "H" || $wband == "h") then 
     echo " Calculating signal at 1.65 microns (note the MEAN)"
     istat image=bbody xstart=1.64 xend=1.66
  else if ($wband == "K" || $wband == "k") then 
     echo " Calculating signal at 2.20 microns (note the MEAN)"
     istat image=bbody xstart=2.19 xend=2.21
   else if ($wband == "L" || $wband == "l") then 
     echo " Calculating signal at 3.40 microns (note the MEAN)"
     istat image=bbody xstart=3.39 xend=3.41
   else if ($wband == "Lp" || $wband == "lp") then 
     echo " Calculating signal at 3.75 microns (note the MEAN)"
     istat image=bbody xstart=3.74 xend=3.76
   else if ($wband == "M" || $wband == "m") then 
     echo " Calculating signal at 4.80 microns (note the MEAN)"
     istat image=bbody xstart=4.79 xend=4.81
   else 
     echo " Confused about waveband - assuming K ... "
     istat image=bbody xstart=2.19 xend=2.21
   endif 

  echo ""
  echo " >>> Must normalise the BB curve  "
  echo -n " >>> What's the mean (above) at the ref. wavelength? "
  set bbmean="$<"

  echo ""
  echo " Normalising the BB function at the reference wavelength"
  cdiv in=bbody scalar=$bbmean out=bbody_norm
  echo " ... and plotting the normalised BB function"
  splot spectrum=bbody_norm whole=true autoscale=true label=Normal-BB \\

  echo " "
  echo "              *** Hit return to continue *** "
  set pause="$<"

  echo " Dividing standard spec. by normalised BB function"
  echo " Plotting BEFORE (blue) and AFTER (white) division by BB "
  div in1=stdspec_ed in2=bbody_norm out=stdspec_bb 
  linplot stdspec_ed  style="color(curves)=5" clear  > /dev/null
   sleep 1
  linplot stdspec_bb  style="color(curves)=1" noclear  > /dev/null
   sleep 1
  rm bbody*
  
##
##
  echo ""
  echo " The standard spectrum after BB division is called stdspec_bb.sdf "
  echo " Must now flux-scale... "

# Now do maths in 1-D spectrum so that when divide standard
# cube into science target cube, science target data in Jy.

## DATA already scaled so in counts/second as first step (above)

   set factor = `calc exp="10**($mag/(-2.5))"`

  if      ($wband == "J" || $wband == "j") then 
        set flux = `calc exp="(1600*$factor)"`
  else if ($wband == "H" || $wband == "h") then 
        set flux = `calc exp="(1020*$factor)"`
  else if ($wband == "K" || $wband == "k") then 
        set flux = `calc exp="(657*$factor)"`   
   else if ($wband == "L" || $wband == "l") then 
        set flux = `calc exp="(290*$factor)"`
   else if ($wband == "Lp" || $wband == "lp") then 
        set flux = `calc exp="(252*$factor)"`
   else if ($wband == "M" || $wband == "m") then 
        set flux = `calc exp="(163*$factor)"`
   else 
        echo " Confused about waveband - assuming K ... "
        set flux = `calc exp="(657*$factor)"`
   endif 


   echo " "
   echo " Flux at reference wavelength is (Zeroth Flux * 10^(-mag/2.5) =  " $flux "Jy " 
   
   echo " "
   echo " DIVIDING standard star spectrum by " $flux " so that"
   echo " division of target data by standard star spectrum will give "
   echo " cube in Jy/pixel."

   cdiv in=stdspec_bb scalar=$flux out=stdspec_cal
#   splot spectrum=stdspec_cal whole=true autoscale=true erase=true \\

  echo ""
  echo " The standard spectrum after flux calibration is called stdspec_cal.sdf "
  echo ""

  echo "              *** Hit return to continue *** "
  set pause="$<"
 

##
#
###################################################
#
# 4. Create data cube from standard star spectrum
#
###################################################
#
                                                                                
echo " "
echo "         4. CREATE STANDARD STAR DATA CUBE "
echo " "

echo " >>> What's are the Spatial Dimensions of Target DATA CUBE ?"
echo " "
echo "     Normal IFU dimensions are 14 slices (in x) by 54 pixels (in y)"
echo "     However, if you have mapped or jittered (using "
echo "     MAP_EXTENDED_SOURCE), your frame will be larger. "
echo "     The standard star cube must be the same size on the sky, so..."

echo " "
echo -n " Have you JITTERED or MAPPED an extended source (y/n)? "
   set answr="$<"
     if ($answr == "Y" || $answr == "y") then
       echo " Ok, using HDSTRACE to display cube dimensions for the science target"
       echo "  -->> note the x and y size (under DATA_ARRAY)! "
       echo " "
       hdstrace object=objcube | grep DATA
       echo " "
       echo -n " Enter the x and y dimensions of the cube from above (x y): "
       set a = ( $< )
       growx spectrum=stdspec_cal new=true image=tempstd ystart=1 yend=${a[2]} ysize=${a[2]}
       growyt image=tempstd new=true cube=stdcube xstart=1 xend=${a[1]} xsize=${a[1]}
     else
       echo " Ok, assuming cube dimensions are 14 by 54 pixels "
       growx spectrum=stdspec_cal new=true image=tempstd ystart=1 yend=54 ysize=54
       growyt image=tempstd new=true cube=stdcube xstart=1 xend=14 xsize=14
     endif
 
rm tempstd.sdf
                                                                      
  echo ""                                               
  echo " The standard spec. has been grown into cube stdcube.sdf  "
  echo ""                                               
 
  echo "              *** Hit return to continue *** "
  set pause="$<"
 
                                                                              
##
#
###################################################
#
# 5. Divide standard cube into source data cube
#
###################################################
#
                                                                                
last:
echo " "
echo "         5. DIVIDE SOURCE CUBE BY STAND CUBE"
echo " "

echo " Target was called ${f1} ... "
echo " "
echo -n " >>> Name of science target (for the file name): "
set name="$<"

echo " "
echo " Dividing source cube by standard cube.  "

div in1=objcube in2=stdcube out=${name}_cube

  echo ""                                                                         
  echo " ................................................................... "
  echo "  "
  echo "  Reduced object cube, in Jy, is called " ${name}_cube.sdf
  echo ""                                                                         
  echo " ................................................................... "
  echo ""                                                                         
 
ls *cube*

echo " "
echo -n " >>> Calibrate another CUBE with the same standard star data (y/n)? "
   set answr="$<"
     if ($answr == "Y" || $answr == "y") then
       echo -n " Name of target DATA CUBE (gu<date>_cube): "
       set f1="$<"  
          set expos1 = `fitsval $f1 EXP_TIME`
          cdiv in=$f1 scalar=$expos1 out=objcube
       goto last
     endif

echo ""