Initial commit

1 files changed, 234 insertions, 0 deletions

diff --git a/noao/onedspec/sensfunc/sfimage.x b/noao/onedspec/sensfunc/sfimage.x
new file mode 100644
index 00000000..71edc213
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfimage.x
@@ -0,0 +1,234 @@
+include	<gset.h>
+include	<math/curfit.h>
+include	"sensfunc.h"
+include	<smw.h>
+
+
+# SF_IMAGE -- Graph fluxed image data and possible standard flux points.
+# For efficiency the IMIO pointer, buffer, and associated data are kept
+# since a redraw is a common occurence and generating the data is slow.
+
+procedure sf_image (gp, wc, stds, nstds, cv, wextn, extn, nextn, ecv)
+
+pointer	gp			# Graphics structure
+int	wc			# WC of graph
+pointer	stds[nstds]		# Standard star data for flux points
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction table values
+pointer	ecv			# Residual extinction curve
+
+int	scale[SF_NGRAPHS], log[SF_NGRAPHS]
+
+bool	newobs, obshead
+int	i, j, n, err
+real	a, t, w, dw, e, sens, latitude, smin, smax, xmin, xmax
+pointer	im, mw, sh, skyim, skymw, skysh, std, gio, sp, str, x, y, z, obs
+pointer	immap(), smw_openim()
+real	cveval(), obsgetr(), cvstatr()
+double	shdr_lw()
+bool	streq(), strne()
+errchk	immap, smw_openim, obsimopen
+
+define	plot_	99
+
+begin
+	# Return if no image name.
+	if (Memc[GP_IMAGES(gp,wc)] == EOS)
+	    return
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the spectrum and sky subtract if necessary.
+	sh = GP_SHDR(gp,wc)
+	if (sh != NULL) {
+	    if (streq (Memc[GP_IMAGES(gp,wc)], IMNAME(sh))) {
+		if (GP_LOG(gp) == log[wc])
+		    goto plot_
+		else
+		    call shdr_close (sh)
+	    }
+	}
+
+	# Determine a valid standard star to get aperture number.
+	do i = 1, nstds
+	    if (STD_FLAG(stds[i]) != SF_EXCLUDE) {
+		std = stds[i]
+		break
+	    }
+
+	im = immap (Memc[GP_IMAGES(gp,wc)], READ_ONLY, 0)
+	mw = smw_openim (im)
+	call shdr_open (im, mw, 1, 1, STD_BEAM(std), SHDATA, sh)
+
+	# Check for dispersion correction
+	if (DC(sh) == DCNO) {
+	    call shdr_close (sh)
+	    call smw_close (mw)
+	    call imunmap (im)
+	    GP_SHDR(gp,wc) = NULL
+	    call sfree (sp)
+	    call printf ("-%s must be dispersion corrected-")
+		call pargstr (Memc[GP_IMAGES(gp,wc)])
+	    return
+	}
+
+	# Sky subtract if necessary
+	if (Memc[GP_SKYS(gp,wc)] != EOS) {
+	    skyim = immap (Memc[GP_SKYS(gp,wc)], READ_ONLY, 0)
+	    skymw = smw_openim (skyim)
+	    call shdr_open (skyim, skymw, 1, 1, STD_BEAM(std), SHDATA, skysh)
+	    call shdr_rebin (skysh, sh)
+	    call asubr (Memr[SY(sh)], Memr[SY(skysh)], Memr[SY(sh)], SN(sh))
+	    call shdr_close (skysh)
+	    call smw_close (skymw)
+	    call imunmap (skyim)
+	}
+
+	# Set airmass and exposure time
+	if (IS_INDEF (AM(sh))) {
+	    obs = NULL
+	    call clgstr ("observatory", Memc[str], SZ_LINE)
+	    call obsimopen (obs, im, Memc[str], NO, newobs, obshead)
+	    latitude = obsgetr (obs, "latitude")
+	    call obsclose (obs)
+	    call get_airm (RA(sh), DEC(sh), HA(sh), ST(sh), latitude,
+		AM(sh))
+	}
+	a = AM(sh)
+	if (IS_INDEF (IT(sh)))
+	    t = 1.
+	else
+	    t = IT(sh)
+
+	# Apply extinction correction if needed
+	if (EC(sh) == ECNO) {
+	    if (ecv != NULL) {
+	        xmin = cvstatr (ecv, CVXMIN)
+	        xmax = cvstatr (ecv, CVXMAX)
+	    }
+	    do i = 1, SN(sh) {
+		w = Memr[SX(sh)+i-1]
+	        call intrp (1, wextn, extn, nextn, w, e, err)
+	        if (ecv != NULL)
+		    e = e + cveval (ecv, min (xmax, max (w, xmin)))
+	        Memr[SY(sh)+i-1] = Memr[SY(sh)+i-1] * 10. ** (0.4 * a * e)
+	    }
+	} else {
+	    call printf ("-%s already extinction corrected-")
+	       call pargstr (Memc[GP_IMAGES(gp,wc)])
+	}
+
+	# Apply flux calibration if needed
+	if (FC(sh) == FCNO) {
+	    do i = 1, SN(sh) {
+	        w = Memr[SX(sh)+i-1]
+		dw = abs (shdr_lw (sh, double (i+0.5)) -
+		    shdr_lw (sh, double (i-0.5)))
+	        sens = cveval (cv, w)
+	        Memr[SY(sh)+i-1] = Memr[SY(sh)+i-1] / t / dw / 10.**(0.4*sens)
+	    }
+	} else {
+	    call printf ("-%s already flux calibrated-")
+		call pargstr (Memc[GP_IMAGES(gp,wc)])
+	}
+
+	# Set flux scaling
+	call alimr (Memr[SY(sh)], SN(sh), smin, smax)
+	if (smax < 0.)
+	    scale[wc] = 0.
+	else if (GP_LOG(gp) == NO) {
+	    scale[wc] = -log10 (smax) + 1
+	    w = 10. ** scale[wc]
+	    call amulkr (Memr[SY(sh)], w, Memr[SY(sh)], SN(sh))
+	} else {
+	    scale[wc] = INDEFI
+	    smin = smax / 1000.
+	    w = smax
+	    y = SY(sh)
+	    do i = 1, SN(sh) {
+		if (Memr[y] > smin)
+		    w = min (w, Memr[y])
+		y = y + 1
+	    }
+	    y = SY(sh)
+	    do i = 1, SN(sh) {
+		Memr[y] = log10 (max (Memr[y], w))
+		y = y + 1
+	    }
+	}
+	log[wc] = GP_LOG(gp)
+
+	# Save the spectrum for future redraw.
+	call smw_close (MW(sh))
+	call imunmap (im)
+	GP_SHDR(gp,wc) = sh
+
+plot_
+	# Plot scaled graph.
+	smin = GP_FMIN(gp)
+	smax = GP_FMAX(gp)
+	if (IS_INDEFI(scale[wc])) {
+	    call sprintf (Memc[str], SZ_LINE, "%s: Log Flux")
+	        call pargstr (Memc[GP_IMAGES(gp,wc)])
+	    if (!IS_INDEF(smin)) {
+		if (smin > 0.)
+		    smin = log10 (smin)
+		else
+		    smin = INDEF
+	    }
+	    if (!IS_INDEF(smax)) {
+		if (smax > 0.)
+		    smax = log10 (smax)
+		else
+		    smax = INDEF
+	    }
+	} else if (scale[wc] != 0) {
+	    call sprintf (Memc[str], SZ_LINE, "%s: Flux x 1E%d")
+	        call pargstr (Memc[GP_IMAGES(gp,wc)])
+	        call pargi (scale[wc])
+	    w = 10. ** scale[wc]
+	    if (!IS_INDEF(smin))
+		smin = w * smin
+	    if (!IS_INDEF(smax))
+		smax = w * smax
+	} else {
+	    call sprintf (Memc[str], SZ_LINE, "%s: Flux")
+	        call pargstr (Memc[GP_IMAGES(gp,wc)])
+	    w = 1.
+	}
+
+	gio = GP_GIO(gp)
+	call gascale (gio, Memr[SX(sh)], SN(sh), 1)
+	call gascale (gio, Memr[SY(sh)], SN(sh), 2)
+	call gswind (gio, INDEF, INDEF, smin, smax)
+	call glabax (gio, Memc[str], "", "")
+	call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
+	call gpline (gio, Memr[SX(sh)], Memr[SY(sh)], SN(sh))
+
+	call sfree (sp)
+
+	# Check if image is one of the standard stars and plot flux points.
+	do i = 1, nstds {
+	    if (strne (Memc[GP_IMAGES(gp,wc)], STD_IMAGE(stds[i])))
+		next
+	    n = STD_NWAVES(stds[i])
+	    x = STD_WAVES(stds[i])
+	    y = STD_FLUXES(stds[i])
+	    z = STD_DWAVES(stds[i])
+	    call gseti (gio, G_PMLTYPE, 1)
+	    call gseti (gio, G_PLCOLOR, GP_CMARK(gp))
+	    if (IS_INDEFI(scale[wc])) {
+	        do j = 0, n-1
+		    call gmark (gio, Memr[x+j], log10 (Memr[y+j]), GM_HEBAR,
+			-Memr[z+j], 1.)
+	    } else {
+	        do j = 0, n-1
+		    call gmark (gio, Memr[x+j], w * Memr[y+j], GM_HEBAR,
+			-Memr[z+j], 1.)
+	    }
+	}
+end