diff options
Diffstat (limited to 'noao/imred/src/fibers')
24 files changed, 3177 insertions, 0 deletions
diff --git a/noao/imred/src/fibers/Revisions b/noao/imred/src/fibers/Revisions new file mode 100644 index 00000000..a21892a2 --- /dev/null +++ b/noao/imred/src/fibers/Revisions @@ -0,0 +1,223 @@ +.help revisions Nov90 noao.imred.src.fibers +.nf + +srcfibers$proc.cl +srcfibers$listonly.cl + For reasons that are now lost in the mists of time, the shortening of + filenames can cause problems. This has been removed. + (2/11/11, Valdes) (See bug log 577) + +======= +V2.15.1 +======= + +======= +V2.12.3 +======= + +srcfibers$batch.cl +srcfibers$proc.cl +srcfibers$fibresponse.cl + Error messages now hint to check imtype setting. + (4/15/05, Valdes) + + +======= +V2.12.2 +======= + +srcfibers$proc.cl +srcfibers$batch.cl + Because dispcor.samedisp=yes when doing the objects any dispersion change + applied to the reference arc was not being applied to the objects. + (5/21/03, Valdes) + +======= +V2.12.1 +======= + +srcfibers$skysub.cl + Added a flpr to workaround a FITS kernel caching problem. + (6/21/02, Valdes) + +===== +V2.12 +===== + +srcfibers$proc.cl + Modified code to eliminate goto. This is for use with pyraf. + (11/21/00, Valdes) + +======== +V2.11.3a +======== + +srcfibers$arcrefs.cl + The test for crval/cdelt both INDEF was not working. (9/24/99, Valdes) + +srcfibers$mkfibers.cl + Had to define a dummy variable 'ar' to get rid of ambiguous parameter + error. (9/14/99, Valdes) + +======= +V2.11.2 +======= + +doslit$arcrefs.cl + The test on CRVAL and CDELT would not work with header keywords. + (9/22/98, Valdes) + +srcfibers$arcrefs.cl +srcfibers$batch.cl +srcfibers$doarcs.cl +srcfibers$fibresponse.cl +srcfibers$getspec.cl +srcfibers$listonly.cl +srcfibers$mkfibers.cl +srcfibers$skysub.cl +srcfibers$proc.cl + Any additional qualifiers in the imtype string are stripped. + (8/14/97, Valdes) + +========= +V2.11BETA +========= + +arcrefs.cl + If both crval and cdelt are INDEF then the autoidentify option is + not used. (5/2/97, Valdes) + +apsript.par + Made changes for the new aperture selection option. (9/5/96, Valdes) + +skysub.cl + Added package name to calls of "match", "sort", "uniq" to avoid + possible conflicts. (5/6/96, Valdes) + +proc.cl +proc.par +arcrefs.cl +arcrefs.par +params.par + Modified to use autoidentify. (4/5/96, Valdes) + +srcfibers$proc.cl +srcfibers$batch.cl + When using subapertures the subapertures were not wavelength + calibrated correctly because the reference arc spectrum which + provides the wavelength scale does not contain the subapertures + and DISPCOR does not use samedisp=yes. Changed the value of + samedisp to yes. (10/27/95, Valdes) + +srcfibers$mkfibers.cl + The calls to mk1dspec did not specify the header file which would + then default to the task parameter which might be invalid. + (10/17/95, Valdes) + +srcfibers$proc.cl + Needed to initialize arcref2 in order work in batch when no dispersion + correction is requested. (10/16/95, Valdes) + +srcfibers$mkfibers.cl + The calls to MK1DSPEC were changed in accordance with parameter changes + to that task. + (7/28/95, Valdes) + +srcfibers$proc.cl + Any image extension is stripped from the apidtable parameter. + (7/24/95, Valdes) + +srcfibers$doalign.cl + +srcfibers$doalign.par + +srcfibers$proc.cl +srcfibers$batch.cl + Added the sky alignment option. (7/19/95, Valdes) + +srcfibers$proc.cl +srcfibers$batch.cl +srcfibers$arcrefs.cl + The wrong range syntax is used with subapertures in SARITH/SCOPY. + Changed all -999 to 1-999. (6/14/95, Valdes) + +======= +V2.10.4 +======= + +srcfibers$proc.cl +srcfibers$fibresponse.cl + 1. Need to add check for the throughput being a file rather than + an image when checking whether to apply a scattered light + correction. + 2. Removed a warning message when using a throughput file containing + fiber values which are not in the flat field (for example, if a fiber + is broken). + (1/25/95, Valdes) + +srcfibers$params.par +srcfibers$doarcs.cl +srcfibers$arcrefs.cl + Added "threshold" as a user parameter. (1/16/95, Valdes) + +srcfibers$response.cl -> imred$src/fibers/fibresponse.cl +srcfibers$response.par -> imred$src/fibers/fibresponse.par +srcfibers$proc.par + Changed the fiber response task name from "response" to "fibresponse" + to avoid conflict with longslit.response. (12/31/94, Valdes) + +srcfibers$proc.cl + The check for arcs2 = " " was done incorrectly. (9/12/94, Valdes) + +srcfibers$proc.cl +srcfibers$batch.cl +srcfibers$doarcs.cl + A check was needed on whether the arc spectra were extracted during + the current execution to avoid reextracting the same arc multiple + times during a "redo" or the initial time. In both those cases + the rextract flag is set causing spectra to be reextracted if they + exist. Previously doarcs could not tell if the arc exists because + of a previous run or during the current run with the same arc + used multiple times. (5/18/94, Valdes) + +=========== +V2.10.3beta +=========== + +srcfibers$skysub.cl +imred$specred/doc/skysub.hlp + 1. The combine option was being ignored. + 2. The help did not mention the reject option and was otherwised + out of date. + (3/31/94, Valdes) + +srcfibers$proc.cl + The scattered light correction is now queried for all images and may + be turned off with NO. (9/1/93, Valdes) + +=========== +V2.10.3beta +=========== + +srcfibers$arcrefs.cl + MOdified to use shift=INDEF in REIDENTIFY. + (2/18/93, Valdes) + +srcfibers$*.cl + Modified to use the "imtype" environment variable to define the + extension type. + (2/18/93, Valdes) + +======= +V2.10.2 +======= + +srcfibers$proc.cl + The aperture reference is redone when a new aperture ID file is seen. + (1/11/93, Valdes) + +srcfibers$* + Updated for new ONEDSPEC. (7/24/91, Valdes) + +srcfibers$* + All occurrences of latitude replaced by observatory as required by + recent changes to setairmass, etc. (11/20/90, Valdes) +.endhelp diff --git a/noao/imred/src/fibers/apscript.par b/noao/imred/src/fibers/apscript.par new file mode 100644 index 00000000..e52248de --- /dev/null +++ b/noao/imred/src/fibers/apscript.par @@ -0,0 +1,145 @@ +# APSCRIPT + +input,s,a,,,,List of input images +output,s,h,"",,,List of output spectra +apertures,s,h,"",,,Apertures +scatter,s,h,"",,,List of scattered light images (optional) +references,s,h,"",,,List of aperture reference images +profiles,s,h,"",,,"List of aperture profile images +" +interactive,b,h,yes,,,Run task interactively? +find,b,h,yes,,,Find apertures? +recenter,b,h,yes,,,Recenter apertures? +resize,b,h,yes,,,Resize apertures? +edit,b,h,yes,,,Edit apertures? +trace,b,h,yes,,,Trace apertures? +fittrace,b,h,yes,,,Fit the traced points interactively? +extract,b,h,yes,,,Extract spectra? +review,b,h,yes,,,Review extractions? +subtract,b,h,yes,,,Subtract scattered light? +smooth,b,h,yes,,,Smooth scattered light along the dispersion? +fitscatter,b,h,yes,,,Fit scattered light interactively? +fitsmooth,b,h,yes,,,"Smooth the scattered light interactively? +" +line,i,h,)params.line,,,>params.line +nsum,i,h,)params.nsum,,,>params.nsum +buffer,r,h,)params.buffer,,,">params.buffer + +# OUTPUT PARAMETERS +" +format,s,h,"multispec",,,Extracted spectra format +extras,b,h,)params.extras,,,"Extract sky, sigma, etc.?" +dbwrite,s,h,"YES",,,Write to database? +initialize,b,h,no,,,Initialize answers? +verbose,b,h,no,,,"Verbose output? + +# DEFAULT APERTURE PARAMETERS +" +lower,r,h,)params.lower,,,>params.lower +upper,r,h,)params.upper,,,>params.upper +apidtable,s,h,"",,,"Aperture ID table (optional) + +# DEFAULT BACKGROUND PARAMETERS +" +b_function,s,h,"legendre","chebyshev|legendre|spline1|spline3",,Background function +b_order,i,h,1,,,Background function order +b_sample,s,h,"-10:-6,6:10",,,Background sample regions +b_naverage,i,h,-3,,,Background average or median +b_niterate,i,h,0,0,,Background rejection iterations +b_low_reject,r,h,3.,0.,,Background lower rejection sigma +b_high_reject,r,h,3.,0.,,Background upper rejection sigma +b_grow,r,h,0.,0.,,"Background rejection growing radius + +# APERTURE CENTERING PARAMETERS +" +width,r,h,5.,0.,,Profile centering width +radius,r,h,10.,,,Profile centering radius +threshold,r,h,10.,0.,,"Detection threshold for profile centering + +# AUTOMATIC FINDING AND ORDERING PARAMETERS +" +nfind,i,h,,,,Number of apertures to be found automatically +minsep,r,h,1.,,,Minimum separation between spectra +maxsep,r,h,100000.,,,Maximum separation between spectra +order,s,h,)params.order,,,"Order of apertures + +# RECENTERING PARAMETERS +" +aprecenter,s,h,"",,,Apertures to use in recentering +npeaks,r,h,0.5,,,Select brightest peaks +shift,b,h,yes,,,"Use average shift instead of recentering? + +# RESIZING PARAMETERS +" +llimit,r,h,INDEF,,,Lower aperture limit relative to center +ulimit,r,h,INDEF,,,Upper aperture limit relative to center +ylevel,r,h,)params.ylevel,,,>params.ylevel +peak,b,h,yes,,,Is ylevel a fraction of the peak? +bkg,b,h,yes,,,"Subtract background in automatic width?" +r_grow,r,h,0.,,,Grow limits by this factor +avglimits,b,h,no,,,"Average limits over all apertures? + +# EDITING PARAMETERS +" +e_output,s,q,,,,Output spectra rootname +e_profiles,s,q,,,,"Profile reference image + +# TRACING PARAMETERS +" +t_nsum,i,h,)params.nsum,,,>params.nsum +t_step,i,h,)params.t_step,,,>params.t_step +t_nlost,i,h,3,1,,Number of consecutive times profile is lost before quitting +t_width,r,h,5.,0.,,Profile centering width +t_function,s,h,)params.t_function,,,>params.t_function +t_sample,s,h,"*",,,Trace sample regions +t_order,i,h,)params.t_order,,,>params.t_order +t_naverage,i,h,1,,,Trace average or median +t_niterate,i,h,)params.t_niterate,,,>params.t_niterate +t_low_reject,r,h,)params.t_low,,,>params.t_low +t_high_reject,r,h,)params.t_high,,,>params.t_high +t_grow,r,h,0.,0.,,"Trace rejection growing radius + +# EXTRACTION PARAMETERS +" +background,s,h,"none","none|average|median|minimum|fit",,Background to subtract +skybox,i,h,1,,,Box car smoothing length for sky +weights,s,h,)params.weights,,,>params.weights +pfit,s,h,)params.pfit,,,>params.pfit +clean,b,h,no,,,Detect and replace bad pixels? +nclean,r,h,0.5,,,Maximum number of pixels to clean +niterate,i,h,5,0,,Number of profile fitting iterations +saturation,r,h,INDEF,,,Saturation level +readnoise,s,h,"0.",,,Read out noise sigma (photons) +gain,s,h,"1.",,,Photon gain (photons/data number) +lsigma,r,h,)params.lsigma,,,>params.lsigma +usigma,r,h,)params.usigma,,,>params.usigma +polysep,r,h,0.95,0.1,0.95,Marsh algorithm polynomial spacing +polyorder,i,h,10,1,,Marsh algorithm polynomial order +nsubaps,i,h,1,1,,"Number of subapertures per aperture + +# ANSWER PARAMETERS +" +ansclobber,s,h,"NO",,," " +ansclobber1,s,h,"NO",,," " +ansdbwrite,s,h,"YES",,," " +ansdbwrite1,s,h,"NO",,," " +ansedit,s,h,"NO",,," " +ansextract,s,h,"NO",,," " +ansfind,s,h,"NO",,," " +ansfit,s,h,"NO",,," " +ansfitscatter,s,h,"NO",,," " +ansfitsmooth,s,h,"NO",,," " +ansfitspec,s,h,"NO",,," " +ansfitspec1,s,h,"NO",,," " +ansfittrace,s,h,"NO",,," " +ansfittrace1,s,h,"NO",,," " +ansflat,s,h,"NO",,," " +ansnorm,s,h,"NO",,," " +ansrecenter,s,h,"NO",,," " +ansresize,s,h,"NO",,," " +ansreview,s,h,"NO",,," " +ansreview1,s,h,"NO",,," " +ansscat,s,h,"NO",,," " +ansskyextract,s,h,"NO",,," " +anssmooth,s,h,"NO",,," " +anstrace,s,h,"NO",,," " diff --git a/noao/imred/src/fibers/arcrefs.cl b/noao/imred/src/fibers/arcrefs.cl new file mode 100644 index 00000000..d0d6fcf1 --- /dev/null +++ b/noao/imred/src/fibers/arcrefs.cl @@ -0,0 +1,326 @@ +# ARCREFS -- Determine dispersion relation for reference arcs. + +procedure arcrefs (arcref1, arcref2, extn, arcreplace, apidtable, response, + crval, cdelt, done, log1, log2) + +file arcref1 +file arcref2 +string extn +file arcreplace +file apidtable +file response +string crval = "INDEF" +string cdelt = "INDEF" +file done +file log1 +file log2 + +struct *fd + +begin + string imtype + string arcref, arcrefms, arc, arcms + string temp, temp1, temp2, str1, str2 + int i, n, nspec, dc + real w + bool log + struct str3 + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + n = strlen (imtype) + + temp = mktemp ("tmp$iraf") + temp1 = mktemp ("tmp$iraf") + temp2 = mktemp ("tmp$iraf") + + # Extract the primary arc reference spectrum. Extract and replace + # any replacement arcs defined in the arcreplace file. Determine the + # dispersion function with IDENTIFY/REIDENTIFY. Set the wavelength + # parameters with MSDISPCOR. + + arcref = arcref1 + arcrefms = arcref1 // extn + if (!access (arcrefms//imtype)) { + print ("Extract arc reference image ", arcref) | tee (log1) + apscript (arcref, output=arcrefms, ansrecenter="NO", + ansresize="NO", ansedit="NO", anstrace="NO", + nsubaps=params.nsubaps, background="none", clean=no, + weights="none") + sapertures (arcrefms, 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 (arcrefms, "/", response, arcrefms, 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 (arcrefms, "/", temp, arcrefms, 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) + } + } + + if (arcreplace != "") { + if (!access (arcreplace)) + error (1, "Can't access file "//arcreplace) + fd = arcreplace + while (fscan (fd, arc, str1, str2) != EOF) { + i = strlen (arc) + if (i > n && substr (arc, i-n+1, i) == imtype) + arc = substr (arc, 1, i-n) + if (arc != arcref) + next + arc = str1 + if (i > n && substr (arc, i-n+1, i) == imtype) + arc = substr (arc, 1, i-n) + arcms = arc // extn + + if (access (arcms//imtype)) + imdelete (arcms, verify=no) + + print ("Extract arc reference image ", arc) | tee (log1) + apscript (arc, output=arcms, ansrecenter="NO", + ansresize="NO", ansedit="NO", anstrace="NO", + nsubaps=params.nsubaps, background="none", clean=no, + weights="none") + sapertures (arcms, 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 (arcms, "/", response, arcms, 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 (arcms, "/", temp, arcms, 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) + } + } + scopy (arcms, arcrefms, w1=INDEF, w2=INDEF, apertures=str2, + bands="", beams="", apmodulus=1000, offset=0, + format="multispec", clobber=yes, merge=yes, renumber=no, + verbose=yes) | tee (log1, > log2) + imdelete (arcms, verify=no) + } + fd = "" + } + } + + # Get the dispersion parameters from the header. These are + # used for all further spectra and also flag whether this + # spectrum has been processed. If the parameters are missing + # the spectrum needs to have the dispersion function and + # wavelength scale determined. The HEDIT is needed because + # in some cases the user may exit IDENTIFY without updating + # the database (if the image was deleted but the database + # entry was not). + + hselect (arcrefms, "dclog1", yes) | scan (str1) + if (nscan () != 1) { + print ("Determine dispersion solution for ", arcref) | tee (log1) + #delete (database//"/id"//arcrefms//"*", verify=no) + printf ("%s %s\n", crval, cdelt) | scan (str3) + if (str3 == "INDEF INDEF") + identify (arcrefms, section="middle line", database=database, + coordlist=params.coordlist, nsum=1, match=params.match, + maxfeatures=50, zwidth=100., ftype="emission", + fwidth=params.fwidth, cradius=params.cradius, + threshold=params.threshold, minsep=2., + function=params.i_function, order=params.i_order, + sample="*", niterate=params.i_niterate, + low_reject=params.i_low, high_reject=params.i_high, + grow=0., autowrite=yes) + else + autoidentify (arcrefms, crval, cdelt, + coordlist=params.coordlist, + interactive="YES", section="middle line", nsum="1", + ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., match=params.match, function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, low_reject=params.i_low, + high_reject=params.i_high, grow=0., dbwrite="YES", + overwrite=yes, database="database", verbose=yes, + logfile=logfile, plotfile=plotfile, + reflist="", refspec="", crpix="INDEF", cddir="unknown", + crsearch="-0.5", cdsearch="INDEF", aidpars="") + + hedit (arcrefms, "refspec1", arcrefms, add=yes, + show=no, verify=no) + + nspec = 1 + hselect (arcrefms, "naxis2", yes) | scan (nspec) + if (nspec > 1) + reidentify (arcrefms, "", interactive=yes, + section="middle line", shift=0., step=1, nsum=1, + cradius=params.cradius, threshold=params.threshold, + nlost=100, refit=params.refit, trace=no, override=yes, + addfeatures=params.addfeatures, newaps=no, + coordlist=params.coordlist, match=params.match, + maxfeatures=50, minsep=2., database=database, + plotfile=plotfile, logfiles=logfile, verbose=yes) + + # Dispersion correct the reference arc. This step is required to + # use the confirm option of MSDISPCOR to set the wavelength scale + # for all further spectra. Set the newdisp flag. + + print ("Dispersion correct ", arcref) | tee (log1) + dispcor (arcrefms, "", linearize=params.linearize, + database=database, table="", w1=INDEF, w2=INDEF, dw=INDEF, + nw=INDEF, log=params.log, flux=params.flux, samedisp=yes, + global=no, ignoreaps=no, confirm=yes, listonly=no, verbose=no, + logfile=logfile) + if (params.nsubaps > 1) { + imrename (arcrefms, temp, verbose=no) + scopy (temp, arcrefms, 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, arcrefms//"[*,*]", verbose=no) + imdelete (temp, verify=no) + } + proc.newdisp = yes + } + if (extn == ".ms") + print (arcref, >> done) + + # Extract the alternate shift arc reference. Transfer the dispersion + # function from the primary arc reference and then identify shift + # lines. + + if (arcref2 != "") { + arcref = arcref2 + arcrefms = arcref2 // extn + if (proc.newdisp && access (arcrefms//imtype)) + imdelete (arcrefms, verify=no) + if (!access (arcrefms)) { + print ("Extract arc reference image ", arcref) | tee (log1) + apscript (arcref, output=arcrefms, ansrecenter="NO", + ansresize="NO", ansedit="NO", anstrace="NO", + nsubaps=params.nsubaps, background="none", clean=no, + weights="none") + sapertures (arcrefms, 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 (arcrefms, "/", response, arcrefms, 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 (arcrefms, "/", temp, arcrefms, 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) + } + } + } + + hselect (arcrefms, "dclog1", yes) | scan (str1) + if (nscan () != 1) { + print ("Determine dispersion solution for ", arcref) | + tee (log1) + #delete (database//"/id"//arcrefms//"*", verify=no) + + print (":r ", arcref1//extn, "\na\nd") | + identify (arcrefms, section="middle line", database=database, + coordlist="", nsum=1, match=params.match, maxfeatures=50, + zwidth=100., ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, low_reject=params.i_low, + high_reject=params.i_high, grow=0., autowrite=yes, + cursor="STDIN", >G "dev$null", >& "dev$null") + identify (arcrefms, section="middle line", database=database, + coordlist="", nsum=1, match=params.match, maxfeatures=50, + zwidth=100., ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, low_reject=params.i_low, + high_reject=params.i_high, grow=0., autowrite=yes) + print (":feat ", temp) | + identify (arcrefms, section="middle line", database=database, + coordlist="", nsum=1, match=params.match, maxfeatures=50, + zwidth=100., ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, low_reject=params.i_low, + high_reject=params.i_high, grow=0., autowrite=yes, + cursor="STDIN", >G "dev$null", >& "dev$null") + print (":r ", arcref1//extn, "\na\nd", > temp1) + fd = temp + while (fscan (fd, i, w, w, w) != EOF) { + if (nscan() == 4) { + print (w, 1, 1, "m", >> temp1) + print (w, >> temp2) + } + } + print ("g", >> temp1) + fd = ""; delete (temp, verify=no) + + nspec = 1 + hselect (arcrefms, "naxis2", yes) | scan (nspec) + for (i = 1; i <= nspec; i+=1) + identify (arcrefms, section="line "//i, + database=database, coordlist="", nsum=1, + match=params.match, maxfeatures=50, zwidth=100., + ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, + low_reject=params.i_low, high_reject=params.i_high, + grow=0., autowrite=yes, cursor=temp1, < temp2, + >G "dev$null", >>& temp) + delete (temp1, verify=no); delete (temp2, verify=no) + system.match ("Coordinate shift", temp, stop=no, print_file_n=yes, + metacharacte=yes) | tee (log1, > log2) + delete (temp, verify=no) + + dispcor (arcrefms, "", linearize=params.linearize, + database=database, table="", w1=INDEF, w2=INDEF, + dw=INDEF, nw=INDEF, log=params.log, flux=params.flux, + samedisp=yes, global=no, ignoreaps=no, confirm=no, + listonly=no, verbose=yes, logfile=logfile, > log2) + if (params.nsubaps > 1) { + imrename (arcrefms, temp, verbose=no) + scopy (temp, arcrefms, 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, arcrefms//"[*,*]", verbose=no) + imdelete (temp, verify=no) + } + } + if (extn == ".ms") + print (arcref, >> done) + } +end diff --git a/noao/imred/src/fibers/arcrefs.par b/noao/imred/src/fibers/arcrefs.par new file mode 100644 index 00000000..56f69145 --- /dev/null +++ b/noao/imred/src/fibers/arcrefs.par @@ -0,0 +1,13 @@ +arcref1,f,a,"",,, +arcref2,f,a,"",,, +extn,s,a,,,, +arcreplace,f,a,"",,, +apidtable,f,a,"",,, +response,f,a,"",,, +crval,s,a,INDEF,,, +cdelt,s,a,INDEF,,, +done,f,a,"",,, +log1,f,a,"",,, +log2,f,a,"",,, +fd,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/batch.cl b/noao/imred/src/fibers/batch.cl new file mode 100644 index 00000000..84b0a7b7 --- /dev/null +++ b/noao/imred/src/fibers/batch.cl @@ -0,0 +1,297 @@ +# BATCH -- Process spectra in batch. +# This task is called in batch mode. It only processes objects +# not previously processed unless the update or redo flags are set. + +procedure batch () + +string objects {prompt="Object spectra"} +real datamax {prompt="Max data value / cosmic ray threshold"} + +file response {prompt="Response spectrum"} +string arcs1 {prompt="List of arc spectra"} +string arcs2 {prompt="List of shift arc spectra"} +file arcref1 {prompt="Arc reference for dispersion solution"} +file arcref2 {prompt="Arc reference for dispersion solution"} +file arcreplace {prompt="Special aperture replacements"} +string arcrefs {prompt="Arc references"} +string extn {prompt="Extraction extension"} + +file apidtable {prompt="Aperture identifications"} +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"} + +file done {prompt="File of spectra already done"} +file logfile {prompt="Logfile"} + +bool redo {prompt="Redo operations?"} +bool update {prompt="Update spectra?"} +bool scattered {prompt="Subtracted scattered light?"} +bool arcap {prompt="Use object apertures for arcs?"} +bool dispcor {prompt="Dispersion correct spectra?"} +bool savearcs {prompt="Save internal arcs?"} +bool skyalign {prompt="Align sky lines?"} +bool skysubtract {prompt="Subtract sky?"} +bool saveskys {prompt="Save sky spectra?\n"} + +bool newaps, newresp, newdisp, newarcs + +struct *fd1, *fd2, *fd3 + +begin + file objs, temp, temp1, spec, specms, arc + bool reextract, extract, scat, disp, sky, log + string imtype, mstype, str, str2, str3, str4 + int i + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + mstype = ".ms" // imtype + + objs = mktemp ("tmp$iraf") + temp = mktemp ("tmp$iraf") + temp1 = mktemp ("tmp$iraf") + + # Initialize extraction to be noninteractive. + if (apscript.ansrecenter == "yes") + apscript.ansrecenter = "YES" + else if (apscript.ansrecenter == "no") + apscript.ansrecenter = "NO" + apscript.ansedit = "NO" + if (apscript.anstrace == "yes") { + apscript.anstrace = "YES" + apscript.ansfittrace = "NO" + } else if (apscript.anstrace == "no") + apscript.anstrace = "NO" + + reextract = redo || (update && (newaps || newresp || newdisp)) + + getspec (objects, objs) + fd1 = objs + while (fscan (fd1, spec) != EOF) { + 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\nCheck setting of imtype\n", spec, imtype, >> logfile) + next + } + specms = spec // mstype + + 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 (str) + if (nscan() == 0) + scat = yes + } + if (reextract || !access (specms) || (update && scat)) + extract = yes + else { + hselect (specms, "dclog1", yes) | scan (str) + 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 (str) + if (nscan() == 0) + sky = skysubtract + } + + if (extract) { + disp = dispcor + sky = skysubtract + } + + if (extract) { + if (access (specms)) + imdelete (specms, verify=no) + if (scat) { + print ("Subtract scattered light from ", spec, >> logfile) + imrename (spec, spec//"noscat") + apscript (spec//"noscat", output=spec, ansextract="NO", + ansscat="YES", anssmooth="YES", verbose=no) + } + print ("Extract object spectrum ", spec, >> logfile) + hselect (spec, "date-obs,ut,exptime", yes, > temp1) + hselect (spec, "ra,dec,epoch,st", yes, >> temp1) + fd3 = temp1 + if (fscan (fd3, str, 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, + >> logfile) + if (fscan (fd3, str, str2, str3, str4) == 4) + setairmass (spec, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, update=yes, + override=yes, >> logfile) + } + fd3 = ""; delete (temp1, verify=no) + apscript (spec, nsubaps=params.nsubaps, saturation=datamax, + verbose=no) + 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) + } + } + } + + 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, str, 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, + >> logfile) + if (fscan (fd3, str, str2, str3, str4) == 4) + setairmass (arc, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, update=yes, + override=yes, >> logfile) + } + 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, str, 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, >> logfile) + if (fscan (fd3, str, str2, str3, str4) == 4) + setairmass (arc, intype="beginning", + outtype="effective", exposure="exptime", + observatory=observatory, show=no, + update=yes, override=yes, >> logfile) + + } + 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, >> logfile) + 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, + >> logfile) + + doarcs (spec, response, arcref1, arcref2, extn, arcreplace, + apidtable, arcaps, arcbeams, savearcs, reextract, arcap, + logfile, yes, 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, >> logfile) + else { + if (skyalign) + doalign (spec, specms, "align"//extn//imtype, + arcref1//extn, logfile, yes) + print ("Dispersion correct ", spec, >> logfile) + dispcor (specms, "", linearize=params.linearize, + database=database, table=arcref1//extn, 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=no, + 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) + } + disp = no + } + } + + if (sky && !disp) { + str = "" + if (skyaps != "") + str = "skyaps=" // skyaps + if (skybeams != "") + str = str // " skybeams=" // skybeams + print ("Sky subtract ", spec, ": ", str, >> logfile) + skysub (specms, output="", objaps=objaps, skyaps=skyaps, + objbeams=objbeams, skybeams=skybeams, skyedit=no, + combine=params.combine, reject=params.reject, + scale=params.scale, saveskys=saveskys, logfile=logfile) + hedit (specms, "skysub", str, add=yes, show=no, verify=no, + update=yes) + } + } + fd1 = ""; delete (objs, verify=no) + + if (access (done)) + delete (done, verify=no) + + flprcache (0) +end diff --git a/noao/imred/src/fibers/batch.par b/noao/imred/src/fibers/batch.par new file mode 100644 index 00000000..54575594 --- /dev/null +++ b/noao/imred/src/fibers/batch.par @@ -0,0 +1,38 @@ +objects,s,h,,,,"Object spectra" +datamax,r,h,,,,"Max data value / cosmic ray threshold" +response,f,h,"",,,"Response spectrum" +arcs1,s,h,,,,"List of arc spectra" +arcs2,s,h,,,,"List of shift arc spectra" +arcref1,f,h,"",,,"Arc reference for dispersion solution" +arcref2,f,h,"",,,"Arc reference for dispersion solution" +arcreplace,f,h,"",,,"Special aperture replacements" +arcrefs,s,h,,,,"Arc references" +extn,s,h,,,,"Extraction extension" +apidtable,f,h,"",,,"Aperture identifications" +objaps,s,h,,,,"Object apertures" +skyaps,s,h,,,,"Sky apertures" +arcaps,s,h,,,,"Arc apertures" +objbeams,s,h,,,,"Object beam numbers" +skybeams,s,h,,,,"Sky beam numbers" +arcbeams,s,h,,,,"Arc beam numbers +" +done,f,h,"",,,"File of spectra already done" +logfile,f,h,"",,,"Logfile" +redo,b,h,,,,"Redo operations?" +update,b,h,,,,"Update spectra?" +scattered,b,h,,,,"Subtracted scattered light?" +arcap,b,h,,,,"Use object apertures for arcs?" +dispcor,b,h,,,,"Dispersion correct spectra?" +savearcs,b,h,,,,"Save internal arcs?" +skyalign,b,h,,,,"Align sky lines?" +skysubtract,b,h,,,,"Subtract sky?" +saveskys,b,h,,,,"Save sky spectra? +" +newaps,b,h,,,, +newresp,b,h,,,, +newdisp,b,h,,,, +newarcs,b,h,,,, +fd1,*struct,h,"",,, +fd2,*struct,h,"",,, +fd3,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/doalign.cl b/noao/imred/src/fibers/doalign.cl new file mode 100644 index 00000000..e3b8d4db --- /dev/null +++ b/noao/imred/src/fibers/doalign.cl @@ -0,0 +1,78 @@ +# DOALIGN -- Align sky lines in objects. +# If there is no database of features for alignment have user identify +# them interactively. + +procedure doalign (spec, specms, align, table, logfile, batch) + +file spec +file specms +file align +file table +file logfile +bool batch + +begin + file temp + bool log, verbose1 + + if (batch) + verbose1 = no + else + verbose1 = verbose + + if (!access (align)) { + print ("Identify alignment features") + dispcor (specms, align, linearize=no, + database=database, table=table, 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=no, logfile="") + identify (align, section="middle line", database=database, + coordlist="", nsum=1, match=params.match, maxfeatures=50, + zwidth=100., ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, low_reject=params.i_low, + high_reject=params.i_high, grow=0., autowrite=yes) + print ("g") | + identify (align, section="middle line", database=database, + coordlist="", nsum=1, match=params.match, maxfeatures=50, + zwidth=100., ftype="emission", fwidth=params.fwidth, + cradius=params.cradius, threshold=params.threshold, + minsep=2., function=params.i_function, + order=params.i_order, sample="*", + niterate=params.i_niterate, low_reject=params.i_low, + high_reject=params.i_high, grow=0., autowrite=yes, + cursor="STDIN", >G "dev$null", >& "dev$null") + reidentify (align, "", + interactive=no, section="middle line", shift=0., + step=1, nsum=1, cradius=params.cradius, + threshold=params.threshold, nlost=100, newaps=no, + refit=no, trace=no, override=yes, addfeatures=no, + database=database, plotfile=plotfile, + logfiles=logfile, verbose=verbose1) + } + + # Set arc dispersion function in image header. + if (!batch) + print ("Identify alignment features in ", spec) + print ("Identify alignment features in ", spec, >> logfile) + dispcor (specms, "", linearize=no, + database=database, table=table, 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=no, logfile="") + hedit (specms, "refspec1", align, add=yes, + verify=no, show=no, update=yes) + delete (database//"/id"//spec//".ms", verify=no, >& "dev$null") + reidentify (align, specms, + interactive=no, section="middle line", shift=0., + step=1, nsum=1, cradius=params.cradius, + threshold=params.threshold, nlost=100, newaps=no, + refit=no, trace=no, override=no, addfeatures=no, + database=database, plotfile=plotfile, + logfiles=logfile, verbose=verbose1) +end diff --git a/noao/imred/src/fibers/doalign.par b/noao/imred/src/fibers/doalign.par new file mode 100644 index 00000000..0ddd375f --- /dev/null +++ b/noao/imred/src/fibers/doalign.par @@ -0,0 +1,7 @@ +spec,f,a,"",,, +specms,f,a,"",,, +align,f,a,"",,, +table,f,a,"",,, +logfile,f,a,"",,, +batch,b,a,,,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/doarcs.cl b/noao/imred/src/fibers/doarcs.cl new file mode 100644 index 00000000..bd4e3748 --- /dev/null +++ b/noao/imred/src/fibers/doarcs.cl @@ -0,0 +1,264 @@ +# DOARCS -- Determine dispersion relation for spectrum based on reference arcs. +# This procedure is complicated by: +# 1. The need to reextract arcs if the objects spectra are being +# recentered or retraced. +# 2. The use of shift spectra to track shifts in the dispersion from +# the reference arc spectrum. +# 3. The use of multiple exposures to correct for illumination problems +# in taking the arcs. + +procedure doarcs (spec, response, arcref1, arcref2, extn, arcreplace, apidtable, + arcaps, arcbeams, savearcs, reextract, arcap, logfile, batch, done) + +file spec +file response +file arcref1 +file arcref2 +string extn +file arcreplace +file apidtable +string arcaps +string arcbeams +bool savearcs +bool reextract +bool arcap +file logfile +bool batch +file done + +struct *fd + +begin + string imtype + int i, j, k, n + file temp, arc1, arc2, str1, str2, arctype, apref, arc, arcms + bool verbose1 + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + n = strlen (imtype) + + temp = mktemp ("tmp$iraf") + + if (batch) + verbose1 = no + else + verbose1 = verbose + + for (j=1; j<=2; j+=1) { + # The reference spectra refer initially to the 2D image. At the + # end we will reset them to refer to the 1D spectra. + + hselect (spec, "refspec"//j, yes, > temp) + fd = temp + k = fscan (fd, arc1, str1) + fd = ""; delete (temp, verify=no) + if (k < 1) + break + + # Strip possible image extension. + i = strlen (arc1) + if (i > n && substr (arc1, i-n+1, i) == imtype) + arc1 = substr (arc1, 1, i-n) + + # Set extraction output and aperture reference depending on whether + # the arcs are to be rextracted using recentered or retraced object + # apertures. + + if (arcap && + (apscript.ansrecenter=="yes" || apscript.anstrace=="yes" || + apscript.ansrecenter=="YES" || apscript.anstrace=="YES")) { + arc2 = spec // arc1 // ".ms" + apref = spec + if (access (arc2//imtype)) + imdelete (arc2//imtype, verify=no) + delete (database//"/id"//arc2//"*", verify = no) + } else { + arc2 = arc1 // extn + apref = apscript.references + if (reextract && access (arc2//imtype)) { + if (arc2 != arcref1 // extn && arc2 != arcref2 // extn) { + if (access (done)) { + fd = done + while (fscan (fd, arcms) != EOF) + if (arcms == arc2) + break + fd = "" + } else + arcms = "" + if (arcms != arc2) + imdelete (arc2, verify=no) + } + } + } + + # SHIFT arcs are reidentified with only a shift. + # HENEAR arcs are reidentified using the user refit option. + # Also internal arcs are checked if HENEAR. + + hselect (arc1, "arctype", yes, > temp) + fd = temp + i = fscan (fd, arctype) + fd = ""; delete (temp, verify=no) + + # Extract and determine dispersion function if necessary. + if (!access (arc2//imtype)) { + delete (database//"/id"//arc2//"*", verify = no) + if (!batch) + print ("Extract and reidentify arc spectrum ", arc1) + print ("Extract and reidentify arc spectrum ", arc1, >> logfile) + apscript (arc1, output=arc2, references=apref, + ansrecenter="NO", ansresize="NO", ansedit="NO", + anstrace="NO", nsubaps=params.nsubaps, background="none", + clean=no, weights="none", verbose=verbose1) + sapertures (arc2, 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 (arc2, "/", response, arc2, 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 (arc2, "/", temp, arc2, 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) + } + } + print (arc2, >> done) + + if (arctype == "shift") { + reidentify (arcref2//extn, arc2, + interactive=no, section="middle line", shift=0., + step=1, nsum=1, cradius=params.cradius, + threshold=params.threshold, nlost=100, newaps=no, + refit=no, trace=no, override=no, addfeatures=no, + database=database, plotfile=plotfile, + logfiles=logfile, verbose=verbose1) + } else { + if (arcreplace != "") { + fd = arcreplace + while (fscan (fd, arc, arcms, str2) != EOF) { + i = strlen (arc) + if (i > n && substr (arc, i-n+1, i) == imtype) + arc = substr (arc, 1, i-n) + if (arc != arc1) + next + arc = arcms + if (i > n && substr (arc, i-n+1, i) == imtype) + arc = substr (arc, 1, i-n) + arcms = arc // extn // imtype + + if (access (arcms)) + imdelete (arcms, verify=no) + + if (!batch) + print ("Extract arc spectrum ", arc) + print ("Extract arc spectrum ", arc, >> logfile) + apscript (arc, references=apref, + ansrecenter="NO", ansresize="NO", ansedit="NO", + anstrace="NO", nsubaps=params.nsubaps, + background="none", clean=no, + weights="none", verbose=verbose1) + sapertures (arcms, 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 (arcms, "/", response, arcfms, + 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 (arcms, "/", temp, arcfms, + 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) + } + } + scopy (arcms, arc2, w1=INDEF, w2=INDEF, + apertures=str2, bands="", beams="", + apmodulus=1000, offset=0, format="multispec", + clobber=yes, merge=yes, renumber=no, + verbose=yes, >> logfile) + imdelete (arcms, verify=no) + } + fd = "" + } + reidentify (arcref1//extn, arc2, + interactive=!batch, section="middle line", + shift=0., step=1, nsum=1, cradius=params.cradius, + threshold=params.threshold, nlost=100, + refit=params.refit, trace=no, override=no, + addfeatures=params.addfeatures, + coordlist=params.coordlist, match=params.match, + maxfeatures=50, minsep=2., database=database, + plotfile=plotfile, logfiles=logfile, + verbose=verbose1) + } + + # If not reextracting arcs based on object apertures + # then save the extracted arc to avoid doing it again. + + if (arc1//extn != arc2) + imdelete (arc2, verify=no) + } + + # Set the REFSPEC parameters for multispec spectrum. + if (k == 1) + hedit (spec//".ms", "refspec"//j, arc2, add=yes, + verify=no, show=no, update=yes) + else + hedit (spec//".ms", "refspec"//j, arc2//" "//str1, + add=yes, verify=no, show=no, update=yes) + + # Check for arc fibers in object spectra. + if (arctype != "shift" && (arcaps != "" || arcbeams != "")) { + scopy (spec//".ms", spec//"arc.ms", w1=INDEF, w2=INDEF, + apertures=arcaps, bands="", beams=arcbeams, apmodulus=1000, + offset=0, format="multispec", clobber=yes, merge=no, + renumber=no, verbose=no, >& "dev$null") + if (access (spec//"arc.ms"//imtype)) { + if (!batch) + print ("Reidentify arc fibers in ", spec, + " with respect to ", arc1) + print ("Reidentify arc fibers in ", spec, + " with respect to ", arc1, >> logfile) + delete (database//"/id"//spec//"arc.ms*", verify = no) + reidentify (arc2, spec//"arc.ms", interactive=no, + section="middle line", shift=0., step=1, nsum=1, + cradius=params.cradius, threshold=params.threshold, + nlost=100, refit=no, trace=no, override=no, + addfeatures=no, database=database, + plotfile=plotfile, logfiles=logfile, + verbose=verbose1) + imdelete (spec//"arc.ms", verify=no) + hedit (spec//".ms", "refshft"//j, spec//"arc.ms interp", + add=yes, verify=no, show=no, update=yes) + if (!savearcs) + scopy (spec//".ms", "", w1=INDEF, w2=INDEF, + apertures="!"//arcaps, bands="", beams=arcbeams, + apmodulus=1000, offset=0, format="multispec", + clobber=yes, merge=no, renumber=no, + verbose=yes, >> logfile) + } + } + } +end diff --git a/noao/imred/src/fibers/doarcs.par b/noao/imred/src/fibers/doarcs.par new file mode 100644 index 00000000..a93b16c6 --- /dev/null +++ b/noao/imred/src/fibers/doarcs.par @@ -0,0 +1,17 @@ +spec,f,a,"",,, +response,f,a,"",,, +arcref1,f,a,"",,, +arcref2,f,a,"",,, +extn,s,a,,,, +arcreplace,f,a,"",,, +apidtable,f,a,"",,, +arcaps,s,a,,,, +arcbeams,s,a,,,, +savearcs,b,a,,,, +reextract,b,a,,,, +arcap,b,a,,,, +logfile,f,a,"",,, +batch,b,a,,,, +done,f,a,"",,, +fd,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/fibresponse.cl b/noao/imred/src/fibers/fibresponse.cl new file mode 100644 index 00000000..2379fcc4 --- /dev/null +++ b/noao/imred/src/fibers/fibresponse.cl @@ -0,0 +1,261 @@ +# FIBRESPONSE -- Make a aperture response spectrum using a flat field +# and a throughput file or image. + +procedure fibresponse (flat, throughput, apreference, response) + +string flat {prompt="Flat field spectrum"} +string throughput {prompt="Throughput file or image"} +string apreference {prompt="Aperture reference spectrum"} +string response {prompt="Response spectrum"} + +bool recenter = no {prompt="Recenter apertures?"} +bool edit = yes {prompt="Edit/review apertures?"} +bool trace = no {prompt="Trace spectra?"} +bool clean = no {prompt="Detect and replace bad pixels?"} +bool fitflat = no {prompt="Fit and ratio flat field spectrum?"} +bool interactive = yes {prompt="Fit flat field interactively?"} +string function = "spline3" {prompt="Fitting function", + enum="spline3|legendre|chebyshev|spline1"} +int order = 20 {prompt="Fitting function order", min=1} + +begin + string imtype, mstype + file flat2d, skyflat2d, apref, resp + file temp1, temp2, log1, log2 + int n, ap, naxis + real respval + struct err + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + mstype = ".ms" // imtype + n = strlen (imtype) + + flat2d = flat + skyflat2d = throughput + apref = apreference + resp = response + temp1 = mktemp ("tmp") + temp2 = mktemp ("tmp") + + # Check required input and output. + if (resp == "" || resp == flat2d || resp == skyflat2d) + error (1, "Bad response image name") + if (flat2d == "" && skyflat2d == "") + error (1, "No flat field or throughput specified") + + if (flat2d != "") { + i = strlen (flat2d) + if (i > n && substr (flat2d, i-n+1, i) == imtype) + flat2d = substr (flat2d, 1, i-n) + if (!access (flat2d // imtype)) { + printf ("Flat field spectrum not found - %s%s\n", + flat2d, imtype) | scan (err) + error (1, err // "\nCheck settting of imtype") + } + } + if (skyflat2d != "") { + i = strlen (skyflat2d) + if (i > n && substr (skyflat2d, i-n+1, i) == imtype) + skyflat2d = substr (skyflat2d, 1, i-n) + if (!access (skyflat2d // imtype)) { + if (!access (skyflat2d)) { + printf ("Throughput file or image not found - %s%s\n", + skyflat2d, imtype) | scan (err) + error (1, err // "\nCheck settting of imtype") + } + + if (flat2d == "") { + i = strlen (apref) + if (i > n && substr (apref, i-n+1, i) == imtype) + apref = substr (apref, 1, i-n) + if (!access (apref // imtype)) + error (1, "Aperture reference image required") + } + } + } + + # Set logging + tee.append = yes + if (logfile == "") + log1 = "dev$null" + else + log1 = logfile + if (verbose) + log2 = "STDOUT" + else + log2 = "dev$null" + + # Initialize APSCRIPT + apscript.references = apref + 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.ansextract = "YES" + + # If using a flat field extract it if necessary and possibly fit it + # and ratio the individual apertures by an overall smooth function + + if (flat2d != "") { + if (!access (flat2d // mstype)) { + print ("Extract flat field ", flat2d) | tee (log1) + if (flat2d != apref) + apscript (flat2d, output=resp, background="none", + clean=clean, extras=no) + else + apscript (flat2d, output=resp, ansrecenter="NO", + ansrecenter="NO", ansresize="NO", ansedit="NO", + anstrace="NO", background="none", + clean=clean, extras=no) + } else + imcopy (flat2d//".ms", resp, verbose=no) + + if (fitflat) { + print ("Fit and ratio flat field ", flat2d) | tee (log1) + blkavg (resp, temp1, option="average", b1=1, b2=10000) + imcopy (temp1//"[*,1]", temp1, verbose=no) + fit1d (temp1, temp1, "fit", axis=1, interactive=interactive, + sample="*", naverage=1, function=function, order=order, + low_reject=0., high_reject=0., niterate=1, grow=0., + graphics="stdgraph") + sarith (resp, "/", temp1, resp, 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 (temp1, verify=no) + } + } + + # If using a throughput image extract it if necessary. + # Apply it to the flat field if given and otherwise only + # compute the throughput through each aperture. + + if (skyflat2d != "") { + if (access (skyflat2d // imtype)) { + if (!access (skyflat2d // mstype)) { + print ("Extract throughput image ", skyflat2d) | tee (log1) + apscript (skyflat2d, output=temp1, background="none", + clean=clean, extras=no) + temp2 = temp1 + } else + temp2 = skyflat2d // ".ms" + + if (flat2d != "") { + print ("Correct flat field to throughput image") | + tee (log1) + sarith (temp2, "/", resp, temp1, 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) + fit1d (temp1, temp1, type="fit", axis=1, + interactive=no, sample="*", naverage=1, + function="legendre", order=1, niterate=0) + sarith (resp, "*", temp1, resp, 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 (temp1, verify=no) + } else { + print ("Compute aperture throughput from image") | + tee (log1) + fit1d (temp2, resp, type="fit", axis=1, + interactive=no, sample="*", naverage=1, + function="legendre", order=1, niterate=0) + if (temp2 == temp1) + imdelete (temp2, verify=no) + } + + # If a flat field and throughput file are used scale the average + # flat field in each aperture to those values + + } else if (flat2d != "") { + print ("Correct flat field with throughput file ", skyflat2d) | + tee (log1) + fit1d (resp, resp, type="ratio", axis=1, + interactive=no, sample="*", naverage=1, + function="legendre", order=1, niterate=0) + + list = skyflat2d + while (fscan (list, ap, respval) != EOF) { + sarith (resp, "*", respval, resp, w1=INDEF, w2=INDEF, + apertures=ap, bands="", beams="", apmodulus=0, + reverse=no, ignoreaps=no, format="multispec", + renumber=no, offset=0, clobber=yes, merge=yes, + errval=0., verbose=no, >& "dev$null") + } + list = "" + + # If only a throughput file is given create the response from the + # aperture reference and set the aperture response to the specified + # values. + + } else { + print ("Set aperture throughput using ", skyflat2d) | tee (log1) + if (!access (apref // mstype)) { + apscript (apref, output=resp, ansrecenter="NO", + ansrecenter="NO", ansresize="NO", ansedit="NO", + anstrace="NO", background="none", + clean=no, extras=no) + sarith (resp, "replace", "0", resp, w1=INDEF, + w2=INDEF, apertures="", bands="", beams="", apmodulus=0, + reverse=no, ignoreaps=no, format="multispec", + renumber=no, offset=0, clobber=yes, merge=yes, + errval=0., verbose=no) + } else + sarith (apref//".ms", "replace", "0", resp, w1=INDEF, + w2=INDEF, apertures="", bands="", beams="", apmodulus=0, + reverse=no, ignoreaps=no, format="multispec", + renumber=no, offset=0, clobber=yes, merge=yes, + errval=0., verbose=no) + + list = skyflat2d + while (fscan (list, ap, respval) != EOF) { + sarith (resp, "replace", respval, resp, w1=INDEF, w2=INDEF, + apertures=ap, bands="", beams="", apmodulus=0, + reverse=no, ignoreaps=no, format="multispec", + renumber=no, offset=0, clobber=yes, merge=yes, + errval=0., verbose=no) + } + list = "" + } + } + + # The final response is normalized to overall unit mean and the + # average aperture response is printed. + + print ("Create the normalized response ", resp) | tee (log1) + bscale (resp, resp, bzero="0.", bscale="mean", section="", + step=1, upper=INDEF, lower=INDEF, verbose=yes) | tee (log1, >log2) + blkavg (resp, temp1, option="average", b1=10000, b2=1) + print ("Average aperture response:") | tee (log1, >log2) + naxis = 5 + #hselect (temp1, "naxis", yes) | scan (naxis) + #if (naxis == 1) + # listpixels (temp1) | tee (log1, >log2) + #else + # listpixels (temp1//"[1,*]") | tee (log1, >log2) + + hselect (temp1, "naxis", yes, > temp2) + list = temp2; ap = fscan (list, naxis) + if (naxis == 1) + listpixels (temp1) | tee (log1, >log2) + else + listpixels (temp1//"[1,*]") | tee (log1, >log2) + list = ""; delete (temp2, verify=no) + + imdelete (temp1, verify=no) +end diff --git a/noao/imred/src/fibers/fibresponse.par b/noao/imred/src/fibers/fibresponse.par new file mode 100644 index 00000000..9a59eb46 --- /dev/null +++ b/noao/imred/src/fibers/fibresponse.par @@ -0,0 +1,13 @@ +flat,s,a,,,,"Flat field spectrum" +throughput,s,a,,,,"Throughput file or image" +apreference,s,a,,,,"Aperture reference spectrum" +response,s,a,,,,"Response spectrum" +recenter,b,h,no,,,"Recenter apertures?" +edit,b,h,yes,,,"Edit/review apertures?" +trace,b,h,no,,,"Trace spectra?" +clean,b,h,no,,,"Detect and replace bad pixels?" +fitflat,b,h,no,,,"Fit and ratio flat field spectrum?" +interactive,b,h,yes,,,"Fit flat field interactively?" +function,s,h,"spline3",spline3|legendre|chebyshev|spline1,,"Fitting function" +order,i,h,20,1,,"Fitting function order" +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/getspec.cl b/noao/imred/src/fibers/getspec.cl new file mode 100644 index 00000000..84ce9a1c --- /dev/null +++ b/noao/imred/src/fibers/getspec.cl @@ -0,0 +1,49 @@ +# GETSPEC -- Get spectra which are processed but are not extracted. +# Strip the imtype extension. + +procedure getspec (images, output) + +string images {prompt="List of images"} +file output {prompt="Output file of images"} +bool ccdproc {prompt="Add CCDPROC keyword and continue?", + mode="q"} +struct *fd + +begin + string imtype, temp, image, system="" + int n, n1 + + imtype = "." // envget ("imtype") + n = stridx (",", imtype) + if (n > 0) + imtype = substr (imtype, 1, n-1) + n1 = strlen (imtype) + + # Initialize files + set clobber=yes + sleep (> output) + set clobber=no + + temp = mktemp ("tmp$iraf") + sections (images, option="fullname", > temp) + fd = temp + while (fscan (fd, image) != EOF) { + hselect (image, "ccdproc", yes) | scan (system) + if (nscan() == 0) { + printf ("%s: CCDPROC keyword not found.\n", image) + printf (" Either run CCDPROC or add CCDPROC keyword with HEDIT.\n") + if (!ccdproc) + error (1, "Exit") + hedit (image, "ccdproc", "DONE", add=yes, update=yes, + verify=no, show=no) + } + hselect (image, "wat0_001", yes) | scanf ("system=%s", system) + if (system=="equispec" || system=="multispec") + next + n = strlen (image) + if (n > n1 && substr (image, n-n1+1, n) == imtype) + image = substr (image, 1, n-n1) + print (image, >> output) + } + fd = ""; delete (temp, verify=no) +end diff --git a/noao/imred/src/fibers/getspec.par b/noao/imred/src/fibers/getspec.par new file mode 100644 index 00000000..e676c18e --- /dev/null +++ b/noao/imred/src/fibers/getspec.par @@ -0,0 +1,5 @@ +images,s,a,,,,"List of images" +output,f,a,,,,"Output file of images" +ccdproc,b,q,,,,"Add CCDPROC keyword and continue?" +fd,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/listonly.cl b/noao/imred/src/fibers/listonly.cl new file mode 100644 index 00000000..ac215892 --- /dev/null +++ b/noao/imred/src/fibers/listonly.cl @@ -0,0 +1,237 @@ +# LISTONLY -- List processing to be done. +# +# This follows pretty much the same logic as the full procedure but doesn't +# do anything but list the operations. + +procedure listonly (objects, apidtable, apref, flat, throughput, arcs1, arcs2, + scattered, dispcor, skysubtract, redo, update) + +string objects = "" {prompt="List of object spectra"} +file apidtable = "" {prompt="Aperture ID table"} +file apref = "" {prompt="Aperture reference spectrum"} +file flat = "" {prompt="Flat field spectrum"} +file throughput = "" {prompt="Throughput file or image"} +string arcs1 = "" {prompt="List of arc spectra"} +string arcs2 = "" {prompt="List of shift arc spectra"} + +bool scattered {prompt="Subtract scattered light?"} +bool dispcor {prompt="Dispersion correct spectra?"} +bool skysubtract {prompt="Subtract sky?"} +bool redo = no {prompt="Redo operations if previously done?"} +bool update = yes {prompt="Update spectra if cal data changes?"} + +struct *fd1 +struct *fd2 + +begin + string imtype, mstype, extn + string spec, arcref1, arcref2 + string specms, arcref1ms, arcref2ms, response + string objs, temp, done, str + bool reextract, newaps, newresp, newdisp, scat, extract, disp, sky + int i, j, n, dc + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + mstype = ".ms" // imtype + n = strlen (imtype) + + objs = mktemp ("tmp$iraf") + temp = mktemp ("tmp$iraf") + done = mktemp ("tmp$iraf") + + if (apidtable != "") { + j = strlen (apidtable) + for (i=1; i<=j && substr(apidtable,i,i)==" "; i+=1); + apidtable = substr (apidtable, i, j) + } + 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" + } + + newaps = no + newresp = no + newdisp = no + + i = strlen (apref) + if (i > n && substr (apref, i-n+1, i) == imtype) + apref = substr (apref, 1, i-n) + + reextract = redo + if (reextract || !access (database // "/ap" // apref // extn)) { + print ("Set reference aperture for ", apref) + newaps = yes + } + + i = strlen (flat) + if (i > n && substr (flat, i-n+1, i) == imtype) + flat = substr (flat, 1, i-n) + if (flat != "") { + scat = no + if (scattered) { + if (redo && access (flat//"noscat"//imtype)) + hselect (flat//"noscat", "apscatte", yes) | scan (str) + else + hselect (flat, "apscatte", yes) | scan (str) + if (nscan() == 0) + scat = yes + } + if (scat) + print ("Subtract scattered light from ", flat) + } + + spec = throughput + i = strlen (spec) + if (i > n && substr (spec, i-n+1, i) == imtype) + spec = substr (spec, 1, i-n) + if (spec != "") { + scat = no + if (scattered) { + if (redo && access (flat//"noscat"//imtype)) + hselect (flat//"noscat", "apscatte", yes) | scan (str) + else + hselect (flat, "apscatte", yes) | scan (str) + if (nscan() == 0) + scat = yes + } + if (scat) + print ("Subtract scattered light from ", spec) + } + + 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) + newresp = yes + } + } + + if (dispcor) { + getspec (arcs1, temp) + fd1 = temp + if (fscan (fd1, arcref1) == EOF) + error (1, "No reference arcs") + fd1 = ""; delete (temp, verify=no) + arcref1ms = arcref1 // extn + + getspec (arcs2, temp) + fd1 = temp + if (fscan (fd1, arcref2) == EOF) + arcref2 = "" + fd1 = ""; delete (temp, verify=no) + arcref2ms = arcref2 // extn + + reextract = redo || (update && newaps) + if (reextract || !access (arcref1ms//imtype)) { + print ("Extract arc reference image ", arcref1) + print ("Determine dispersion solution for ", arcref1) + newdisp = yes + } else { + hselect (arcref1ms, "dclog1", yes, > temp) + fd1 = temp + dc = -1 + i = fscan (fd1, dc) + fd1 = ""; delete (temp, verify=no) + if (i < 1) { + print ("Determine dispersion solution for ", arcref1) + newdisp = yes + } + } + print (arcref1, > done) + + if (arcref2 != "") { + if (reextract || !access (arcref2ms//imtype) || newdisp) { + print ("Extract shift arc reference image ", arcref2) + print ("Determine dispersion solution for ", arcref2) + } else { + hselect (arcref2ms, "dclog1", yes, > temp) + fd1 = temp + dc = -1 + i = fscan (fd1, dc) + fd1 = ""; delete (temp, verify=no) + if (i < 1) + print ("Determine dispersion solution for ", arcref2) + } + } + print (arcref2, >> done) + } + + reextract = redo || (update && (newaps || newresp || newdisp)) + getspec (objects, objs) + fd1 = objs + while (fscan (fd1, spec) != EOF) { + if (access (done)) { + fd2 = done + while (fscan (fd2, specms) != EOF) + if (spec == specms) + break + if (spec == specms) + next + fd2 = "" + } + + specms = spec // mstype + + scat = no + extract = no + disp = no + sky = no + if (scattered) { + if (redo && access (spec//"noscat"//imtype)) + hselect (spec//"noscat", "apscatte", yes) | scan (str) + else + hselect (spec, "apscatte", yes) | scan (str) + if (nscan() == 0) + scat = yes + } + if (reextract || !access (specms) || (redo && scat)) + extract = yes + else { + hselect (specms, "dclog1", yes) | scan (str) + if (nscan() == 0) + disp = yes + else + extract = update && newdisp + hselect (specms, "skysub", yes) | scan (str) + if (nscan() == 0) + sky = skysubtract + } + + if (extract) { + disp = dispcor + sky = skysubtract + } + + if (scat) + print ("Subtract scattered light from ", spec) + if (extract) + print ("Extract object spectrum ", spec) + if (disp) + print ("Dispersion correct ", spec) + if (sky) + print ("Sky subtract ", spec) + } + fd1 = ""; delete (objs, verify=no) + + if (access (done)) + delete (done, verify=no) +end diff --git a/noao/imred/src/fibers/listonly.par b/noao/imred/src/fibers/listonly.par new file mode 100644 index 00000000..aa52691e --- /dev/null +++ b/noao/imred/src/fibers/listonly.par @@ -0,0 +1,15 @@ +objects,s,a,"",,,"List of object spectra" +apidtable,f,a,"",,,"Aperture ID table" +apref,f,a,"",,,"Aperture reference spectrum" +flat,f,a,"",,,"Flat field spectrum" +throughput,f,a,"",,,"Throughput file or image" +arcs1,s,a,"",,,"List of arc spectra" +arcs2,s,a,"",,,"List of shift arc spectra" +scattered,b,a,,,,"Subtract scattered light?" +dispcor,b,a,,,,"Dispersion correct spectra?" +skysubtract,b,a,,,,"Subtract sky?" +redo,b,a,no,,,"Redo operations if previously done?" +update,b,a,yes,,,"Update spectra if cal data changes?" +fd1,*struct,h,"",,, +fd2,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/mkfibers.cl b/noao/imred/src/fibers/mkfibers.cl new file mode 100644 index 00000000..b7ffdac2 --- /dev/null +++ b/noao/imred/src/fibers/mkfibers.cl @@ -0,0 +1,167 @@ +# MKFIBERS - Make multifiber examples + +procedure mkfibers (image) + +file image {prompt="Image name"} +string type="object" {prompt="Object type", + enum="object|objnosky|sky|flat|henear|ehenear|ohenear|mercury"} +file fibers="" {prompt="Fiber data file"} +string title="Multifiber artificial image" {prompt="Title"} +file header="artdata$stdhdr.dat" {prompt="Header keyword file"} +int ncols=400 {prompt="Number of columns"} +int nlines=512 {prompt="Number of lines"} +real wstart=4210. {prompt="Starting wavelength"} +real wend=7362. {prompt="Ending wavelength"} +int seed=1 {prompt="Noise seed"} + +begin + int i, ap, beam + real ar + file out, obj, sky, arc, dat + string htype, imtype + + out = image + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + if (access (out) || access (out//imtype)) + return + + print ("Creating image ", out, " ...") + + obj = mktemp ("art") + sky = mktemp ("art") + arc = mktemp ("art") + dat = mktemp ("art") + + list = fibers + if (type == "object") { # Object spectrum + sky + htype = "object" + mk1dspec (obj, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=1000., slope=0., + temperature=7000., lines="", nlines=50, peak=-0.5, + profile="gaussian", gfwhm=24, seed=2, comments=no, header="") + mk1dspec (sky, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=1000., slope=0., + temperature=5800., lines="", nlines=20, peak=1., + profile="gaussian", gfwhm=12, seed=1, comments=no, header="") + imarith (obj, "+", sky, obj, verbose=no, noact=no) + mk1dspec (arc, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=0.8, slope=0., + temperature=0., lines="mkexamples$henear2.dat", + profile="gaussian", gfwhm=14, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, continuum=20, + slope=0., temperature=0., lines="", nlines=0, + comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) { + if (beam == 0) + print (sky, " ", line, >> dat) + else if (beam == 1) + print (obj, " ", line, >> dat) + else if (beam == 2) + print (arc, " ", line, >> dat) + } + } else if (type == "objnosky") { # Object spectrum + htype = "object" + mk1dspec (obj, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=1000., slope=0., + temperature=7000., lines="", nlines=50, peak=-0.5, + profile="gaussian", gfwhm=24, seed=2, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=0.8, slope=0., + temperature=0., lines="mkexamples$henear2.dat", + profile="gaussian", gfwhm=14, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, continuum=20, + slope=0., temperature=0., lines="", nlines=0, + comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) { + if (beam == 1) + print (obj, " ", line, >> dat) + else if (beam == 2) + print (arc, " ", line, >> dat) + } + } else if (type == "sky") { # Sky only + htype = "object" + mk1dspec (sky, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=1000., slope=0., + temperature=5800., lines="", nlines=20, peak=1., + profile="gaussian", gfwhm=12, seed=1, comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) + print (sky, " ", line, >> dat) + } else if (type == "flat") { # Flat field + htype = "flat" + mk1dspec (obj, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=10000., slope=0., + temperature=8000., lines="", nlines=0, comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) + print (obj, " ", line, >> dat) + } else if (type == "henear") { # HE-NE-AR + htype = "comp" + mk1dspec (arc, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=0.8, slope=0., + temperature=0., lines="mkexamples$henear2.dat", + profile="gaussian", gfwhm=14, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, continuum=20, + slope=0., temperature=0., lines="", nlines=0, + comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) + print (arc, " ", line, >> dat) + } else if (type == "ehenear") { # HE-NE-AR Even fibers + htype = "comp" + mk1dspec (arc, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=0.8, slope=0., + temperature=0., lines="mkexamples$henear2.dat", + profile="gaussian", gfwhm=14, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, continuum=20, + slope=0., temperature=0., lines="", nlines=0, + comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) { + if (mod (ap, 2) == 0) { + print (arc, " ", line, >> dat) + } + } + } else if (type == "ohenear") { # HE-NE-AR Odd fibers + htype = "comp" + mk1dspec (arc, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=0.8, slope=0., + temperature=0., lines="mkexamples$henear2.dat", + profile="gaussian", gfwhm=14, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, continuum=20, + slope=0., temperature=0., lines="", nlines=0, + comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) { + if (mod (ap, 2) == 1) { + print (arc, " ", line, >> dat) + } + } + } else if (type == "mercury") { # Emission lines + htype = "comp" + mk1dspec (arc, output="", ap=1, rv=0., z=no, ncols=nlines, naps=1, + wstart=wstart, wend=wend, continuum=0., slope=0., + temperature=0., lines="", nlines=30, peak=10000., + profile="gaussian", gfwhm=7, seed=i, comments=no, header="") + mk1dspec (arc, output="", ap=1, rv=0., z=no, continuum=20, + slope=0., temperature=0., lines="", nlines=0, + comments=no, header="") + while (fscan (list, ap, beam, line) != EOF) { + print (arc, " ", line, >> dat) + } + } + list = "" + + mk2dspec (out, output="", model=dat, ncols=ncols, nlines=nlines, + title=title, header=header, comments=no) + hedit (out, "imagetyp", htype, update=yes, add=no, delete=no, + show=no, verify=no) + + + mknoise (out, output="", background=0., gain=1., rdnoise=3., + poisson=yes, seed=seed, cosrays="", ncosrays=0, energy=30000., + radius=0.5, ar=1., pa=0., comments=no) + + imdelete (obj, verify=no, >& "dev$null") + imdelete (sky, verify=no, >& "dev$null") + imdelete (arc, verify=no, >& "dev$null") + delete (dat, verify=no, >& "dev$null") +end diff --git a/noao/imred/src/fibers/mkfibers.par b/noao/imred/src/fibers/mkfibers.par new file mode 100644 index 00000000..9ba33353 --- /dev/null +++ b/noao/imred/src/fibers/mkfibers.par @@ -0,0 +1,11 @@ +image,f,a,"",,,"Image name" +type,s,h,"object",object|objnosky|sky|flat|henear|ehenear|ohenear|mercury,,"Object type" +fibers,f,h,"",,,"Fiber data file" +title,s,h,"Multifiber artificial image",,,"Title" +header,f,h,"artdata$stdhdr.dat",,,"Header keyword file" +ncols,i,h,400,,,"Number of columns" +nlines,i,h,512,,,"Number of lines" +wstart,r,h,4210.,,,"Starting wavelength" +wend,r,h,7362.,,,"Ending wavelength" +seed,i,h,1,,,"Noise seed" +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/params.par b/noao/imred/src/fibers/params.par new file mode 100644 index 00000000..bf41e8ec --- /dev/null +++ b/noao/imred/src/fibers/params.par @@ -0,0 +1,75 @@ +line,i,h,INDEF,,,Default dispersion line +nsum,i,h,10,,,Number of dispersion lines to sum or median +width,r,h,,,,Width of profiles +minsep,r,h,1.,,,Minimum separation between spectra +maxsep,r,h,100000.,,,Maximum separation between spectra +order,s,h,"increasing","increasing|decreasing",,"Order of apertures" +extras,b,h,no,,,"Extract sky, sigma, etc.? + +-- DEFAULT APERTURE LIMITS --" +lower,r,h,-3.,,,"Lower aperture limit relative to center" +upper,r,h,3.,,,"Upper aperture limit relative to center + +-- AUTOMATIC APERTURE RESIZING PARAMETERS --" +ylevel,r,h,0.05,,,"Fraction of peak or intensity for resizing" +peak,b,h,yes,,,"Is ylevel a fraction of the peak?" +bkg,b,h,yes,,,"Subtract background for resizing?" +avglimits,b,h,no,,,"Average limits over all apertures? + +-- TRACE PARAMETERS --" +t_step,i,h,10,,,"Tracing step" +t_function,s,h,"spline3","chebyshev|legendre|spline1|spline3",,"Trace fitting function" +t_order,i,h,2,,,"Trace fitting function order" +t_niterate,i,h,1,0,,"Trace rejection iterations" +t_low,r,h,3.,0.,,"Trace lower rejection sigma" +t_high,r,h,3.,0.,,"Trace upper rejection sigma + +-- SCATTERED LIGHT PARAMETERS --" +buffer,r,h,1.,0.,,Buffer distance from apertures +apscat1,pset,h,"",,,Fitting parameters across the dispersion +apscat2,pset,h,"",,,"Fitting parameters along the dispersion + +-- APERTURE EXTRACTION PARAMETERS --" +weights,s,h,"none","none|variance",,Extraction weights (none|variance) +pfit,s,h,"fit1d","fit1d|fit2d",,Profile fitting algorithm (fit1d|fit2d) +readnoise,s,h,0.,,,Read out noise sigma (photons) +gain,s,h,1.,,,Photon gain (photons/data number) +lsigma,r,h,3.,,,Lower rejection threshold +usigma,r,h,3.,,,Upper rejection threshold +nsubaps,i,h,1,1,,"Number of subapertures + +-- FLAT FIELD FUNCTION FITTING PARAMETERS --" +f_interactive,b,h,yes,,,"Fit flat field interactively?" +f_function,s,h,"spline3",spline3|legendre|chebyshev|spline1,,"Fitting function" +f_order,i,h,20,1,,"Fitting function order + +-- ARC DISPERSION FUNCTION PARAMETERS --" +threshold,r,h,10.,0.,,"Minimum line contrast threshold" +coordlist,f,h,linelist$idhenear.dat,,,"Line list" +match,r,h,-3.,,,"Line list matching limit in Angstroms" +fwidth,r,h,4.,,,"Arc line widths in pixels" +cradius,r,h,10.,,,Centering radius in pixels +i_function,s,h,"chebyshev","legendre|chebyshev|spline1|spline3",,"Coordinate function" +i_order,i,h,3,1,,"Order of dispersion function" +i_niterate,i,h,2,0,,"Rejection iterations" +i_low,r,h,3.,0.,,"Lower rejection sigma" +i_high,r,h,3.,0.,,"Upper rejection sigma" +refit,b,h,yes,,,"Refit coordinate function when reidentifying?" +addfeatures,b,h,no,,,"Add features when reidentifying? + +-- AUTOMATIC ARC ASSIGNMENT PARAMETERS --" +select,s,h,"interp",,,"Selection method for reference spectra" +sort,s,h,"jd",,,"Sort key" +group,s,h,"ljd",,,"Group key" +time,b,h,no,,,"Is sort key a time?" +timewrap,r,h,17.,0.,24.,"Time wrap point for time sorting + +-- DISPERSION CORRECTION PARAMETERS --" +linearize,b,h,yes,,,Linearize (interpolate) spectra? +log,b,h,no,,,"Logarithmic wavelength scale?" +flux,b,h,yes,,,"Conserve flux? + +-- SKY SUBTRACTION PARAMETERS --" +combine,s,h,"average","average|median",,Type of combine operation +reject,s,h,"avsigclip","none|minmax|avsigclip",,"Sky rejection option" +scale,s,h,"none","none|mode|median|mean",,"Sky scaling option" 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 diff --git a/noao/imred/src/fibers/proc.par b/noao/imred/src/fibers/proc.par new file mode 100644 index 00000000..1b3961ea --- /dev/null +++ b/noao/imred/src/fibers/proc.par @@ -0,0 +1,52 @@ +objects,s,a,,,,"List of object spectra" +apref,f,a,"",,,"Aperture reference spectrum" +flat,f,a,"",,,"Flat field spectrum" +throughput,f,a,"",,,"Throughput file or image (optional)" +arcs1,s,a,,,,"List of arc spectra" +arcs2,s,a,,,,"List of shift arc spectra" +arcreplace,f,a,"",,,"Special aperture replacements" +arctable,f,a,"",,,"Arc assignment table (optional) +" +fibers,i,a,,,,"Number of fibers" +apidtable,f,a,"",,,"Aperture identifications" +crval,s,a,INDEF,,,"Approximate wavelength" +cdelt,s,a,INDEF,,,"Approximate dispersion" +objaps,s,a,,,,"Object apertures" +skyaps,s,a,,,,"Sky apertures" +arcaps,s,a,,,,"Arc apertures" +objbeams,s,a,,,,"Object beam numbers" +skybeams,s,a,,,,"Sky beam numbers" +arcbeams,s,a,,,,"Arc beam numbers +" +scattered,b,a,,,,"Subtract scattered light?" +fitflat,b,a,,,,"Fit and ratio flat field spectrum?" +recenter,b,a,,,,"Recenter object apertures?" +edit,b,a,,,,"Edit/review object apertures?" +trace,b,a,,,,"Trace object spectra?" +arcap,b,a,,,,"Use object apertures for arcs?" +clean,b,a,,,,"Detect and replace bad pixels?" +dispcor,b,a,,,,"Dispersion correct spectra?" +savearcs,b,a,,,,"Save internal arcs?" +skyalign,b,a,,,,"Align sky lines?" +skysubtract,b,a,,,,"Subtract sky?" +skyedit,b,a,,,,"Edit the sky spectra?" +saveskys,b,a,,,,"Save sky spectra?" +splot,b,a,,,,"Plot the final spectrum?" +redo,b,a,,,,"Redo operations if previously done?" +update,b,a,,,,"Update spectra if cal data changes?" +batch,b,a,,,,"Extract objects in batch?" +listonly,b,a,,,,"List steps but don\'t process? +" +datamax,r,h,INDEF,,,"Max data value / cosmic ray threshold" +ansskyedit,s,q,"yes",,,"Edit the sky spectra?" +newaps,b,h,,,, +newresp,b,h,,,, +newdisp,b,h,,,, +newarcs,b,h,,,, +dobatch,b,h,,,, +anssplot,s,q,"yes",no|yes|NO|YES,,"Splot spectrum?" +extn,s,h,".ms",,,"Extraction extension" +fd1,*struct,h,"",,, +fd2,*struct,h,"",,, +fd3,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/skysub.cl b/noao/imred/src/fibers/skysub.cl new file mode 100644 index 00000000..b7522591 --- /dev/null +++ b/noao/imred/src/fibers/skysub.cl @@ -0,0 +1,145 @@ +# SKYSUB -- Sky subtract fiber spectra. +# Subtract the selected sky apertures from the selected object apertures. +# The object apertures may include the sky apertures if desired. +# The sky apertures are combined to a single master sky which is subtracted +# from each selected object aperture. match. The subtracted sky may be +# saved in an image with the prefix "sky" and the same output name. Note +# that existing output images are clobbered. + +procedure skysub (input) + +string input = "" {prompt="Input spectra to sky subtract"} + +string output = "" {prompt="Output sky subtracted spectra"} +string objaps = "" {prompt="Object apertures"} +string skyaps = "" {prompt="Sky apertures"} +string objbeams = "" {prompt="Object beam numbers"} +string skybeams = "" {prompt="Sky beam numbers"} +bool skyedit = yes {prompt="Edit the sky spectra?"} +string combine = "average" {prompt="Combining option", + enum="average|median|sum"} +string reject = "avsigclip" {prompt="Sky rejection option", + enum="none|minmax|avsigclip"} +string scale = "none" {prompt="Sky scaling option", + enum="none|mode|median|mean"} +bool saveskys = yes {prompt="Save sky spectra?"} +file logfile = "" {prompt="Logfile"} + +struct *fd1 +struct *fd2 +struct *fd3 + +begin + string imtype, mstype + string in, out, out1, sky, log, aps, str, str2 + file temp1, temp2, temp3, temp4 + int i, j, n + + imtype = "." // envget ("imtype") + i = stridx (",", imtype) + if (i > 0) + imtype = substr (imtype, 1, i-1) + n = strlen (imtype) + + temp1 = mktemp ("tmp$iraf") + temp2 = mktemp ("tmp$iraf") + temp3 = mktemp ("tmp$iraf") + temp4 = mktemp ("tmp$iraf") + + if (logfile == "") + log = "dev$null" + else + log = logfile + + sections (input, option="fullname", > temp1) + sections (output, option="fullname", > temp2) + fd1 = temp1 + fd2 = temp2 + while (fscan (fd1, in) != EOF) { + i = strlen (in) + if (i > n && substr (in, i-n+1, i) == imtype) + in = substr (in, 1, i-n) + if (fscan (fd2, out) < 1) + out = in + out1 = out + i = strlen (out1) + if (i > 3 && substr (out1, i-2, i) == ".ms") + out1 = substr (out1, 1, i-3) + + aps = skyaps + sky = "sky" // out1 + if (access (sky // imtype)) + imdelete (sky, verify=no) + if (skyedit) { + scopy (in, sky, w1=INDEF, w2=INDEF, apertures=aps, bands="1", + beams=skybeams, apmodulus=0, offset=0, clobber=yes, + format="multispec", merge=no, renumber=no, + verbose=yes, >> "dev$null") + specplot (sky, apertures="", bands="1", autolayout=no, + autoscale=yes, fraction=1., scale=1., offset=0., + step=0., ptype="1", labels="user", ulabels="", + sysid=yes, yscale=yes, xlpos=1.02, ylpos=0., + title="Edit sky spectra from "//in, xlabel="", + ylabel="", xmin=INDEF, xmax=INDEF, ymin=INDEF, + ymax=INDEF, logfile=temp4, graphics="stdgraph") + imdelete (sky, verify=no) + system.match (sky, temp4, stop=no) | + fields (fields="2", lines="1-9999") | + system.sort (column=0, ignore=yes, numeric=no, + reverse_sort=no) | + lists.unique (> temp3) + delete (temp4, verify=no) + aps = "@" // temp4 + fd3 = temp3 + while (fscan (fd3, str) != EOF) { + i = stridx ("(", str) + j = stridx (")", str) + if (i > 0 && j > i) + str = substr(str,i+1,j-1) + else + str = "" + print (str, >> temp4) + } + fd3 = ""; delete (temp3, verify=no) + + reject.p_mode="q" + str = reject + reject.p_mode="h" + } + + if (skybeams == "") { + scombine (in, sky, noutput="", logfile=logfile, + apertures=aps, group="all", combine=combine, + reject=reject, first=yes, scale=scale, zero="none", + weight="none", sample="", lthreshold=INDEF, + hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, + lsigma=3., hsigma=2., rdnoise="0.", gain="1.", snoise="0.", + sigscale=0., pclip=-0.5, grow=0, blank=0.) + } else { + temp3 = mktemp ("sky") + scopy (in, sky, w1=INDEF, w2=INDEF, apertures=aps, bands="", + beams=skybeams, apmodulus=0, offset=0, clobber=yes, + format="multispec", merge=no, renumber=no, + verbose=yes, >> log) + scombine (sky, temp3, noutput="", logfile=logfile, + apertures=aps, group="all", combine=combine, + reject=reject, first=yes, scale=scale, zero="none", + weight="none", sample="", lthreshold=INDEF, + hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, + lsigma=3., hsigma=2., rdnoise="0.", gain="1.", snoise="0.", + sigscale=0., pclip=-0.5, grow=0, blank=0.) + flpr + imdelete (sky, verify=no) + imrename (temp3, sky, verbose=yes, >> log) + } + sarith (in, "-", sky, out, w1=INDEF, w2=INDEF, apertures=objaps, + bands="", beams=objbeams, reverse=no, ignoreaps=yes, + format="multispec", renumber=no, offset=0, clobber=yes, + merge=no, errval=0., verbose=yes, >> log) + if (!saveskys) + imdelete (sky, verify=no) + } + fd1 = ""; delete (temp1, verify=no) + fd2 = ""; delete (temp2, verify=no) + delete (temp4, verify=no, >>& "dev$null") +end diff --git a/noao/imred/src/fibers/skysub.par b/noao/imred/src/fibers/skysub.par new file mode 100644 index 00000000..bf1f4c2c --- /dev/null +++ b/noao/imred/src/fibers/skysub.par @@ -0,0 +1,16 @@ +input,s,a,"",,,"Input spectra to sky subtract" +output,s,h,"",,,"Output sky subtracted spectra" +objaps,s,h,"",,,"Object apertures" +skyaps,s,h,"",,,"Sky apertures" +objbeams,s,h,"",,,"Object beam numbers" +skybeams,s,h,"",,,"Sky beam numbers" +skyedit,b,h,yes,,,"Edit the sky spectra?" +combine,s,h,"average",average|median|sum,,"Combining option" +reject,s,h,"avsigclip",none|minmax|avsigclip,,"Sky rejection option" +scale,s,h,"none",none|mode|median|mean,,"Sky scaling option" +saveskys,b,h,yes,,,"Save sky spectra?" +logfile,f,h,"",,,"Logfile" +fd1,*struct,h,"",,, +fd2,*struct,h,"",,, +fd3,*struct,h,"",,, +mode,s,h,"ql",,, diff --git a/noao/imred/src/fibers/temp b/noao/imred/src/fibers/temp new file mode 100644 index 00000000..2205597e --- /dev/null +++ b/noao/imred/src/fibers/temp @@ -0,0 +1,16 @@ +Prototype Data Manager Keywords + +COMMAND +STATUS +EOF + +LABEL + +IMAGEID +MJD-OBS + +FILTER +XTALK +BPM +ZERO +FLAT |