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