diff options
Diffstat (limited to 'noao/imred/src/dofoe/proc.cl')
| -rw-r--r-- | noao/imred/src/dofoe/proc.cl | 464 |
1 files changed, 464 insertions, 0 deletions
diff --git a/noao/imred/src/dofoe/proc.cl b/noao/imred/src/dofoe/proc.cl new file mode 100644 index 00000000..75740bda --- /dev/null +++ b/noao/imred/src/dofoe/proc.cl @@ -0,0 +1,464 @@ +# PROC -- Process echelle fiber spectra +# This program combines the operations of extraction, flat fielding, and +# dispersion correction in as simple and noninteractive way as possible. +# It supports a second simultaneous arc fiber. The data must all share +# the same position on the 2D image and the same dispersion solution +# apart from small instrumental changes which can be tracked +# automatically. The apertures must be identified sequentially and must +# be properly paired if a arc fiber is used. +# +# If every needed on could add sky subtraction (with a sky fiber) and +# fluxing following the model of the multifiber packages. + +procedure proc (objects, apref, flat, arcs, arctable, naps, objaps, arcaps, + objbeams, arcbeams, fitflat, recenter, scattered, edit, trace, arcap, + clean, dispcor, splot, redo, update, batch, listonly) + +string objects {prompt="List of object spectra"} + +file apref {prompt="Aperture reference spectrum"} +file flat {prompt="Flat field spectrum"} +string arcs {prompt="List of arc spectra"} +file arctable {prompt="Arc assignment table (optional)\n"} + +int naps {prompt="Number of apertures"} +string objaps {prompt="Object apertures"} +string arcaps {prompt="Arc apertures"} +string objbeams {prompt="Object beam numbers"} +string arcbeams {prompt="Arc beam numbers\n"} + +bool fitflat {prompt="Fit and ratio flat field spectrum?"} +bool recenter {prompt="Recenter object apertures?"} +bool scattered {prompt="Subtract scattered light?"} +bool edit {prompt="Edit/review object apertures?"} +bool trace {prompt="Trace object spectra?"} +bool arcap {prompt="Use object apertures for arcs?"} +bool clean {prompt="Detect and replace bad pixels?"} +bool dispcor {prompt="Dispersion correct spectra?"} +bool splot {prompt="Plot the final spectrum?"} +bool redo {prompt="Redo operations if previously done?"} +bool update {prompt="Update spectra if cal data changes?"} +bool batch {prompt="Extract objects in batch?"} +bool listonly {prompt="List steps but don't process?\n"} + +real datamax = INDEF {prompt="Max data value / cosmic ray threshold"} + +bool newaps, newresp, newdisp, newarcs, dobatch + +string anssplot = "yes" {prompt="Splot spectrum?", mode="q", + enum="no|yes|NO|YES"} + +struct *fd1, *fd2 + +begin + string imtype, ectype + string arcref, spec, arc + string arcrefec, specec, arcec, response + string temp1, temp2, done + string str1, objs, arcrefs, log1, log2 + bool reextract, extract, scat, disp, disperr, log + bool splot1, splot2 + int i, j, n, nspec + struct err + + # Call a separate task to do the listing to minimize the size of + # this script and improve it's readability. + + dobatch = no + if (listonly) { + listonly (objects, apref, flat, arcs, scattered, dispcor, + redo, update) + bye + } + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + ectype = ".ec" // imtype + n = strlen (imtype) + + # Get query parameter. + objs = objects + if (arctable == "") + arcrefs = arcs + else + arcrefs = arctable + arcref = "" + + # Temporary files used repeatedly in this script. Under some + # abort circumstances these files may be left behind. + + temp1 = mktemp ("tmp$iraf") + temp2 = mktemp ("tmp$iraf") + done = mktemp ("tmp$iraf") + + # Rather than always have switches on the logfile and verbose flags + # we use TEE and set a file to "dev$null" if output is not desired. + # We must check for the null string to signify no logfile. + + tee.append = yes + if (logfile == "") + log1 = "dev$null" + else + log1 = logfile + if (verbose) + log2 = "STDOUT" + else + log2 = "dev$null" + + # If the update switch is used changes in the calibration data + # can cause images to be reprocessed (if they are in the object + # list). Possible changes are in the aperture definitions, + # response function, dispersion solution, and sensitivity + # function. The newarcs flag is used to only go through the arc + # image headers once setting the reference spectrum, airmass, and + # UT. + + newaps = no + newresp = no + newdisp = no + newarcs = yes + + # Check if there are aperture definitions in the database and + # define them if needed. This is usually somewhat interactive. + # Delete the database entry to start fresh if we enter this + # because of a redo. Set the newaps flag in case an update is + # desired. + + i = strlen (apref) + if (i > n && substr (apref, i-n+1, i) == imtype) + apref = substr (apref, 1, i-n) + + # Initialize + apscript.saturation = INDEF + apscript.references = apref + apscript.profiles = "" + apscript.nfind = naps + apscript.clean = clean + if (splot) { + splot1 = yes + splot2 = yes + } else { + splot1 = no + splot2 = no + } + + reextract = redo + if (reextract || !access (database // "/ap" // apref)) { + if (!access (apref // imtype)) { + printf ("Aperture reference spectrum not found - %s%s\n", + apref, imtype) | scan (err) + error (1, err // "\nCheck setting of imtype") + } + print ("Set reference apertures for ", apref) | tee (log1) + if (access (database // "/ap" // apref)) + delete (database // "/ap" // apref, verify=no) + apscript.ansresize = "yes" + apscript.ansedit = "YES" + apscript.ansfittrace = "yes" + apscript (apref, references="", ansfind="YES", ansrecenter="NO", + anstrace="YES", ansextract="NO") + newaps = yes + } + + if (recenter) + apscript.ansrecenter = "YES" + else + apscript.ansrecenter = "NO" + apscript.ansresize = "NO" + if (edit) + apscript.ansedit = "yes" + else + apscript.ansedit = "NO" + if (trace) + apscript.anstrace = "YES" + else + apscript.anstrace = "NO" + apscript.ansfittrace = "NO" + apscript.ansextract = "YES" + apscript.ansreview = "NO" + + # The next step is to setup the scattered light correction if needed. + # We use the flat field image for the interactive setting unless + # one is not used an then we use the aperture reference. + # If these images have been scattered light corrected we assume the + # scattered light functions parameters are correctly set. + + i = strlen (flat) + if (i > n && substr (flat, i-n+1, i) == imtype) + flat = substr (flat, 1, i-n) + + if (flat != "") + spec = flat + else + spec = apref + + scat = no + if (scattered) { + hselect (spec, "apscatte", yes, > temp1) + fd1 = temp1 + if (fscan (fd1, str1) < 1) + scat = yes + fd1 = ""; delete (temp1, verify=no) + } + if (scat) { + print ("Subtract scattered light in ", spec) | tee (log1) + apscript.ansfitscatter = "yes" + apscript.ansfitsmooth = "yes" + apscript (spec, output="", ansextract="NO", ansscat="YES", + anssmooth="YES") + apscript.ansfitscatter = "NO" + apscript.ansfitsmooth = "NO" + } + + response = "" + if (flat != "") { + response = flat // "norm.ec" + reextract = redo || (update && newaps) + if (reextract || !access (response // imtype) || (update && scat)) { + print ("Create response function ", response) | tee (log1) + + if (access (response // imtype)) + imdelete (response, verify=no) + if (access (flat //ectype)) + imdelete (flat//ectype, verify=no) + + response (flat, apref, response, recenter=recenter, + edit=edit, trace=trace, clean=clean, fitflat=fitflat, + interactive=params.f_interactive, + function=params.f_function, order=params.f_order) + + newresp = yes + } + } + + # If not dispersion correcting we can go directly to extracting + # the object spectra. The reference arcs are the first on + # the arc lists. The processing of the reference arcs is done + # by the task ARCREFS. + + if (dispcor) { + hselect (arcs, "$I", yes, >temp1) + fd1 = temp1 + i = fscan (fd1, arcref) + if (i < 1) + error (1, "No reference arcs") + fd1 = ""; delete (temp1, verify=no) + i = strlen (arcref) + if (i > n && substr (arcref, i-n+1, i) == imtype) + arcref = substr (arcref, 1, i-n) + if (!access (arcref // imtype)) { + printf ("Arc reference spectrum not found - %s%s\n", + arcref, imtype) | scan (err) + error (1, err // "\nCheck setting of imtype") + } + arcrefec = arcref // ectype + reextract = redo || (update && newaps) + if (reextract && access (arcrefec)) + imdelete (arcrefec, verify=no) + + arcrefs (arcref, arcaps, arcbeams, response, done, log1, log2) + } + + # Now we are ready to process the object spectra. + + reextract = redo || (update && (newaps || newresp || newdisp)) + hselect (objs, "$I", yes, > temp1) + fd1 = temp1 + while (fscan (fd1, spec) != EOF) { + i = strlen (spec) + if (i > n && substr (spec, i-n+1, i) == imtype) + spec = substr (spec, 1, i-n) + + # Check if previously done; i.e. arc. + if (access (done)) { + fd2 = done + while (fscan (fd2, specec) != EOF) + if (spec == specec) + break + if (spec == specec) + next + fd2 = "" + } + if (!access (spec // imtype)) { + printf ("Object spectrum not found - %s%s\n", + spec, imtype) | scan (err) + print (err) | tee (log1) + print ("Check setting of imtype") + next + } + specec = spec // ectype + + # Determine required operations from the flags and image header. + scat = no + extract = no + disp = no + if (scattered) { + hselect (spec, "apscatte", yes, > temp2) + fd2 = temp2 + if (fscan (fd2, str1) < 1) + scat = yes + fd2 = ""; delete (temp2, verify=no) + } + if (reextract || !access (specec) || (update && scat)) + extract = yes + else { + hselect (specec, "dc-flag", yes, > temp2) + fd2 = temp2 + if (fscan (fd2, str1) == 1) { + extract = update && newdisp + if (update && !newdisp) + # We really should check if REFSPEC will assign + # different reference spectra. + ; + } else + disp = dispcor + + fd2 = ""; delete (temp2, verify=no) + } + + if (extract) + disp = dispcor + + # If fully processed go to the next object. + if (!extract && !disp) + next + + # If not interactive and the batch flag is set submit rest to batch. + if (batch && !splot1 && !splot2 && apscript.ansedit == "NO") { + fd1 = ""; delete (temp1, verify=no) + flprcache + batch.objects = objs + batch.datamax = datamax + batch.response = response + batch.arcs = arcs + batch.arcref = arcref + batch.arcrefs = arcrefs + batch.objaps = objaps + batch.arcaps = arcaps + batch.objbeams = objbeams + batch.arcbeams = arcbeams + batch.done = done + batch.logfile = log1 + batch.redo = reextract + batch.update = update + batch.scattered = scattered + batch.arcap = arcap + batch.dispcor = dispcor + batch.newaps = newaps + batch.newresp = newresp + batch.newdisp = newdisp + batch.newarcs = newarcs + dobatch = yes + return + } + + # Process the spectrum in foreground. + if (extract) { + if (access (specec)) + imdelete (specec, verify=no) + + if (scat) { + print ("Subtract scattered light in ", spec) | tee (log1) + apscript (spec, output="", ansextract="NO", + ansscat="YES", anssmooth="YES") + } + + print ("Extract object spectrum ", spec) | tee (log1) + setjd (spec, observatory=observatory, date="date-obs", + time="ut", exposure="exptime", jd="jd", hjd="", + ljd="ljd", utdate=yes, uttime=yes, listonly=no, + >> log1) + setairmass (spec, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, update=yes, + override=yes, >> log1) + apscript (spec, saturation=datamax) + if (response != "") + imarith (specec, "/", response, specec) + } + + disperr = no + if (disp) { + # Fix arc headers if necessary. + if (newarcs) { + sections (arcs, option="fullname", >temp2) + setjd ("@"//temp2, observatory=observatory, date="date-obs", + time="ut", exposure="exptime", jd="jd", hjd="", + ljd="ljd", utdate=yes, uttime=yes, listonly=no, + >> log1) + setairmass ("@"//temp2, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, update=yes, + override=yes, >> log1) + delete (temp2, verify=no) + hselect (arcs, "$I", yes, >temp2) + fd2 = temp2 + while (fscan (fd2, arc) != EOF) { + i = strlen (arc) + if (i > n && substr (arc, i-n+1, i) == imtype) + arc = substr (arc, 1, i-n) + hedit (arc, "refspec1", arc, add=yes, verify=no, + show=no, update=yes) + hedit (arc, "arctype", "henear", add=yes, verify=no, + show=no, update=yes) + } + fd2 = ""; delete (temp2, verify=no) + newarcs = no + } + + print ("Assign arc spectra for ", spec) | tee (log1) + refspectra (spec, references=arcrefs, + apertures="", refaps="", ignoreaps=no, + select=params.select, sort=params.sort, + group=params.group, time=params.time, + timewrap=params.timewrap, override=yes, confirm=no, + assign=yes, logfiles="STDOUT", verbose=no) | + tee (log1, > log2) + + doarcs (spec, response, arcref, arcaps, arcbeams, reextract, + arcap, log1, no) + + hselect (specec, "refspec1", yes, > temp2) + fd2 = temp2 + i = fscan (fd2, arc) + fd2 = ""; delete (temp2, verify=no) + if (i < 1) { + print ("No arc reference assigned for ", spec) | tee (log1) + disperr = yes + } else { + print ("Dispersion correct ", spec) | tee (log1) + dispcor (specec, "", linearize=params.linearize, + database=database, table=arcref//ectype, + w1=INDEF, w2=INDEF, dw=INDEF, nw=INDEF, + log=params.log, samedisp=no, flux=params.flux, + global=no, ignoreaps=no, confirm=no, listonly=no, + verbose=verbose, logfile=logfile) + hedit (specec, "dc-flag", 0, add=yes, verify=no, + show=no, update=yes) + } + } + + if (!disperr && (extract || disp)) { + if (splot1) { + print (specec, ":") + str1 = anssplot + if (str1 == "NO" || str1 == "YES") + splot1 = no + if (str1 == "no" || str1 == "NO") + splot2 = no + else + splot2 = yes + } + if (splot2) + splot (specec) + } + + print (spec, >> done) + } + fd1 = ""; delete (temp1, verify=no) + + if (access (done)) + delete (done, verify=no) +end |