diff options
Diffstat (limited to 'noao/imred/src/fibers/proc.cl')
| -rw-r--r-- | noao/imred/src/fibers/proc.cl | 707 |
1 files changed, 707 insertions, 0 deletions
diff --git a/noao/imred/src/fibers/proc.cl b/noao/imred/src/fibers/proc.cl new file mode 100644 index 00000000..a55039d8 --- /dev/null +++ b/noao/imred/src/fibers/proc.cl @@ -0,0 +1,707 @@ +# PROC -- Process spectra from 2D to wavelength calibrated 1D. +# This program combines the operations of extraction, flat fielding, +# fiber throughput correction, dispersion correction, and sky subtraction +# in as simple and noninteractive way as possible. Certain assumptions +# are made about the data and the output. A blank sky image, called a +# sky flat, may be used to determine the instrument throughput. 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. + +procedure proc (objects, apref, flat, throughput, arcs1, arcs2, arcreplace, + arctable, fibers, apidtable, crval, cdelt, objaps, skyaps, arcaps, + objbeams, skybeams, arcbeams, scattered, fitflat, recenter, edit, + trace, arcap, clean, dispcor, savearcs, skyalign, skysubtract, + skyedit, saveskys, splot, redo, update, batch, listonly) + +string objects {prompt="List of object spectra"} + +file apref {prompt="Aperture reference spectrum"} +file flat {prompt="Flat field spectrum"} +file throughput {prompt="Throughput file or image (optional)"} +string arcs1 {prompt="List of arc spectra"} +string arcs2 {prompt="List of shift arc spectra"} +file arcreplace {prompt="Special aperture replacements"} +file arctable {prompt="Arc assignment table (optional)\n"} + +int fibers {prompt="Number of fibers"} +file apidtable {prompt="Aperture identifications"} +string crval = "INDEF" {prompt="Approximate wavelength"} +string cdelt = "INDEF" {prompt="Approximate dispersion"} +string objaps {prompt="Object apertures"} +string skyaps {prompt="Sky apertures"} +string arcaps {prompt="Arc apertures"} +string objbeams {prompt="Object beam numbers"} +string skybeams {prompt="Sky beam numbers"} +string arcbeams {prompt="Arc beam numbers\n"} + +bool scattered {prompt="Subtract scattered light?"} +bool fitflat {prompt="Fit and ratio flat field spectrum?"} +bool recenter {prompt="Recenter object apertures?"} +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 savearcs {prompt="Save internal arcs?"} +bool skyalign {prompt="Align sky lines?"} +bool skysubtract {prompt="Subtract sky?"} +bool skyedit {prompt="Edit the sky spectra?"} +bool saveskys {prompt="Save sky 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"} + +string ansskyedit = "yes" {prompt="Edit the sky spectra?", mode="q"} +bool newaps, newresp, newdisp, newarcs, dobatch + +string anssplot = "yes" {prompt="Splot spectrum?", mode="q", + enum="no|yes|NO|YES"} + +string extn = ".ms" {prompt="Extraction extension"} +struct *fd1, *fd2, *fd3 + +begin + string imtype, mstype + string arcref1, arcref2, spec, arc, align="" + string arcref1ms, arcref2ms, specms, arcms, response + string objs, temp, temp1, done + string str1, str2, str3, str4, arcrefs, log1, log2 + bool scat, reextract, extract, disp, disperr, sky, log + bool skyedit1, skyedit2, 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, apidtable, apref, flat, throughput, arcs1, arcs2, + scattered, dispcor, skysubtract, redo, update) + bye + } + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + mstype = ".ms" // imtype + n = strlen (imtype) + + # Temporary files used repeatedly in this script. Under some + # abort circumstances these files may be left behind. + + objs = mktemp ("tmp$iraf") + temp = mktemp ("tmp$iraf") + temp1 = mktemp ("tmp$iraf") + done = mktemp ("tmp$iraf") + + if (apidtable != "") { + i = strlen (apidtable) + if (i > n && substr (apidtable, i-n+1, i) == imtype) { + apidtable = substr (apidtable, 1, i-n) + i = strlen (apidtable) + } + for (j=1; j<=i && substr(apidtable,j,j)==" "; j+=1); + apidtable = substr (apidtable, j, i) + } + i = strlen (apidtable) + if (i == 0) + extn = ".ms" + else { + extn = apidtable + while (yes) { + i = stridx ("/$]", extn) + if (i == 0) + break + j = strlen (extn) + extn = substr (extn, i+1, j) + } + extn = extn // ".ms" + } + + # Get query parameter. + getspec (objects, objs) + if (arctable == "" || arctable == " ") { + if (arcs2 == "" || arcs2 == " ") + arcrefs = arcs1 + else + arcrefs = arcs2 + } else + arcrefs = arctable + arcref1 = "" + arcref2 = "" + + # 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) + + getspec (apref, temp) + fd1 = temp + if (fscan (fd1, apref) == EOF) + error (1, "No aperture reference") + fd1 = ""; delete (temp, verify=no) + + # Initialize + apscript.saturation = INDEF + apscript.references = apref + apscript.profiles = "" + apscript.apidtable = apidtable + apscript.nfind = fibers + apscript.clean = clean + if (splot) { + splot1 = yes + splot2 = yes + } else { + splot1 = no + splot2 = no + } + + reextract = redo + if (reextract || !access (database // "/ap" // apref // extn)) { + if (!access (apref // imtype)) { + printf ("Aperture reference spectrum not found - %s%s\n", + apref, imtype) | scan (err) + error (1, err // "\nCheck settting of imtype") + } + print ("Set reference apertures for ", apref) | tee (log1) + if (access (database // "/ap" // apref)) + delete (database // "/ap" // apref, verify=no) + if (access (database // "/ap" // apref//extn)) + delete (database // "/ap" // apref // extn, verify=no) + apscript.ansresize = "yes" + apscript.ansedit = "yes" + apscript.ansfittrace = "yes" + apscript (apref, references="", ansfind="YES", ansrecenter="NO", + anstrace="YES", ansextract="NO") + newaps = yes + copy (database//"/ap"//apref, database//"/ap"//apref//extn, + verbose=no) + } else { + if (access (database // "/ap" // apref)) + delete (database // "/ap" // apref, verify=no) + copy (database//"/ap"//apref//extn, database//"/ap"//apref, + verbose=no) + } + + 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" + if (scattered) { + apscript.ansfitscatter = "yes" + apscript.ansfitsmooth = "yes" + } + apscript.ansfittrace = "NO" + apscript.ansextract = "YES" + apscript.ansreview = "NO" + if (skyedit) { + skyedit1 = yes + skyedit2 = yes + } else { + skyedit1 = no + skyedit2 = 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 and 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 + + getspec (spec, temp) + fd1 = temp + if (fscan (fd1, spec) == EOF) + error (1, "No flat field") + fd1 = ""; delete (temp, verify=no) + + scat = no + if (scattered) { + if (redo && access (spec//"noscat"//imtype)) { + imdelete (spec, verify=no) + imrename (spec//"noscat", spec) + } + hselect (spec, "apscatte", yes) | scan (str1) + if (nscan() == 0) + scat = yes + } + if (scat) { + print ("Subtract scattered light from ", spec) | tee (log1) + #apscript.ansfitscatter = "yes" + #apscript.ansfitsmooth = "yes" + imrename (spec, spec//"noscat") + apscript (spec//"noscat", output=spec, ansextract="NO", + ansscat="YES", anssmooth="YES") + #apscript.ansfitscatter = "NO" + #apscript.ansfitsmooth = "NO" + } + + # The next step is to process the flat field image which is used + # as a flat field and a throughput correction. + + spec = throughput + i = strlen (spec) + if (i > n && substr (spec, i-n+1, i) == imtype) + spec = substr (spec, 1, i-n) + specms = spec // mstype + + if (spec != "" && access (spec//imtype)) { + getspec (spec, temp) + fd1 = temp + if (fscan (fd1, spec) == EOF) + error (1, "No flat field") + fd1 = ""; delete (temp, verify=no) + + scat = no + if (scattered) { + if (redo && access (spec//"noscat"//imtype)) { + imdelete (spec, verify=no) + imrename (spec//"noscat", spec) + } + hselect (spec, "apscatte", yes) | scan (str1) + if (nscan() == 0) + scat = yes + } + if (scat) { + print ("Subtract scattered light from ", spec) | tee (log1) + imrename (spec, spec//"noscat") + apscript (spec//"noscat", output=spec, ansextract="NO", + ansscat="YES", anssmooth="YES") + } + } + + response = "" + if (flat != "" || spec != "") { + if (extn == ".ms") + response = flat // spec // "norm.ms" + else + response = flat // spec // extn + reextract = redo || (update && newaps) + if (reextract || !access (response // imtype) || (redo && scat)) { + print ("Create response function ", response) | tee (log1) + + if (access (response // imtype)) + imdelete (response, verify=no) + if (access (flat //mstype)) + imdelete (flat//mstype, verify=no) + if (access (specms)) + imdelete (specms, verify=no) + + fibresponse (flat, spec, 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) { + getspec (arcs1, temp) + fd1 = temp + if (fscan (fd1, arcref1) == EOF) + error (1, "No reference arcs") + fd1 = ""; delete (temp, verify=no) + if (!access (arcref1 // imtype)) { + printf ("Arc reference spectrum not found - %s%s\n", + arcref1, imtype) | scan (err) + error (1, err // "\nCheck settting of imtype") + } + arcref1ms = arcref1 // extn + reextract = redo || (update && newaps) + if (reextract && access (arcref1ms//imtype)) + imdelete (arcref1ms, verify=no) + + getspec (arcs2, temp) + fd1 = temp + if (fscan (fd1, arcref2) == EOF) + arcref2 = "" + else { + if (!access (arcref2 // imtype)) { + printf ("Arc reference spectrum not found - %s%s\n", + arcref2, imtype) | scan (err) + error (1, err // "\nCheck settting of imtype") + } + arcref2ms = arcref2 // extn + if (reextract && access (arcref2ms//imtype)) + imdelete (arcref2ms, verify=no) + } + fd1 = ""; delete (temp, verify=no) + + arcrefs (arcref1, arcref2, extn, arcreplace, apidtable, response, + crval, cdelt, done, log1, log2) + + # Define alignment if needed. + if (skyalign) { + align = "align" // extn // imtype + if (reextract) + imdelete (align, verify=no, >& "dev$null") + } + } + + # Now we are ready to process the object spectra. + + reextract = redo || (update && (newaps || newresp || newdisp)) + fd1 = objs + while (fscan (fd1, spec) != EOF) { + # Check if previously done; i.e. arc. + if (access (done)) { + fd2 = done + while (fscan (fd2, specms) != EOF) + if (spec == specms) + break + if (spec == specms) + next + fd2 = "" + } + if (!access (spec // imtype)) { + printf ("Object spectrum not found - %s%s\n", + spec, imtype) | tee (log1) + print ("Check settting of imtype") + } + specms = spec // mstype + + # Determine required operations from the flags and image header. + scat = no + extract = no + disp = no + sky = no + if (scattered) { + if (redo && access (spec//"noscat"//imtype)) { + imdelete (spec, verify=no) + imrename (spec//"noscat", spec) + } + hselect (spec, "apscatte", yes) | scan (str1) + if (nscan() == 0) + scat = yes + } + if (reextract || !access (specms) || (redo && scat)) + extract = yes + else { + hselect (specms, "dclog1", yes) | scan (str1) + if (nscan () == 1) { + extract = update && newdisp + if (update && !newdisp) + # We really should check if REFSPEC will assign + # different reference spectra. + ; + } else + disp = dispcor + + hselect (specms, "skysub", yes) | scan (str1) + if (nscan() == 0) + sky = skysubtract + } + + if (extract) { + disp = dispcor + sky = skysubtract + } + + # If fully processed go to the next object. + if (!extract && !disp && !sky) + next + + # If not interactive and the batch flag is set submit rest to batch. + if (batch && !skyedit1 && !skyedit2 && !splot1 && !splot2 && + apscript.ansedit == "NO" && (!skyalign || access (align))) { + fd1 = ""; delete (objs, verify=no) + apscript.ansfitscatter = "NO" + apscript.ansfitsmooth = "NO" + + flprcache + batch.objects = objects + batch.datamax = datamax + batch.response = response + batch.arcs1 = arcs1 + batch.arcs2 = arcs2 + batch.arcref1 = arcref1 + batch.arcref2 = arcref2 + batch.arcreplace = arcreplace + batch.apidtable = apidtable + batch.arcrefs = arcrefs + batch.extn = extn + batch.objaps = objaps + batch.skyaps = skyaps + batch.arcaps = arcaps + batch.objbeams = objbeams + batch.skybeams = skybeams + 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.savearcs = savearcs + batch.skyalign = skyalign + batch.skysubtract = skysubtract + batch.saveskys = saveskys + batch.newaps = newaps + batch.newresp = newresp + batch.newdisp = newdisp + batch.newarcs = newarcs + dobatch = yes + return + } + + # Process the spectrum in foreground. + if (extract) { + if (access (specms)) + imdelete (specms, verify=no) + + if (scat) { + print ("Subtract scattered light from ", spec) | tee (log1) + imrename (spec, spec//"noscat") + apscript (spec//"noscat", output=spec, ansextract="NO", + ansscat="YES", anssmooth="YES") + } + + print ("Extract object spectrum ", spec) | tee (log1) + hselect (spec, "date-obs,ut,exptime", yes, > temp1) + hselect (spec, "ra,dec,epoch,st", yes, >> temp1) + fd3 = temp1 + if (fscan (fd3, str1, str2, str3) == 3) { + setjd (spec, observatory=observatory, date="date-obs", + time="ut", exposure="exptime", jd="jd", hjd="", + ljd="ljd", utdate=yes, uttime=yes, listonly=no, + >> log1) + if (fscan (fd3, str1, str2, str3, str4) == 4) + setairmass (spec, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, update=yes, + override=yes, >> log1) + } + fd3 = ""; delete (temp1, verify=no) + apscript (spec, nsubaps=params.nsubaps, saturation=datamax) + sapertures (specms, apertures="", apidtable=apidtable, + wcsreset=no, verbose=no, beam=INDEF, dtype=INDEF, w1=INDEF, + dw=INDEF, z=INDEF, aplow=INDEF, aphigh=INDEF, title=INDEF) + if (response != "") { + if (params.nsubaps == 1) + sarith (specms, "/", response, specms, w1=INDEF, + w2=INDEF, apertures="", bands="", beams="", + apmodulus=0, reverse=no, ignoreaps=no, + format="multispec", renumber=no, offset=0, + clobber=yes, merge=no, errval=0, verbose=no) + else { + blkrep (response, temp, 1, params.nsubaps) + sarith (specms, "/", temp, specms, w1=INDEF, + w2=INDEF, apertures="", bands="", beams="", + apmodulus=0, reverse=no, ignoreaps=yes, + format="multispec", renumber=no, offset=0, + clobber=yes, merge=no, errval=0, verbose=no) + imdelete (temp, verify=no) + } + } + } + + disperr = no + if (disp) { + # Fix arc headers if necessary. + if (newarcs) { + getspec (arcs1, temp) + fd2 = temp + while (fscan (fd2, arc) != EOF) { + hselect (arc, "date-obs,ut,exptime", yes, > temp1) + hselect (arc, "ra,dec,epoch,st", yes, >> temp1) + fd3 = temp1 + if (fscan (fd3, str1, str2, str3) == 3) { + setjd (arc, observatory=observatory, date="date-obs", + time="ut", exposure="exptime", jd="jd", hjd="", + ljd="ljd", utdate=yes, uttime=yes, listonly=no, + >> log1) + if (fscan (fd3, str1, str2, str3, str4) == 4) + setairmass (arc, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, update=yes, + override=yes, >> log1) + } + fd3 = ""; delete (temp1, verify=no) + 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 (temp, verify=no) + getspec (arcs2, temp) + fd2 = temp + while (fscan (fd2, arc) != EOF) { + hselect (arc, "date-obs,ut,exptime", yes, > temp1) + hselect (arc, "ra,dec,epoch,st", yes, >> temp1) + fd3 = temp1 + if (fscan (fd3, str1, str2, str3) == 3) { + setjd (arc, observatory=observatory, + date="date-obs", time="ut", exposure="exptime", + jd="jd", hjd="", ljd="ljd", utdate=yes, + uttime=yes, listonly=no, >> log1) + if (fscan (fd3, str1, str2, str3, str4) == 4) + setairmass (arc, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, + update=yes, override=yes, >> log1) + + } + fd3 = ""; delete (temp1, verify=no) + hedit (arc, "refspec1", arc, add=yes, verify=no, + show=no, update=yes) + hedit (arc, "arctype", "shift", add=yes, verify=no, + show=no, update=yes) + } + fd2 = ""; delete (temp, 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, arcref1, arcref2, extn, arcreplace, + apidtable, arcaps, arcbeams, savearcs, reextract, arcap, + log1, no, done) + + hselect (specms, "refspec1", yes, > temp) + fd2 = temp + i = fscan (fd2, arc) + fd2 = ""; delete (temp, verify=no) + if (i < 1) { + print ("No arc reference assigned for ", spec) | tee (log1) + disperr = yes + } else { + if (skyalign) + doalign (spec, specms, align, arcref1ms, log1, no) + print ("Dispersion correct ", spec) | tee (log1) + dispcor (specms, "", linearize=params.linearize, + database=database, table=arcref1ms, w1=INDEF, + w2=INDEF, dw=INDEF, nw=INDEF, log=params.log, + flux=params.flux, samedisp=no, global=no, + ignoreaps=no, confirm=no, listonly=no, + verbose=verbose, logfile=logfile) + if (params.nsubaps > 1) { + imrename (specms, temp, verbose=no) + scopy (temp, specms, w1=INDEF, w2=INDEF, + apertures="1-999", bands="", beams="", apmodulus=0, + offset=0, format="multispec", clobber=no, merge=no, + renumber=no, verbose=no) + blkavg (temp, temp, 1, params.nsubaps, option="sum") + imcopy (temp, specms//"[*,*]", verbose=no) + imdelete (temp, verify=no) + } + } + } + + if (sky && !disperr) { + str1 = "" + if (skyaps != "") + str1 = "skyaps=" // skyaps + if (skybeams != "") + str1 = str1 // " skybeams=" // skybeams + print ("Sky subtract ", spec, ": ", str1) | tee (log1) + if (skyedit1) { + str1 = substr (ansskyedit, 1, 1) + if (str1 == "N" || str1 == "Y") + skyedit1 = no + if (str1 == "n" || str1 == "N") + skyedit2 = no + else + skyedit2 = yes + } + skysub.reject = params.reject + skysub (specms, output="", objaps=objaps, skyaps=skyaps, + objbeams=objbeams, skybeams=skybeams, skyedit=skyedit2, + combine=params.combine, scale=params.scale, + saveskys=saveskys, logfile=logfile) + params.reject = skysub.reject + hedit (specms, "skysub", yes, add=yes, show=no, verify=no, + update=yes) + } + + if (!disperr && (extract || disp || sky)) { + if (splot1) { + print (specms, ":") + str1 = anssplot + if (str1 == "NO" || str1 == "YES") + splot1 = no + if (str1 == "no" || str1 == "NO") + splot2 = no + else + splot2 = yes + } + if (splot2) + splot (specms) + } + + print (spec, >> done) + } + fd1 = ""; delete (objs, verify=no) + + if (access (done)) + delete (done, verify=no) +end |