Initial commit

1 files changed, 437 insertions, 0 deletions

diff --git a/noao/onedspec/t_calibrate.x b/noao/onedspec/t_calibrate.x
new file mode 100644
index 00000000..5df62b45
--- /dev/null
+++ b/noao/onedspec/t_calibrate.x
@@ -0,0 +1,437 @@
+include	<error.h>
+include	<imset.h>
+include	<imhdr.h>
+include	<math/iminterp.h>
+include	<smw.h>
+
+define	EXTN_LOOKUP	10	# Interp index for de-extinction
+define	VLIGHT	2.997925e18	# Speed of light, Angstroms/sec
+
+# T_CALIBRATE -- Apply extinction correction and flux calibration to spectra.
+# The sensitivity function derived from the tasks STANDARD and SENSFUNC
+# are applied to the given spectra.  The output may be the same as the
+# input or new spectra may be created.
+#
+# The sensitivity function is contained in an image having its aperture
+# number indicated by the trailing integer of the image filename.
+# An option, "ignoreaps", can be set to override the appending of the
+# aperture number on those cases where no aperture correspondence is
+# appropriate.
+
+procedure t_calibrate ()
+
+pointer	inlist		# Input list
+pointer	outlist		# Output list
+pointer	sens		# Sensitivity image root name
+pointer	ob		# Observatory
+bool	ignoreaps	# Ignore aperture numbers?
+bool	extinct		# Apply extinction correction?
+bool	flux		# Apply flux calibration?
+bool	fnu		# Calibration flux in FNU?
+
+bool	doextinct, doflux, newobs, obshead
+int	i, j, k, l, n, enwaves, nout, ncal
+real	a, latitude, time, ext, fcor, ical, w, dw
+pointer	sp, input, output, temp
+pointer	obs, in, smw, sh, out, ewaves, emags, pcal, cal, asi, x, y, data
+
+int	imtgetim(), imtlen()
+bool	clgetb(), streq()
+real	clgetr(), obsgetr(), asieval()
+double	shdr_lw(), shdr_wl()
+pointer	imtopenp(), immap(), smw_openim(), imgl3r(), impl3r()
+errchk	immap, smw_openim, shdr_open, imgl3r, impl3r
+errchk	obsimopen, get_airm, ext_load, cal_getflux, cal_extn, cal_flux
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (sens, SZ_FNAME, TY_CHAR)
+	call salloc (ob, SZ_FNAME, TY_CHAR)
+	call salloc (temp, SZ_LINE, TY_CHAR)
+
+	# Get task parameters.
+	inlist = imtopenp ("input")
+	outlist = imtopenp ("output")
+	call clgstr ("records", Memc[temp], SZ_LINE)
+	call odr_openp (inlist, Memc[temp])
+	call odr_openp (outlist, Memc[temp])
+	call clgstr ("sensitivity", Memc[sens], SZ_FNAME)
+	call clgstr ("observatory", Memc[ob], SZ_FNAME)
+	extinct = clgetb ("extinct")
+	flux = clgetb ("flux")
+	fnu = clgetb ("fnu")
+	ignoreaps = clgetb ("ignoreaps")
+
+	if (!extinct && !flux)
+	    call error (0, "No calibration correction specified")
+
+	# Loop over all input images.
+	sh = NULL
+	obs = NULL
+	enwaves = 0
+	ncal = 0
+	while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+
+	    # Set output image.  Use a temporary image when output=input.
+	    if (imtlen (outlist) > 0) {
+		if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF)
+		    break
+	    } else
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+
+	    # Map the input image.
+	    iferr (in = immap (Memc[input], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+	    smw = smw_openim (in)
+
+	    # Check the input image calibration status.
+	    call shdr_open (in, smw, 1, 1, INDEFI, SHHDR, sh)
+	    if (DC(sh) == DCNO) {
+		call eprintf ("WARNING: [%s] has no dispersion function\n")
+		    call pargstr (Memc[input])
+		call smw_close (MW(sh))
+		call imunmap (in)
+		next
+	    }
+	    call shdr_units (sh, "angstroms")
+
+	    doextinct = extinct && (EC(sh) == ECNO)
+	    doflux = flux && (FC(sh) == FCNO)
+	    if (!(doextinct || doflux)) {
+		call eprintf ("WARNING: [%s] is already calibrated\n")
+		    call pargstr (Memc[input])
+		call smw_close  (MW(sh))
+		call imunmap (in)
+		next
+	    }
+
+	    # Map the output image.
+	    if (streq (Memc[input], Memc[output]))
+		call mktemp ("temp", Memc[temp], SZ_LINE)
+	    else
+		call strcpy (Memc[output], Memc[temp], SZ_LINE)
+	    out = immap (Memc[temp], NEW_COPY, in)
+	    if (IM_PIXTYPE(out) != TY_DOUBLE)
+		IM_PIXTYPE(out) = TY_REAL
+
+	    # Log the operation.
+	    call printf ("%s: %s\n")
+		call pargstr (Memc[output])
+		call pargstr (IM_TITLE(out))
+	    call flush (STDOUT)
+
+	    # Initialize the extinction correction.
+	    if (doextinct) {
+		EC(sh) = ECYES
+
+		# Load extinction function.
+		if (enwaves == 0) {
+		    call ext_load (ewaves, emags, enwaves)
+		    call intrp0 (EXTN_LOOKUP)
+		}
+
+		# Determine airmass if needed.
+		if (IS_INDEF(AM(sh))) {
+		    call obsimopen (obs, in, Memc[ob], NO, newobs, obshead)
+		    if (newobs)
+			call obslog (obs, "CALIBRATE", "latitude", STDOUT)
+		    latitude = obsgetr (obs, "latitude")
+		    iferr (call get_airm (RA(sh), DEC(sh), HA(sh), ST(sh),
+			latitude, AM(sh))) {
+			call printf ("%s: ")
+			    call pargstr (Memc[input])
+			call flush (STDOUT)
+			AM(sh) = clgetr ("airmass")
+			call imunmap (in)
+			ifnoerr (in = immap (Memc[input], READ_WRITE, 0)) {
+			    IM(sh) = in
+			    call imseti (IM(sh), IM_WHEADER, YES)
+			    call imaddr (IM(sh), "airmass", AM(sh))
+			} else {
+			    in = immap (Memc[input], READ_ONLY, 0)
+			    IM(sh) = in
+			}
+		    }
+		}
+		a = AM(sh)
+	    } else
+		ext = 1.
+
+	    # Initialize the flux correction.
+	    nout = 0
+	    if (doflux) {
+		FC(sh) = FCYES
+
+		if (IS_INDEF (IT(sh)) || IT(sh) <= 0.) {
+		    call printf ("%s: ")
+			call pargstr (Memc[input])
+		    call flush (STDOUT)
+		    IT(sh) = clgetr ("exptime")
+		    call imunmap (in)
+		    ifnoerr (in = immap (Memc[input], READ_WRITE, 0)) {
+			IM(sh) = in
+			call imseti (IM(sh), IM_WHEADER, YES)
+			call imaddr (IM(sh), "exptime", IT(sh))
+			call imaddr (out, "exptime", IT(sh))
+		    } else {
+		        in = immap (Memc[input], READ_ONLY, 0)
+			IM(sh) = in
+		    }
+		}
+		time = IT(sh)
+	    } else
+		fcor = 1.
+
+	    # Calibrate.
+	    do j = 1, IM_LEN(in,3) {
+		do i = 1, IM_LEN(in,2) {
+		    data = impl3r (out, i, j)
+		    switch (SMW_FORMAT(smw)) {
+		    case SMW_ND:
+			if (doflux) {
+			    call cal_getflux (Memc[sens], INDEFI, fnu,
+				pcal, ncal, cal)
+
+			    asi = IM(cal)
+			    n = SN(cal)
+			}
+			y = imgl3r (in, i, j)
+			switch (SMW_LAXIS(smw,1)) {
+			case 1:
+			    do k = 1, IM_LEN(out,1) {
+				w = shdr_lw (sh, double(k))
+				if (doextinct) {
+				    call intrp (EXTN_LOOKUP, Memr[ewaves],
+					Memr[emags], enwaves, w, ext, l)
+				    ext = 10.0 ** (0.4 * a * ext)
+				}
+				if (doflux) {
+				    ical = shdr_wl (cal, double(w))
+				    if (ical < 1. || ical > n) {
+					if (ical < 0.5 || ical > n + 0.5)
+					    nout = nout + 1
+					ical = max (1., min (real(n), ical))
+				    }
+				    dw = abs (shdr_lw (sh, double(k+0.5)) -
+					shdr_lw (sh, double(k-0.5)))
+				    fcor = asieval (asi, ical) / dw / time
+				}
+				Memr[data] = Memr[y] * ext * fcor
+				y = y + 1
+				data = data + 1
+			    }
+			case 2, 3:
+			    if (SMW_LAXIS(smw,1) == 2)
+				k = i
+			    else
+				k = j
+			    w = shdr_lw (sh, double(k))
+			    if (doextinct) {
+				call intrp (EXTN_LOOKUP, Memr[ewaves],
+				    Memr[emags], enwaves, w, ext, l)
+				ext = 10.0 ** (0.4 * a * ext)
+			    }
+			    if (doflux) {
+				ical = shdr_wl (cal, double(w))
+				if (ical < 1. || ical > n) {
+				    if (ical < 0.5 || ical > n + 0.5)
+					nout = nout + 1
+				    ical = max (1., min (real(n), ical))
+				}
+				dw = abs (shdr_lw (sh, double(k+0.5)) -
+				    shdr_lw (sh, double(k-0.5)))
+				fcor = asieval (asi, ical) / dw / time
+			    }
+			    call amulkr (Memr[y], ext * fcor, Memr[data],
+				IM_LEN(out,1))
+			}
+		    case SMW_ES, SMW_MS:
+			call shdr_open (in, smw, i, j, INDEFI, SHDATA, sh)
+			call shdr_units (sh, "angstroms")
+			if (doflux) {
+			    if (ignoreaps)
+				call cal_getflux (Memc[sens], INDEFI, fnu,
+				    pcal, ncal, cal)
+			    else
+				call cal_getflux (Memc[sens], AP(sh), fnu,
+				    pcal, ncal, cal)
+
+			    asi = IM(cal)
+			    n = SN(cal)
+			}
+			x = SX(sh)
+			y = SY(sh)
+			do k = 1, SN(sh) {
+			    w = Memr[x]
+			    if (doextinct) {
+				call intrp (EXTN_LOOKUP, Memr[ewaves],
+				    Memr[emags], enwaves, w, ext, l)
+				ext = 10.0 ** (0.4 * a * ext)
+			    }
+			    if (doflux) {
+				ical = shdr_wl (cal, double(w))
+				if (ical < 1. || ical > n) {
+				    if (ical < 0.5 || ical > n + 0.5)
+					nout = nout + 1
+				    ical = max (1., min (real(n), ical))
+				}
+				dw = abs (shdr_lw (sh, double(k+0.5)) -
+				    shdr_lw (sh, double(k-0.5)))
+				fcor = asieval (asi, ical) / dw / time
+			    }
+			    Memr[data] = Memr[y] * ext * fcor
+			    x = x + 1
+			    y = y + 1
+			    data = data + 1
+			}
+			do k = SN(sh)+1, IM_LEN(out,1) {
+			    Memr[data] = 0
+			    data = data + 1
+			}
+		    }
+		}
+	    }
+
+	    # Log the results.
+	    if (doflux && (IS_INDEF (IT(sh)) || IT(sh) <= 0.)) {
+		call printf (
+		    "  WARNING: No exposure time found.  Using a time of %g.\n")
+		    call pargr (time)
+	    }
+	    if (nout > 0) {
+		call printf (
+		    "  WARNING: %d pixels outside of flux calibration limits\n")
+		    call pargi (nout)
+	    }
+	    if (doextinct)
+		call printf ("  Extinction correction applied\n")
+	    if (doflux)
+		call printf ("  Flux calibration applied\n")
+	    call flush (STDOUT)
+
+	    call imaddr (out, "AIRMASS", AM(sh))
+	    call imaddi (out, "EX-FLAG", EC(sh))
+	    call imaddi (out, "CA-FLAG", FC(sh))
+	    if (doflux) {
+		if (fnu)
+		    call imastr (out, "BUNIT", "erg/cm2/s/Hz")
+		else
+		    call imastr (out, "BUNIT", "erg/cm2/s/A")
+	    }
+
+	    # Close the input and output images.
+	    call smw_close (MW(sh))
+	    call imunmap (in)
+	    call imunmap (out)
+	    if (streq (Memc[input], Memc[output])) {
+		call imdelete (Memc[input])
+		call imrename (Memc[temp], Memc[output])
+	    }
+	}
+
+	# Finish up.
+	if (enwaves > 0) {
+	    call mfree (ewaves, TY_REAL)
+	    call mfree (emags, TY_REAL)
+	}
+	if (ncal > 0) {
+	    do i = 0, ncal-1 {
+		cal = Memi[pcal+i]
+		call asifree (IM(cal))
+		call smw_close (MW(cal))
+		call shdr_close (cal)
+	    }
+	    call mfree (pcal, TY_POINTER)
+	}
+	if (obs != NULL)
+	    call obsclose (obs)
+	call shdr_close (sh)
+	call imtclose (inlist)
+	call imtclose (outlist)
+	call sfree (sp)
+end
+
+
+# CAL_GETFLUX -- Get flux calibration data
+# The sensitivity spectrum is in peculiar magnitudish units of 2.5 log10
+# [counts/sec/A / ergs/cm2/s/A].  This is converted back to reasonable
+# numbers to be multiplied into the data spectra.  An interpolation function
+# is then fit and stored in the image pointer field.  For efficiency the
+# calibration data is saved by aperture so that additional calls simply
+# return the data pointer.
+
+procedure cal_getflux (sens, ap, fnu, pcal, ncal, cal)
+
+char	sens[ARB]		# Sensitivity function image or rootname
+int	ap			# Aperture
+bool	fnu			# Fnu units?
+pointer	pcal			# Pointer to cal data
+int	ncal			# Number of active cal data structures
+pointer	cal			# Calibration data structure
+
+int	i, j, n, clgwrd()
+pointer	sp, fname, im, smw, x, y, immap(), smw_openim()
+errchk	immap, smw_openim, shdr_open, asifit
+
+begin
+	# Check for previously saved calibration
+	for (i=0; i<ncal; i=i+1) {
+	    cal = Memi[pcal+i]
+	    if (AP(cal) == ap)
+		return
+	}
+
+	# Allocate space for a new data pointer, get the calibration data,
+	# and convert to calibration array.
+
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	if (ncal == 0)
+	    call malloc (pcal, 10, TY_POINTER)
+	else if (mod (ncal, 10) == 0)
+	    call realloc (pcal, ncal+10, TY_POINTER)
+
+	if (IS_INDEFI(ap))
+	    call strcpy (sens, Memc[fname], SZ_FNAME)
+	else {
+	    call sprintf (Memc[fname], SZ_FNAME, "%s.%04d")
+		call pargstr (sens)
+		call pargi (ap)
+	}
+
+	im = immap (Memc[fname], READ_ONLY, 0)
+	smw = smw_openim (im)
+	cal = NULL
+	call shdr_open (im, smw, 1, 1, ap, SHDATA, cal)
+	call shdr_units (cal, "angstroms")
+	AP(cal) = ap
+	Memi[pcal+ncal] = cal
+	ncal = ncal + 1
+	call imunmap (im)
+
+	x = SX(cal)
+	y = SY(cal)
+	n = SN(cal)
+	do j = 1, n {
+	    Memr[y] = 10.0 ** (-0.4 * Memr[y])
+	    if (fnu) {
+		Memr[y] = Memr[y] * Memr[x] ** 2 / VLIGHT
+		x = x + 1
+	    }
+	    y = y + 1
+	}
+
+	call asiinit (im, clgwrd ("interp", Memc[fname], SZ_FNAME,II_FUNCTIONS))
+	call asifit (im, Memr[SY(cal)], n)
+	IM(cal) = im
+
+	call mfree (SX(cal), TY_REAL)
+	call mfree (SY(cal), TY_REAL)
+
+	call sfree (sp)
+end