diff options
Diffstat (limited to 'noao/onedspec/sensfunc/sfimage.x')
| -rw-r--r-- | noao/onedspec/sensfunc/sfimage.x | 234 |
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 |