#!/bin/csh
alias echo "echo > /dev/null"
kappa
figaro
unalias echo
gdset xw
gdclear
lutheat
#
echo " "
echo "************************************************************** "
echo " "
echo "         EXTRACT IMAGES/SPECTRA FROM SOURCE IFU CUBE   "
echo "             (I.E. OUTPUT FROM IFU_CALIB.PRG)          "
echo " "
echo "                  CJD - 9 Dec 2003                             "
echo " "
echo " 1. Extract images over given wavelenth                        "
echo "    (look at a scrunched spectral image to aid in step 1)      "
echo " 2. Extract spectra over given area                            "
echo "    (use displayed image to aid in step 2)                     "
echo " "
echo "        *** NEED ONLY 1 FRAME:                                 "
echo "        ***          DATA CUBE FROM ifu_calib.prg              "
echo " "
echo "        BUT -  frame name must be name_cube.sdf                "
echo " "
echo " Note:  ifu_cal.prg  uses the cube from oracdr, so can use the "
echo " scrunched spec image, white light image, or an orac-extracted "
echo " image to choose spatial and/or spectral ranges for extraction " 
echo " with this script.                                             " 
echo " "
echo "************************************************************** "
echo " "


##
# History
# -------
#
# Aug 04 - change pixel-scale calibration.

##
#
###################################################
#
# 0.  Input file
#
###################################################
#

ls *cub* *gu*
                                                                                
echo " "
echo "         0. SELECT CUBE  "
echo " "
                                                                                
echo    " Name of target DATA CUBE (from IFU_CALIB.PRG)  "
echo -n " DON'T append the _cube.sdf to the name       : "
  set f1="$<"  
                                 
echo    " Name of target SCRUNCHED IMAGE                 "
echo -n " Can use this to get wavelength ranges        : "
  set f2="$<"  
        
echo " "
echo " >>> Displaying scrunched image in gaia... "
gaia $f2 &
                         
##
#
###################################################
#
# 1.  Extract images and display, one at a time
#
###################################################
#

echo " "
echo " ********************************************"
echo "         1. EXTRACT IMAGES in XY SPACE  "
echo " ********************************************"
echo " "

echo -n " >>> Want to extract an IMAGE in XY plane (y/n)? "
set answr="$<"

  if ($answr == "Y" || $answr == "y") then    
       
  loop1:                   
  echo -n " >>> Ok, file extension (e.g. 10s1, feii, etc.): "
  set name="$<"
 
  echo " "
  echo " Need to select wavelength range across line of interest"
  echo " Examine the scrunched file from oracdr.       "
  echo " Probably need to zoom-in in Gaia to 7x or so. "
  echo " Enter wavelengths as TSTART and TEND below... "
  echo " "
  echo " Since the spectral resolution is 2-pixels with the IFU, "
  echo " a wavelength range over 4 pixels is probably appropriate. "
  echo " "

  xyplane cube=${f1}_cube image=temp
  cdiv in=temp scalar=0.0288 out=${f1}_${name}
  rm temp.sdf

echo " Have DIVIDED by 0.0288 arcsec2, so that flux in final image"
echo " is in Jy/arcsec. "
echo " NB: Multiply by spectral (frequency) resolution, assuming "
echo " R~1000, to get surface brightness in 10{-26}W/m2/arcsec2."
echo " "
echo " Displaying extracted image " ${f1}_${name} " in Gaia "
  gaiadisp ${f1}_${name}
   
  echo " "
  echo -n " >>> Want to subtract adjacent continuum (y/n)? " 
    set answr="$<"
      if ($answr == "Y" || $answr == "y") then    

         echo " OK, need to define sky regions on either side of line "
         echo " OVER THE SAME WAVELENGTH RANGES.  Images will be AVERAGED "
         echo " and subtracted from the on-line image. "
         echo " "
        
           echo "LEFT SKY..." 
           xyplane cube=${f1}_cube image=sky1
           echo "RIGHT SKY..."
           xyplane cube=${f1}_cube image=sky2

           echo " "
           echo " To examine the two sky frames, look at sky1 and sky2. "

           add in1=sky1 in2=sky2 out=skyadd
           cdiv in=skyadd scalar=0.0576 out=sky    ## div by 2x0.0288 = 0.0576
           sub in1=${f1}_${name} in2=sky out=${f1}_${name}_ss

	   echo ""
	   echo " Blinking images before and after continuum-subtraction... "
           echo " Compare counts before and after the continuum subtraction ... "
             set n = 1
                while ($n <= 5)
                gaiadisp ${f1}_${name}
                sleep 1
                gaiadisp ${f1}_${name}_ss
                @ n = $n + 1
             end
 
           rm sky1.sdf sky2.sdf skyadd.sdf sky.sdf
           
      else

         echo " OK, no adjacent continuum subtraction "

      endif

    echo ""
    echo " .....................................................  "
    echo " Stage 1, results:    "
    echo " Extracted image                - " ${f1}_${name}.sdf  
    echo " Extracted cont subt. (if done) - " ${f1}_${name}_ss.sdf
    echo "     - Surface brightness (in both) in Jy/arcsec2      "
    echo " ......................................................"
    echo ""

      echo -n " >>> Extract another IMAGE (y/n)? "
      set answr="$<"
        if ($answr == "Y" || $answr == "y") then    
           goto loop1
        endif

  endif 


##
#
###################################################
#
# 2.  Extract spectra and display, one at a time
#
###################################################
#
   
echo " "                                                                             
echo " *************************************************** "
echo "      2. Extract SPECTRA along T (wavelength) axis  "
echo " *************************************************** "
echo " " 

echo -n " >>> Want to extract a SPECTRUM along T axis (y/n)? "
set answr="$<"

  if ($answr == "Y" || $answr == "y") then    
                              
  loop2:
  echo -n " >>> Ok, subscript for spectrum (e.g. knota, posn1, etc.): "
  set name2="$<"
 
  echo " "
  echo " First - extract plane in x-t by specifying YSTART/YEND "
  echo " with the figaro routine xtplane.  "
  echo " Second - ystract a spectrum by specifying XSTART/XEND "
  echo " "
  echo " Note:  X-axis and Y-axis ranges can be taken from white-light "
  echo " or extracted image (from above) - use X,Y coords in gaia "
  echo " " 
  echo " IMPORTANT: pixel 3,3 is actual 2.5,2.5! " 
  echo " " 
  echo " Note also; spectra seem to be added, so flux scale still in Jy, "
  echo " with area on source specified by area in pixels.       "
  echo " "

  echo "First extract xt plane;  specify Y-range from gaia..."
  xtplane cube=${f1}_cube image=tempplane

  echo " "
  echo "Now ystract spectrum;  specify X-range from gaia..."
  ystract image=tempplane spectrum=${f1}_${name2}_spec

  rm tempplane.sdf

  echo " Plotting extracted spectrum " ${f1}_${name2}_spec " in GWM display "
  splot spectrum=${f1}_${name2}_spec whole=true autoscale=true hardcopy=false \\

  echo ""
  echo " .......................................................  "
  echo " Stage 2, results:    "
  echo " Extracted spectrum             - " ${f1}_${name2}_spec.sdf
  echo "       - Flux density in Jy      "
  echo " ........................................................"
  echo ""

    echo -n " >>> Extract another SPECTRUM along T axis (y/n)? "
    set answr="$<"
       if ($answr == "Y" || $answr == "y") then    
          goto loop2
       endif

  endif 

echo " "
ls ${f1}*
echo " "