diff options
Diffstat (limited to 'noao/onedspec/sensfunc')
32 files changed, 3675 insertions, 0 deletions
diff --git a/noao/onedspec/sensfunc/mkpkg b/noao/onedspec/sensfunc/mkpkg new file mode 100644 index 00000000..4fdd11c0 --- /dev/null +++ b/noao/onedspec/sensfunc/mkpkg @@ -0,0 +1,38 @@ +# SENSFUNC + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + sfadd.x sensfunc.h <gset.h> + sfapertures.x sensfunc.h + sfcgraph.x sensfunc.h <gset.h> + sfcolon.x sensfunc.h <error.h> <gset.h> + sfcolors.x sensfunc.h <gset.h> + sfcomposite.x sensfunc.h + sfdata.x sensfunc.h + sfdelete.x sensfunc.h <gset.h> + sfeout.x <ctype.h> <error.h> <mach.h> <math/curfit.h> + sfextinct.x sensfunc.h <pkg/gtools.h> + sffit.x sensfunc.h <math/curfit.h> + sfginit.x sensfunc.h <gset.h> + sfgraph.x sensfunc.h <error.h> <gset.h> <math/curfit.h> + sfimage.x sensfunc.h <smw.h> <gset.h> <math/curfit.h> + sfmarks.x sensfunc.h <gset.h> + sfmove.x sensfunc.h <gset.h> + sfnearest.x sensfunc.h <gset.h> <mach.h> + sfoutput.x sensfunc.h <imhdr.h> <mach.h> + sfreset.x sensfunc.h + sfrms.x sensfunc.h + sfsensfunc.x sensfunc.h <error.h> <gset.h> <mach.h> + sfshift.x sensfunc.h + sfstats.x sensfunc.h + sfstds.x sensfunc.h + sftitle.x sensfunc.h + sfundelete.x sensfunc.h <gset.h> + sfvstats.x sensfunc.h + sfweights.x sensfunc.h + t_sensfunc.x sensfunc.h + ; diff --git a/noao/onedspec/sensfunc/sensfunc.h b/noao/onedspec/sensfunc/sensfunc.h new file mode 100644 index 00000000..9e3afb77 --- /dev/null +++ b/noao/onedspec/sensfunc/sensfunc.h @@ -0,0 +1,64 @@ +# SENSFUNC definitions. + +define SF_NGRAPHS 4 # Number of graphs per frame +define SF_INCLUDE 1 # Include observation +define SF_EXCLUDE 2 # Exclude observation +define SF_DELETE 3 # Delete observation + +# SENSFUNC Standard Star Data Structure. + +define SZ_STDIMAGE 63 # Length of standard image name +define SZ_STDTITLE 63 # Length of standard title + +define LEN_STD 115 # Length of standard obs. structure + +define STD_IMAGE Memc[P2C($1)] # Standard image name +define STD_SKY Memc[P2C($1+32)] # Standard image sky +define STD_TITLE Memc[P2C($1+64)] # Standard title +define STD_FLAG Memi[$1+96] # Flag +define STD_BEAM Memi[$1+97] # Beam number of spectrum +define STD_NPTS Memi[$1+98] # Number of points in spectrum +define STD_EXPTIME Memr[P2R($1+99)] # Exposure time +define STD_AIRMASS Memr[P2R($1+100)] # Airmass of spectrum +define STD_WSTART Memr[P2R($1+101)] # Starting wavelength of spectrum +define STD_WEND Memr[P2R($1+102)] # Ending wavelength of spectrum +define STD_SHIFT Memr[P2R($1+103)] # Added shift +define STD_NWAVES Memi[$1+104] # Number of calibration wavelengths +define STD_WAVES Memi[$1+105] # Pointer to wavelengths +define STD_FLUXES Memi[$1+106] # Pointer to standard flux values +define STD_DWAVES Memi[$1+107] # Pointer to flux bandwidths +define STD_COUNTS Memi[$1+108] # Pointer to counts +define STD_SENS Memi[$1+109] # Pointer to sensitivities +define STD_FIT Memi[$1+110] # Pointer to fitted sensitivities +define STD_WTS Memi[$1+111] # Pointer to weights +define STD_IWTS Memi[$1+112] # Pointer to weights +define STD_X Memi[$1+114] # Pointer to plotted x values +define STD_Y Memi[$1+115] # Pointer to plotted y values + +# Graphics structure + +define GP_SZTITLE 79 # Size of title string + +define LEN_GP 75 # Length of structure + +define GP_GIO Memi[$1] # GIO pointer +define GP_TITLE Memc[P2C($1+1)] # General title +define GP_GRAPHS Memc[P2C($1+41)+$2-1] # Graphs +define GP_IMAGES Memi[$1+44+$2-1] # Pointer to image names +define GP_SKYS Memi[$1+48+$2-1] # Pointer to sky names +define GP_MARK Memi[$1+52] # Mark type +define GP_SZMARK Memr[P2R($1+53)] # Mark size +define GP_CMARK Memi[$1+54] # Mark color +define GP_MDEL Memi[$1+55] # Deleted mark +define GP_SZMDEL Memr[P2R($1+56)] # Size of deleted mark +define GP_CDEL Memi[$1+57] # Color of deleted mark +define GP_MADD Memi[$1+58] # Mark type +define GP_CADD Memi[$1+59] # Mark color +define GP_PLCOLOR Memi[$1+60] # Line color +define GP_WSTART Memr[P2R($1+61)] # Starting wavelength for plots +define GP_WEND Memr[P2R($1+62)] # Ending wavelength for plots +define GP_LOG Memi[$1+63] # Log flux plots? +define GP_FMIN Memr[P2R($1+64)] # Minimum flux plot limit +define GP_FMAX Memr[P2R($1+65)] # Maximum flux plot limit +define GP_SHDR Memi[$1+65+$2] # SHDR pointer +define GP_AIRMASS Memr[P2R($1+69+$2)] # Airmass range of plots diff --git a/noao/onedspec/sensfunc/sensfunc.key b/noao/onedspec/sensfunc/sensfunc.key new file mode 100644 index 00000000..5f09739f --- /dev/null +++ b/noao/onedspec/sensfunc/sensfunc.key @@ -0,0 +1,81 @@ + SENSFUNC: Determine Sensitivity Function + +SUMMARY: + +? Help a Add data c Composite data d Delete data +e Extinction f Fit (overplot) g Fit (redraw) i Info +m Move data o Original data q Quit r Redraw +s Shift data u Undelete data w Change weights I Interrupt + +:function [type] :graphs [types] :images [images] :marks [types] +:order [value] :skys [images] :stats [file] :vstats [file] +:colors [colors] + +Graph types: a=(resid,airmass), c=(composite,lambda), e=(extinction,lambda) + i=(Fluxed image,lambda), l=(Log of fluxed image, lambda), + r=(resid, lambda), s=(Sensitivity,lambda) + + +CURSOR KEYS: + +? Print help +a Add a point at the cursor position +c Toggle composite points +d Delete point, star, or wavelength nearest the cursor +e Toggle residual extinction correction +f Fit data with a sensitivity function and overplot the fit +g Fit data with a sensitivity function and redraw the graphs +i Print information about point nearest the cursor +m Move point, star, wavelength nearest the cursor to new sensitivity +o Reset to original data +q Quit and write sensitivity function for current aperture +r Redraw graph(s) +s Toggle shift of standard stars to eliminate mean deviations +u Undelete point, star, or wavelength nearest the cursor +w Change weights of point, star, or wavelength nearest the cursor +I Interrupt task immediately + + +COLON COMMANDS AND ARGUMENTS: + +:flux [min] [max] Limits for flux calibrated graphs (INDEF for autoscale) +:function [type] Function to be fit to sensitivity data. The types are: + chebyshev - Chebyshev polynomial + legendre - Legendre polynomial + spline1 - Linear spline + spline3 - Cubic spline +:graphs [types] Graphs to be displayed (up to four). The types are: + a - Residual sensitivity vs airmass + c - Composite residuals and error bars vs wavelength + e - Extinction (and revised extinction) vs wavelength + i - Flux calibrated image vs wavelength + l - Log of flux calibrated image vs wavelength + r - Residual sensitivity vs wavelength + s - Sensitivity vs wavelength +:images [images] Images to flux calibrate and plot (up to four images) +:colors [colors] Line and mark colors to use for line and included, deleted, + and added points. The colors are specified as four + whitespace separated integers between 1 and 9. +:marks [marks] Mark types to use for included, deleted, and added points: + point, box, plus, cross, diamond, hline, vline, + hebar, vebar, circle +:order [order] Order of function (polynomial terms or spline pieces) +:skys [images] Sky images for flux calibration (optional, up to four images) +:stats [file] Statistics about standard stars and sensitivity fit +:vstats [file] Verbose statistics about standard stars and sensitivity fit + + +EXAMPLES: + +:func spline3 Select cubic spline function +:g srae Graph sensitivity, residuals, airmass, and extinction +:g sii Graph sensitivity and two images +:i n1.0004 n1.0008 Set first two images to graph (the defaults are + taken from the standard star list) +:skys n1.0005 Subtract this sky image from first image for calibration +:colors 2 Change only the line color to 2 +:colors 2 5 4 3 Change the line and mark colors +:m plus Change the mark type for included points and don't + change the deleted or added point mark type +:stats Print statistics to terminal +:vstats stdstats Print verbose statistics to file diff --git a/noao/onedspec/sensfunc/sfadd.x b/noao/onedspec/sensfunc/sfadd.x new file mode 100644 index 00000000..c77694d7 --- /dev/null +++ b/noao/onedspec/sensfunc/sfadd.x @@ -0,0 +1,105 @@ +include <gset.h> +include "sensfunc.h" + +# SF_ADD -- Add a point to the added point observation structure. +# The added star is the next to last of the standard stars. + +procedure sf_add (gp, stds, nstds, cv, wx, wy, wc) + +pointer gp # Graphics structure +pointer stds[nstds] # Standard star structures +int nstds # Number of standard stars +pointer cv # Sensitivity function curve +real wx # Cursor X value +real wy # Cursor Y value +int wc # Cursor WCS + +int nwaves +real wave, sen, fit, cveval() +pointer std, waves, sens, fits, wts, iwts, x, y +errchk malloc, realloc + +begin + # Convert from particular WCS to wavelength and sensitivity. In + # order to add a point the graph must be either sensitivity or + # residual verses wavelength. If not then return without adding + # a point. + + switch (GP_GRAPHS(gp,wc)) { + case 's': + wave = wx + fit = cveval (cv, wx) + sen = wy + case 'r': + wave = wx + fit = cveval (cv, wx) + sen = wy + fit + default: + return + } + + # Add the point to the next to last standard star. Allocate + # or reallocate memory as needed. Turn the added star on by + # setting INCLUDE flag. + + std = stds[nstds-1] + nwaves = STD_NWAVES(std) + 1 + waves = STD_WAVES(std) + if (waves == NULL) { + call malloc (waves, nwaves, TY_REAL) + call malloc (sens, nwaves, TY_REAL) + call malloc (fits, nwaves, TY_REAL) + call malloc (wts, nwaves, TY_REAL) + call malloc (iwts, nwaves, TY_REAL) + call malloc (x, nwaves, TY_REAL) + call malloc (y, nwaves, TY_REAL) + } else { + sens = STD_SENS(std) + fits = STD_FIT(std) + wts = STD_WTS(std) + iwts = STD_IWTS(std) + x = STD_X(std) + y = STD_Y(std) + call realloc (waves, nwaves, TY_REAL) + call realloc (sens, nwaves, TY_REAL) + call realloc (fits, nwaves, TY_REAL) + call realloc (wts, nwaves, TY_REAL) + call realloc (iwts, nwaves, TY_REAL) + call realloc (x, nwaves, TY_REAL) + call realloc (y, nwaves, TY_REAL) + } + STD_FLAG(std) = SF_INCLUDE + STD_NWAVES(std) = nwaves + STD_WAVES(std) = waves + STD_SENS(std) = sens + STD_FIT(std) = fits + STD_WTS(std) = wts + STD_IWTS(std) = iwts + STD_X(std) = x + STD_Y(std) = y + + Memr[waves+nwaves-1] = wave + Memr[sens+nwaves-1] = sen + Memr[fits+nwaves-1] = fit + Memr[wts+nwaves-1] = 1 + Memr[iwts+nwaves-1] = 1 + + # Mark the added point on all graphs. + for (wc = 1; GP_GRAPHS(gp,wc) != EOS; wc=wc+1) { + call gseti (GP_GIO(gp), G_WCS, wc) + call gseti (GP_GIO(gp), G_PLCOLOR, GP_CADD(gp)) + switch (GP_GRAPHS(gp,wc)) { + case 's': + call gmark (GP_GIO(gp), wave, sen, GP_MADD(gp), GP_SZMARK(gp), + GP_SZMARK(gp)) + case 'r': + wy = sen - cveval (cv, wave) + call gmark (GP_GIO(gp), wave, wy, GP_MADD(gp), GP_SZMARK(gp), + GP_SZMARK(gp)) + case 'a': + wy = sen - cveval (cv, wave) + call gmark (GP_GIO(gp), STD_AIRMASS(std), wy, GP_MADD(gp), + GP_SZMARK(gp), GP_SZMARK(gp)) + } + } +end diff --git a/noao/onedspec/sensfunc/sfapertures.x b/noao/onedspec/sensfunc/sfapertures.x new file mode 100644 index 00000000..7eb2b6f8 --- /dev/null +++ b/noao/onedspec/sensfunc/sfapertures.x @@ -0,0 +1,27 @@ +include "sensfunc.h" + +# SF_APERTURES -- Determine the apertures in use. + +procedure sf_apertures (stds, nstds, apertures, napertures) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +pointer apertures # Pointer to apertures (returned) +int napertures # Number of apertures (returned) + +int i, j, aperture + +errchk malloc, realloc + +begin + call malloc (apertures, nstds, TY_INT) + napertures = 0 + do i = 1, nstds { + aperture = STD_BEAM(stds[i]) + for (j=1; (j<=napertures)&&(aperture!=Memi[apertures+j-1]); j=j+1) + ; + napertures = max (napertures, j) + Memi[apertures+j-1] = aperture + } + call realloc (apertures, napertures, TY_INT) +end diff --git a/noao/onedspec/sensfunc/sfcgraph.x b/noao/onedspec/sensfunc/sfcgraph.x new file mode 100644 index 00000000..2f689d47 --- /dev/null +++ b/noao/onedspec/sensfunc/sfcgraph.x @@ -0,0 +1,104 @@ +include <gset.h> +include "sensfunc.h" + +# SF_CGRAPH -- Graph of composite points and errors + +procedure sf_cgraph (gp, stds, nstds, cv) + +pointer gp # Graphics structure +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +pointer cv # Sensitivity function curve + +int i, j, n, nwaves +real w, s, ymin, ymax, cveval() +double sum, sum2 +pointer sp, waves, sens, errors, xp, yp, zp, gio + +begin + nwaves = 0 + do i = 1, nstds-2 + if (STD_FLAG(stds[i]) != SF_EXCLUDE) + nwaves = nwaves + STD_NWAVES(stds[i]) + + call smark (sp) + call salloc (waves, nwaves, TY_REAL) + call salloc (sens, nwaves, TY_REAL) + call salloc (errors, nwaves, TY_REAL) + + nwaves = 0 + do i = 1, nstds-2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + n = STD_NWAVES(stds[i]) + xp = STD_WAVES(stds[i]) + yp = STD_SENS(stds[i]) + zp = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[zp] != 0.) { + Memr[waves+nwaves] = Memr[xp] + Memr[sens+nwaves] = Memr[yp] + nwaves = nwaves + 1 + } + xp = xp + 1 + yp = yp + 1 + zp = zp + 1 + } + } + call xt_sort2 (Memr[waves], Memr[sens], nwaves) + + n = 0 + sum = 0. + sum2 = 0. + ymin = 0. + ymax = 0. + j = 0 + do i = 1, nwaves { + w = Memr[waves+i-1] + s = Memr[sens+i-1] + n = n + 1 + sum = sum + s + sum2 = sum2 + s * s + + if ((i < nwaves) && (w == Memr[waves+i])) + next + + if (n > 2) { + sum = sum / n + sum2 = sum2 / n - sum * sum + if (sum2 > 0) + sum2 = sqrt (sum2 / n) + else + sum2 = 0. + sum = sum - cveval (cv, w) + + Memr[waves+j] = w + Memr[sens+j] = sum + Memr[errors+j] = sum2 + j = j + 1 + + if (sum + sum2 > ymax) + ymax = sum + sum2 + if (sum - sum2 < ymin) + ymin = sum - sum2 + } + n = 0 + sum = 0. + sum2 = 0. + } + nwaves = j + + if (j == 0) { + call printf ("No wavelength overlap for composite points") + } else { + gio = GP_GIO(gp) + call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax) + call glabax (gio, "Composite Points vs Wavelength", "", "") + call gseti (gio, G_PLCOLOR, GP_CMARK(gp)) + do i = 0, nwaves-1 + call gmark (gio, Memr[waves+i], Memr[sens+i], GM_VEBAR, + 1., -Memr[errors+i]) + } + + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sfcolon.x b/noao/onedspec/sensfunc/sfcolon.x new file mode 100644 index 00000000..43a056e7 --- /dev/null +++ b/noao/onedspec/sensfunc/sfcolon.x @@ -0,0 +1,193 @@ +include <error.h> +include <gset.h> +include "sensfunc.h" + +# SENSFUNC colon commands +define CMDS "|stats|vstats|function|order|graphs|images|skys|marks\ + |fluxlimits|colors|" +define STATS 1 # Show results +define VSTATS 2 # Show verbose results +define FUNCTION 3 # Sensitivity function type +define ORDER 4 # Function order +define GRAPHS 5 # Select graphs +define IMAGES 6 # Select images +define SKYS 7 # Select skys +define MARKS 8 # Set graph mark types +define FLIMITS 9 # Flux graph limits +define COLORS 10 # Flux graph limits + +# SF_COLON -- Process SENSFUNC colon commands. +# This procedure has so many arguments because of the STATS option. + +procedure sf_colon (cmd, gp, stds, nstds, cv, wextn, extn, nextn, ecv, function, + order, npts, rms, newfit, newgraph) + +char cmd[ARB] # Colon command +pointer gp # Graphics structure +pointer stds[nstds] # Standard star data +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 +char function[ARB] # Function type +int order # Function order +int npts # Number of points in fit +real rms # RMS in fit +int newfit # New function? +int newgraph # New graphs? + +int i, j, ncmd, ival, fd, nscan(), strdic(), open(), stridx() +real rval1, rval2 +bool streq() +pointer sp, str +errchk open + +begin + # Match the command against a dictionary. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call sscan (cmd) + call gargwrd (Memc[str], SZ_LINE) + ncmd = strdic (Memc[str], Memc[str], SZ_LINE, CMDS) + + # Switch on the command. + switch (ncmd) { + case STATS: + call gargwrd (Memc[str], SZ_LINE) + iferr { + # If no argument write to STDOUT otherwise append to file. + if (nscan() == 1) { + call gdeactivate (GP_GIO(gp), AW_CLEAR) + call sf_stats (STDOUT, stds, nstds, function, order, npts, + rms) + call greactivate (GP_GIO(gp), AW_PAUSE) + } else { + fd = open (Memc[str], APPEND, TEXT_FILE) + call sf_stats (fd, stds, nstds, function, order, npts, rms) + call close (fd) + } + } then + call erract (EA_WARN) + case VSTATS: + call gargwrd (Memc[str], SZ_LINE) + iferr { + if (nscan() == 1) { + # If no argument page on STDOUT otherwise append to file. + # A temp file is used in order to page output. + + call mktemp ("tmp$sf", Memc[str], SZ_LINE) + fd = open (Memc[str], NEW_FILE, TEXT_FILE) + call sf_stats (fd, stds, nstds, function, order, npts, rms) + call sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn, + ecv) + call close (fd) + call gpagefile (GP_GIO(gp), Memc[str], "sensfunc") + call delete (Memc[str]) + } else { + fd = open (Memc[str], APPEND, TEXT_FILE) + call sf_stats (fd, stds, nstds, function, order, npts, rms) + call sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn, + ecv) + call close (fd) + } + } then + call erract (EA_WARN) + case FUNCTION: + call gargwrd (Memc[str], SZ_LINE) + if (nscan() == 2) { + call strcpy (Memc[str], function, SZ_FNAME) + newfit = NO + } else { + call printf ("function %s") + call pargstr (function) + } + case ORDER: + call gargi (ival) + if (nscan() == 2) { + order = ival + newfit = NO + } else { + call printf ("order %d") + call pargi (order) + } + case GRAPHS: + call gargstr (Memc[str], SZ_LINE) + j = str + for (i=str; Memc[i] != EOS; i=i+1) { + switch (Memc[i]) { + case 'a','c','e','i','l','r','s': + Memc[j] = Memc[i] + j = j + 1 + } + } + Memc[j] = EOS + if (Memc[str] != EOS) { + call strcpy (Memc[str], GP_GRAPHS(gp,1), SF_NGRAPHS) + newgraph = YES + } else { + call printf ("graphs %s") + call pargstr (GP_GRAPHS(gp,1)) + } + case IMAGES: + # Note that changing the image automatically clears the sky. + do i = 1, SF_NGRAPHS { + call gargwrd (Memc[str], SZ_LINE) + if (nscan() == i + 1) { + call strcpy (Memc[str], Memc[GP_IMAGES(gp,i)], SZ_FNAME) + Memc[GP_SKYS(gp,i)] = EOS + do j = 1, nstds + if (streq (Memc[str], STD_IMAGE(stds[j]))) + call strcpy (STD_SKY(stds[j]), Memc[GP_SKYS(gp,i)], + SZ_FNAME) + } else + break + } + if (nscan() == 1) { + call printf ("images %s %s %s %s") + call pargstr (Memc[GP_IMAGES(gp,1)]) + call pargstr (Memc[GP_IMAGES(gp,2)]) + call pargstr (Memc[GP_IMAGES(gp,3)]) + call pargstr (Memc[GP_IMAGES(gp,4)]) + } + case SKYS: + do i = 1, SF_NGRAPHS { + call gargwrd (Memc[str], SZ_LINE) + if (nscan() == i + 1) + call strcpy (Memc[str], Memc[GP_SKYS(gp,i)], SZ_FNAME) + else + break + } + if (nscan() == 1) { + call printf ("skys %s %s %s %s") + call pargstr (Memc[GP_SKYS(gp,1)]) + call pargstr (Memc[GP_SKYS(gp,2)]) + call pargstr (Memc[GP_SKYS(gp,3)]) + call pargstr (Memc[GP_SKYS(gp,4)]) + } + case MARKS: + call gargstr (Memc[str], SZ_LINE) + call sf_marks (gp, Memc[str]) + case FLIMITS: + call gargr (rval1) + call gargr (rval2) + if (nscan() == 3) { + GP_FMIN(gp) = rval1 + GP_FMAX(gp) = rval2 + if (stridx (GP_GRAPHS(gp,1), "il") != 0) + newgraph = YES + } else { + call printf ("fluxlimits %g %g") + call pargr (GP_FMIN(gp)) + call pargr (GP_FMAX(gp)) + } + case COLORS: + call gargstr (Memc[str], SZ_LINE) + call sf_colors (gp, Memc[str]) + default: + call printf ("Unrecognized or ambiguous command\007") + } + + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sfcolors.x b/noao/onedspec/sensfunc/sfcolors.x new file mode 100644 index 00000000..db3f69af --- /dev/null +++ b/noao/onedspec/sensfunc/sfcolors.x @@ -0,0 +1,28 @@ +include <gset.h> +include "sensfunc.h" + + +# SF_COLORS -- Set colors. + +procedure sf_colors (gp, colors) + +pointer gp +char colors[ARB] + +int i, nscan() + +begin + call sscan (colors) + call gargi (i) + if (nscan() == 1) + GP_PLCOLOR(gp) = i + call gargi (i) + if (nscan() == 2) + GP_CMARK(gp) = i + call gargi (i) + if (nscan() == 3) + GP_CDEL(gp) = i + call gargi (i) + if (nscan() == 4) + GP_CADD(gp) = i +end diff --git a/noao/onedspec/sensfunc/sfcomposite.x b/noao/onedspec/sensfunc/sfcomposite.x new file mode 100644 index 00000000..506416e3 --- /dev/null +++ b/noao/onedspec/sensfunc/sfcomposite.x @@ -0,0 +1,147 @@ +include "sensfunc.h" + +# SF_COMPOSITE -- Create a composite standard structure. +# The composite star is the last of the standard stars. +# When the composite star is created the other stars are turned off. +# The function toggles. + +procedure sf_composite (stds, nstds, cv) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +pointer cv # Sensitivity pointer + +int i, j, k, n, nwaves +pointer std, waves, sens, fit, wts, iwts, x, y, z +errchk malloc, realloc, xt_sort3 + +begin + # If data is already composite toggle back to original data. + # Delete data points if composite point is deleted. + std = stds[nstds] + if (STD_FLAG(std) == SF_INCLUDE) { + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + STD_FLAG(stds[i]) = SF_INCLUDE + } + STD_FLAG(std) = SF_EXCLUDE + + n = STD_NWAVES(std) + x = STD_WAVES(std) + z = STD_WTS(std) + do i = 1, n { + if (Memr[z] == 0.) { + do j = 1, nstds - 2 { + if (STD_FLAG(stds[j]) == SF_EXCLUDE) + next + nwaves = STD_NWAVES(stds[j]) + waves = STD_WAVES(stds[j]) + wts = STD_WTS(stds[j]) + do k = 1, nwaves { + if (Memr[waves] == Memr[x]) + Memr[wts] = 0. + waves = waves + 1 + wts = wts + 1 + } + } + } + x = x + 1 + z = z + 1 + } + call printf ("Individual star data") + return + } + + # Initialize + if (STD_WAVES(std) != NULL) { + call mfree (STD_WAVES(std), TY_REAL) + call mfree (STD_SENS(std), TY_REAL) + call mfree (STD_WTS(std), TY_REAL) + call mfree (STD_IWTS(std), TY_REAL) + call mfree (STD_X(std), TY_REAL) + call mfree (STD_Y(std), TY_REAL) + } + + # To bin the data we collect all the data and then sort by wavelength. + nwaves = 0 + do i = 1, nstds - 2 + if (STD_FLAG(stds[i]) == SF_INCLUDE) + nwaves = nwaves + STD_NWAVES(stds[i]) + + call malloc (waves, nwaves, TY_REAL) + call malloc (sens, nwaves, TY_REAL) + call malloc (wts, nwaves, TY_REAL) + + nwaves = 0 + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + x = STD_WAVES(stds[i]) + y = STD_SENS(stds[i]) + z = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[z] != 0.) { + Memr[waves+nwaves] = Memr[x] + Memr[sens+nwaves] = Memr[y] + Memr[wts+nwaves] = Memr[z] + nwaves = nwaves + 1 + } + x = x + 1 + y = y + 1 + z = z + 1 + } + STD_FLAG(stds[i]) = SF_DELETE + STD_BEAM(std) = STD_BEAM(stds[i]) + STD_WSTART(std) = STD_WSTART(stds[i]) + STD_WEND(std) = STD_WEND(stds[i]) + } +# STD_NWAVES(stds[nstds-1]) = 0 + + call xt_sort3 (Memr[waves], Memr[sens], Memr[wts], nwaves) + + # Go through the wavelength sorted data and composite all points + # with the same wavelength. + + n = 0 + Memr[sens] = Memr[wts] * Memr[sens] + do i = 1, nwaves-1 { + if (Memr[waves+i] == Memr[waves+n]) { + Memr[sens+n] = Memr[sens+n] + Memr[wts+i] * Memr[sens+i] + Memr[wts+n] = Memr[wts+n] + Memr[wts+i] + } else { + n = n + 1 + Memr[waves+n] = Memr[waves+i] + Memr[sens+n] = Memr[wts+i] * Memr[sens+i] + Memr[wts+n] = Memr[wts+i] + } + } + + nwaves = n + 1 + do i = 0, nwaves-1 + Memr[sens+i] = Memr[sens+i] / Memr[wts+i] + + # Store the composite data in the standard star structure. + call realloc (waves, nwaves, TY_REAL) + call realloc (sens, nwaves, TY_REAL) + call realloc (wts, nwaves, TY_REAL) + call malloc (iwts, nwaves, TY_REAL) + call malloc (fit, nwaves, TY_REAL) + call malloc (x, nwaves, TY_REAL) + call malloc (y, nwaves, TY_REAL) + call amovr (Memr[wts], Memr[iwts], nwaves) + call cvvector (cv, Memr[waves], Memr[fit], nwaves) + + STD_FLAG(std) = SF_INCLUDE + STD_NWAVES(std) = nwaves + STD_WAVES(std) = waves + STD_SENS(std) = sens + STD_FIT(std) = fit + STD_WTS(std) = wts + STD_IWTS(std) = iwts + STD_X(std) = x + STD_Y(std) = y + + call printf ("Composite star data") +end diff --git a/noao/onedspec/sensfunc/sfdata.x b/noao/onedspec/sensfunc/sfdata.x new file mode 100644 index 00000000..94140049 --- /dev/null +++ b/noao/onedspec/sensfunc/sfdata.x @@ -0,0 +1,59 @@ +include "sensfunc.h" + +# SF_DATA -- Compute the X and Y values for the particular graph. + +procedure sf_data (stds, nstds, graph) + +pointer stds[nstds] # Standard star structures +int nstds # Number of standard stars +char graph # Graph type + +real a +int i, n +pointer wp, sp, fp, xp, yp + +begin + switch (graph) { + case 's': # Sensitivity vs. Wavelength + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + wp = STD_WAVES(stds[i]) + sp = STD_SENS(stds[i]) + xp = STD_X(stds[i]) + yp = STD_Y(stds[i]) + call amovr (Memr[wp], Memr[xp], n) + call amovr (Memr[sp], Memr[yp], n) + } + case 'a': # Residuals vs. Airmass + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + wp = STD_WAVES(stds[i]) + sp = STD_SENS(stds[i]) + fp = STD_FIT(stds[i]) + xp = STD_X(stds[i]) + yp = STD_Y(stds[i]) + call amovkr (a, Memr[xp], n) + call asubr (Memr[sp], Memr[fp], Memr[yp], n) + } + case 'r': # Residuals vs. Wavelength + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + wp = STD_WAVES(stds[i]) + sp = STD_SENS(stds[i]) + fp = STD_FIT(stds[i]) + xp = STD_X(stds[i]) + yp = STD_Y(stds[i]) + call amovr (Memr[wp], Memr[xp], n) + call asubr (Memr[sp], Memr[fp], Memr[yp], n) + } + } +end diff --git a/noao/onedspec/sensfunc/sfdelete.x b/noao/onedspec/sensfunc/sfdelete.x new file mode 100644 index 00000000..ff2d267c --- /dev/null +++ b/noao/onedspec/sensfunc/sfdelete.x @@ -0,0 +1,127 @@ +include <gset.h> +include "sensfunc.h" + +# SF_DELETE -- Delete point, star, or wavelength identified by the +# star index and index within the array of values. + +procedure sf_delete (gp, stds, nstds, key, istd, ipt) + +pointer gp # GIO pointer +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +int key # Delete point, star, or wavelength +int istd # Index of standard star +int ipt # Index of point + +int i, j, n, wcs, mark, mdel, cdel, stridx() +real wave, szmark, szmdel +pointer x, y, z, w, gio + +begin + gio = GP_GIO(gp) + mdel = GP_MDEL(gp) + cdel = GP_CDEL(gp) + szmdel = GP_SZMDEL(gp) + szmark = GP_SZMARK(gp) + + # Delete the point or points from each displayed graph. + # When deleting multiple points check if point already deleted. + for (wcs = 1; GP_GRAPHS(gp,wcs) != EOS; wcs = wcs + 1) { + if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0) + next + + call gseti (gio, G_WCS, wcs) + call sf_data (stds, nstds, GP_GRAPHS(gp,wcs)) + switch (key) { + case 'p': + if (istd != nstds-1) + mark = GP_MARK(gp) + else + mark = GP_MADD(gp) + x = STD_X(stds[istd])+ipt-1 + y = STD_Y(stds[istd],1)+ipt-1 + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel) + case 's': + if (istd != nstds-1) + mark = GP_MARK(gp) + else + mark = GP_MADD(gp) + n = STD_NWAVES(stds[istd]) + x = STD_X(stds[istd]) + y = STD_Y(stds[istd]) + w = STD_WTS(stds[istd]) + do i = 1, n { + if (Memr[w] != 0.) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel) + } + x = x + 1 + y = y + 1 + w = w + 1 + } + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + if (i != nstds-1) + mark = GP_MARK(gp) + else + mark = GP_MADD(gp) + n = STD_NWAVES(stds[i]) + x = STD_X(stds[i]) + y = STD_Y(stds[i]) + z = STD_WAVES(stds[i]) + w = STD_WTS(stds[i]) + do j = 1, n { + if ((Memr[z] == wave) && (Memr[w] != 0.)) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mark, szmark, + szmark) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, + szmdel) + } + x = x + 1 + y = y + 1 + z = z + 1 + w = w + 1 + } + } + } + } + + # Mark the points as deleted by setting their weights to zero. + switch (key) { + case 'p': + w = STD_WTS(stds[istd])+ipt-1 + Memr[w] = 0. + case 's': + n = STD_NWAVES(stds[istd]) + w = STD_WTS(stds[istd]) + call aclrr (Memr[w], n) + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + z = STD_WAVES(stds[i]) + w = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[z] == wave) + Memr[w] = 0. + w = w + 1 + z = z + 1 + } + } + } +end diff --git a/noao/onedspec/sensfunc/sfeout.x b/noao/onedspec/sensfunc/sfeout.x new file mode 100644 index 00000000..8dae6301 --- /dev/null +++ b/noao/onedspec/sensfunc/sfeout.x @@ -0,0 +1,114 @@ +include <error.h> +include <ctype.h> +include <mach.h> +include <math/curfit.h> + +define NEPOINTS 100 # Number of points in extinction table + +# SF_EOUT -- Output a revised extinction table. This is only done if there +# is at least one residual extinction curve. No assumption is made about +# overlapping extinction curves. In the overlap the extinction corrections +# are averaged. + +procedure sf_eout (wextn, extn, nextn, ecvs, necvs) + +real wextn[nextn] # Standard extinction wavelengths +real extn[nextn] # Standard extinction values +int nextn # Number of standard extinction points +pointer ecvs[necvs] # Residual extinction curves (one for each beam) +int necvs # Number of residual extinction curves + +int i, j, fd, open(), scan(), nscan() +real w, ext, wmin, wmax, dw, xmin, xmax, cvstatr(), cveval() +pointer sp, fname, waves, extns, navg, cv + +define newfile_ 99 + +begin + # If there are no residual extinction values then do nothing. + for (i=1; (i<=necvs) && (ecvs[i]==NULL); i=i+1) + ; + if (i > necvs) + return + + # The output table consists of NEPOINTS. + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (waves, NEPOINTS, TY_REAL) + call salloc (extns, NEPOINTS, TY_REAL) + call salloc (navg, NEPOINTS, TY_INT) + call aclrr (Memr[extns], NEPOINTS) + call aclri (Memi[navg], NEPOINTS) + + # Open the extinction table. If it fails allow the user to + # enter a new name. + + call clgstr ("newextinction", Memc[fname], SZ_FNAME) + for (i=fname; (Memc[i]!=EOS) && IS_WHITE(Memc[i]); i=i+1) + if (Memc[i] == EOS) { + call sfree (sp) + return + } + +newfile_ + iferr (fd = open (Memc[i], NEW_FILE, TEXT_FILE)) { + call printf ("Cannot create %s -- Enter new extinction filename: ") + call pargstr (Memc[fname]) + call flush (STDOUT) + if (scan() != EOF) { + call gargwrd (Memc[fname], SZ_FNAME) + if (nscan() == 1) + goto newfile_ + } + call sfree (sp) + call printf ("No revised extinction file created\n") + return + } + + # Determine the range of the extinction table. + wmin = MAX_REAL + wmax = -MAX_REAL + do i = 1, necvs { + if (ecvs[i] == NULL) + next + wmin = min (wmin, cvstatr (ecvs[i], CVXMIN)) + wmax = max (wmax, cvstatr (ecvs[i], CVXMAX)) + } + dw = (wmax - wmin) / (NEPOINTS - 1) + do i = 1, NEPOINTS + Memr[waves+i-1] = wmin + (i - 1) * dw + + # Average the residual extinctions and add the original extinction. + do j = 1, necvs { + if (ecvs[j] == NULL) + next + cv = ecvs[j] + xmin = cvstatr (cv, CVXMIN) + xmax = cvstatr (cv, CVXMAX) + do i = 1, NEPOINTS { + w = Memr[waves+i-1] + if ((w < xmin) || (w > xmax)) + next + Memr[extns+i-1] = Memr[extns+i-1] + cveval (cv, w) + Memi[navg+i-1] = Memi[navg+i-1] + 1 + } + } + do i = 1, NEPOINTS { + if (Memi[navg+i-1] > 0) + Memr[extns+i-1] = Memr[extns+i-1] / Memi[navg+i-1] + w = Memr[waves+i-1] + call intrp (1, wextn, extn, nextn, w, ext, j) + Memr[extns+i-1] = Memr[extns+i-1] + ext + } + + # Output extinction table. + call fprintf (fd, "# Revised extinction table.\n") + do i = 1, NEPOINTS { + call fprintf (fd, "%7.2f %6.3f\n") + call pargr (Memr[waves+i-1]) + call pargr (Memr[extns+i-1]) + } + call close (fd) + + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sfextinct.x b/noao/onedspec/sensfunc/sfextinct.x new file mode 100644 index 00000000..f7b95326 --- /dev/null +++ b/noao/onedspec/sensfunc/sfextinct.x @@ -0,0 +1,226 @@ +include <pkg/gtools.h> +include "sensfunc.h" + +define RANGE_AIRMASS 0.1 # Minimum airmass range +define SIGMA_AIRMASS 0.05 # Minimum sigma in airmass + +# SF_EXINCT -- Determine a residual extinction curve. At each wavelength +# for which there are multiple observations or neighboring wavelengths +# such that there is a sufficient airmass range determine the slope +# of the sensitivity vs airmass. Residual sensitivity is used to minimize +# wavelength scatter when multiple wavelengths are needed because of +# nonoverlapping standard star data. Each such slope is a measure of the +# residual extinction at that wavelength. To make the residual extinction +# curve fit the extinction vs. wavelength using ICFIT. + +procedure sf_extinct (gp, stds, nstds, cv, ecv, function, order) + +pointer gp # Graphics structure +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +pointer cv # Sensitivity function curve +pointer ecv # Residual extinction curve +char function[ARB] # Fitting function +int order # Function order + +bool ans +int i, j, n, nwaves, sum, npts, scan() +real a, amin, amax, rms, rms1, r2, sig, cveval() +double x, y, z, sumx, sumy, sumz, sumx2, sumxy +pointer sp, waves, sens, airm, xp, yp, fp, wp, ic +pointer gt, gt_init() +errchk salloc, xt_sort3, icg_fit, ic_open + +define cancel_ 99 + +begin + # Cancel previous extinction if defined. + if (ecv != NULL) + goto cancel_ + + # Check for minimum airmass range and determine number of points. + # Ignore added objects and composite data. + amin = 100. + amax = 0. + nwaves = 0 + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + nwaves = nwaves + STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + amin = min (amin, a) + amax = max (amax, a) + } + if (amax - amin < RANGE_AIRMASS) { + call printf ( + "Insufficient airmass coverage for extinction determination") + return + } + + # Extract data to be used and sort by wavelength. + # The data is wavelength, airmass, and residual sensitivity. + call smark (sp) + call salloc (waves, nwaves, TY_REAL) + call salloc (sens, nwaves, TY_REAL) + call salloc (airm, nwaves, TY_REAL) + + nwaves = 0 + do i = 1, nstds-2 { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + xp = STD_WAVES(stds[i]) + yp = STD_SENS(stds[i]) + fp = STD_FIT(stds[i]) + wp = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[wp] != 0.) { + Memr[airm+nwaves] = a + Memr[waves+nwaves] = Memr[xp] + Memr[sens+nwaves] = Memr[yp] - Memr[fp] + nwaves = nwaves + 1 + } + xp = xp + 1 + yp = yp + 1 + fp = fp + 1 + wp = wp + 1 + } + } + + call xt_sort3 (Memr[waves], Memr[sens], Memr[airm], nwaves) + + # Bin points with common wavelengths or at least two points. + sum = 0 + sumx = 0. + sumy = 0. + sumz = 0. + sumx2 = 0. + sumxy = 0. + n = 0 + do i = 0, nwaves-2 { + x = Memr[airm+i] + y = Memr[sens+i] + z = Memr[waves+i] + sum = sum + 1 + sumx = sumx + x + sumy = sumy + y + sumx2 = sumx2 + x * x + sumxy = sumxy + x * y + sumz = sumz + z + + if ((z == Memr[waves+i+1]) || (sum < 2)) + next + + x = sumx2 - sumx * sumx / sum + if (x > SIGMA_AIRMASS) { + Memr[waves+n] = sumz / sum + Memr[sens+n] = (sumx * sumy / sum - sumxy) / x + Memr[airm+n] = 1. + n = n + 1 + sum = 0 + sumx = 0. + sumy = 0. + sumz = 0. + sumx2 = 0. + sumxy = 0. + } + } + if (sum > 1) { + x = sumx2 - sumx * sumx / sum + if (x > SIGMA_AIRMASS) { + Memr[waves+n] = sumz / sum + Memr[sens+n] = (sumx * sumy / sum - sumxy) / x + Memr[airm+n] = 1. + n = n + 1 + } + } + + if (n < 2) { + call printf ("Cannot determine residual extinction") + call sfree (sp) + return + } + + # Fit residual extinction curve using ICFIT. + gt = gt_init() + call gt_sets (gt, GTTYPE, "mark") + call gt_seti (gt, GTCOLOR, GP_PLCOLOR(gp)) + call ic_open (ic) + call ic_putr (ic, "xmin", min (STD_WSTART(stds[1]), STD_WEND(stds[1]))) + call ic_putr (ic, "xmax", max (STD_WSTART(stds[1]), STD_WEND(stds[1]))) + call ic_pstr (ic, "function", "chebyshev") + call ic_puti (ic, "order", 1) + call ic_pstr (ic, "xlabel", "wavelength") + call ic_pstr (ic, "ylabel", "residual extinction") + call ic_pstr (ic, "yunits", "mag") + call icg_fit (ic, GP_GIO(gp), "cursor", gt, ecv, Memr[waves], + Memr[sens], Memr[airm], n) + call ic_closer (ic) + call gt_free (gt) + + # Determine significance of the fit. + call sf_fit (stds, nstds, cv, function, order, + min (GP_WSTART(gp), GP_WEND(gp)), max (GP_WSTART(gp), GP_WEND(gp))) + call sf_rms (stds, nstds, rms1, npts) + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + xp = STD_WAVES(stds[i]) + yp = STD_SENS(stds[i]) + call cvvector (ecv, Memr[xp], Memr[sens], n) + call amulkr (Memr[sens], a, Memr[sens], n) + call aaddr (Memr[yp], Memr[sens], Memr[yp], n) + } + call sf_fit (stds, nstds, cv, function, order, + min (GP_WSTART(gp), GP_WEND(gp)), max (GP_WSTART(gp), GP_WEND(gp))) + call sf_rms (stds, nstds, rms, npts) + do i = 1, SF_NGRAPHS + if (GP_SHDR(gp,i) != NULL) + call shdr_close (GP_SHDR(gp,i)) + + r2 = 1 - rms ** 2 / rms1 ** 2 + sig = r2 * (nwaves - 2) / max (0.01, 1. - r2) + if (sig <= 0.0) + sig = 0. + else + sig = sqrt (sig) + + # Apply to data if desired. + call printf ( + "Significance = %4.1f sigma: Apply residual extinction correction? ") + call pargr (sig) + call flush (STDOUT) + + ans = false + if (scan() != EOF) + call gargb (ans) + + # Undo last fit if not applying correction. + if (!ans) + goto cancel_ + + call printf ("Residual extinction correction applied") + call sfree (sp) + return + +cancel_ + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + n = STD_NWAVES(stds[i]) + a = STD_AIRMASS(stds[i]) + xp = STD_WAVES(stds[i]) + yp = STD_SENS(stds[i]) + do j = 1, n { + Memr[yp] = Memr[yp] - a * cveval (ecv, Memr[xp]) + xp = xp + 1 + yp = yp + 1 + } + } + call cvfree (ecv) + call printf ("Residual extinction correction canceled") + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sffit.x b/noao/onedspec/sensfunc/sffit.x new file mode 100644 index 00000000..3be306ad --- /dev/null +++ b/noao/onedspec/sensfunc/sffit.x @@ -0,0 +1,78 @@ +include <math/curfit.h> +include "sensfunc.h" + +define FUNCTIONS "|chebyshev|legendre|spline3|spline1|" + +procedure sf_fit (stds, nstds, cv, function, order, xmin, xmax) + +pointer stds[nstds] +int nstds +pointer cv +char function[ARB] +int order +real xmin +real xmax + +int i, n, func, ord, strdic() +pointer x, y, w + +int functions[4] +data functions/CHEBYSHEV,LEGENDRE,SPLINE3,SPLINE1/ + +begin + func = strdic (function, function, SZ_FNAME, FUNCTIONS) + func = functions[max (1, func)] + ord = order + + while (ord > 0) { + call cvfree (cv) + call cvinit (cv, func, ord, xmin, xmax) + do i = 1, nstds { + if (STD_FLAG(stds[i]) == SF_INCLUDE) { + n = STD_NWAVES(stds[i]) + x = STD_WAVES(stds[i]) + y = STD_SENS(stds[i]) + w = STD_WTS(stds[i]) + call cvacpts (cv, Memr[x], Memr[y], Memr[w], n, WTS_USER) + } + } + call cvsolve (cv, i) + if (i == OK) + break + + switch (func) { + case CHEBYSHEV, LEGENDRE: + ord = ord - 1 + case SPLINE3: + ord = ord - 1 + if (ord == 0) { + func = CHEBYSHEV + ord = 2 + } + case SPLINE1: + ord = ord - 1 + if (ord == 0) { + func = CHEBYSHEV + ord = 1 + } + } + } + + switch (i) { + case SINGULAR: + call error (0, "Singular solution") + case NO_DEG_FREEDOM: + call error (0, "No degrees of freedom") + } + + # Set fitted values + do i = 1, nstds + if (STD_FLAG(stds[i]) != SF_EXCLUDE) { + n = STD_NWAVES(stds[i]) + if (n < 1) + next + x = STD_WAVES(stds[i]) + y = STD_FIT(stds[i]) + call cvvector (cv, Memr[x], Memr[y], n) + } +end diff --git a/noao/onedspec/sensfunc/sfginit.x b/noao/onedspec/sensfunc/sfginit.x new file mode 100644 index 00000000..0214c7a7 --- /dev/null +++ b/noao/onedspec/sensfunc/sfginit.x @@ -0,0 +1,89 @@ +include <gset.h> +include "sensfunc.h" + +# SF_GINIT -- Initialize graphics structure and open the graphics device. +# This includes CL requests for the starting graphs (default is "sr"), +# the mark types (default is "plus box"), and graphics device. + +procedure sf_ginit (gp) + +pointer gp # Graphics structure (returned) + +int i, j +pointer sp, str, gopen() +errchk malloc, gopen + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + call calloc (gp, LEN_GP, TY_STRUCT) + do i = 1, SF_NGRAPHS { + call malloc (GP_IMAGES(gp,i), SZ_FNAME, TY_CHAR) + Memc[GP_IMAGES(gp,i)] = EOS + call malloc (GP_SKYS(gp,i), SZ_FNAME, TY_CHAR) + Memc[GP_SKYS(gp,i)] = EOS + } + + # Set the starting graph types. + call clgstr ("graphs", Memc[str], SZ_FNAME) + j = str + for (i=str; Memc[i] != EOS; i=i+1) { + switch (Memc[i]) { + case 'a','c','e','i','l','r','s': + Memc[j] = Memc[i] + j = j + 1 + } + } + Memc[j] = EOS + if (Memc[str] != EOS) + call strcpy (Memc[str], GP_GRAPHS(gp,1), SF_NGRAPHS) + else + call strcpy ("sr", GP_GRAPHS(gp,1), SF_NGRAPHS) + + # Set the starting mark types and colors. + GP_MARK(gp) = GM_PLUS + GP_MDEL(gp) = GM_CROSS + GP_MADD(gp) = GM_BOX + GP_PLCOLOR(gp) = 2 + GP_CMARK(gp) = 1 + GP_CDEL(gp) = 3 + GP_CADD(gp) = 4 + call clgstr ("marks", Memc[str], SZ_FNAME) + call sf_marks (gp, Memc[str]) + call clgstr ("colors", Memc[str], SZ_FNAME) + call sf_colors (gp, Memc[str]) + + # Set flux limits + GP_FMIN(gp) = INDEF + GP_FMAX(gp) = INDEF + + # Open the graphics device. + call clgstr ("device", Memc[str], SZ_FNAME) + GP_GIO(gp) = gopen (Memc[str], NEW_FILE, STDGRAPH) + + call sfree (sp) +end + + +# SF_GFREE -- Free the graphics structure and close the graphics device. + +procedure sf_gfree (gp) + +pointer gp # Graphics structure + +int i + +begin + if (gp == NULL) + return + + call gclose (GP_GIO(gp)) + do i = 1, SF_NGRAPHS { + call mfree (GP_IMAGES(gp,i), TY_CHAR) + call mfree (GP_SKYS(gp,i), TY_CHAR) + if (GP_SHDR(gp,i) != NULL) + call shdr_close (GP_SHDR(gp,i)) + } + call mfree (gp, TY_STRUCT) +end diff --git a/noao/onedspec/sensfunc/sfgraph.x b/noao/onedspec/sensfunc/sfgraph.x new file mode 100644 index 00000000..bb6fb26f --- /dev/null +++ b/noao/onedspec/sensfunc/sfgraph.x @@ -0,0 +1,289 @@ +include <gset.h> +include <error.h> +include <math/curfit.h> +include "sensfunc.h" + +define NCURVE 50 # Number of vectors in curve + + +# SF_GRAPH -- Graph the desired data on the output graphics device. +# This entry procedure determines the graph types, sets the device viewports +# for each graph, and calls a procedure to make the graph. + +procedure sf_graph (gp, stds, nstds, cv, wextn, extn, nextn, ecv) + +pointer gp # Graphics structure +pointer stds[nstds] # Standard star data +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 values in extinction table +pointer ecv # Residual extinction curve. + +int i, image, ngraphs, strlen() +real fa[8] +pointer sp, id, gio + +data fa/0.,1.,1.,0.,0.,0.,1.,1./ + +begin + # Clear the graphs, write the title, and set the viewports based on + # the number of graphs. + + call smark (sp) + call salloc (id, SZ_LINE, TY_CHAR) + call sysid (Memc[id], SZ_LINE) + + gio = GP_GIO(gp) + call gclear (gio) + call gseti (gio, G_FACOLOR, 0) + call gseti (gio, G_NMINOR, 0) + ngraphs = strlen (GP_GRAPHS(gp,1)) + switch (ngraphs) { + case 1: + call gseti (gio, G_WCS, 1) + call gsview (gio, .10, .97, .1, .9) + GP_SZMARK(gp) = 2. + GP_SZMDEL(gp) = 2. + case 2: + call gseti (gio, G_WCS, 1) + call gsview (gio, .10, .97, .55, .9) + call gseti (gio, G_WCS, 2) + call gsview (gio, .10, .97, .10, .45) + GP_SZMARK(gp) = 2. + GP_SZMDEL(gp) = 2. + case 3: + call gseti (gio, G_WCS, 1) + call gsview (gio, .10, .97, .55, .9) + call gseti (gio, G_WCS, 2) + call gsview (gio, .10, .50, .1, .45) + call gseti (gio, G_WCS, 3) + call gsview (gio, .57, .97, .1, .45) + GP_SZMARK(gp) = 2. + GP_SZMDEL(gp) = 2. + default: + call gseti (gio, G_WCS, 1) + call gsview (gio, .10, .50, .55, .9) + call gseti (gio, G_WCS, 2) + call gsview (gio, .57, .97, .55, .9) + call gseti (gio, G_WCS, 3) + call gsview (gio, .10, .50, .1, .45) + call gseti (gio, G_WCS, 4) + call gsview (gio, .57, .97, .1, .45) + GP_SZMARK(gp) = .01 + GP_SZMDEL(gp) = .01 + } + + # For each graph select the viewport and call a procedure to make + # the graph. + + image = 0 + for (i = 1; GP_GRAPHS(gp,i) != EOS; i = i + 1) { + call gseti (gio, G_WCS, i) + if (i > 1) + call gfill (gio, fa, fa[5], 4, GF_SOLID) + switch (GP_GRAPHS(gp,i)) { + case 'a', 's', 'r': + call sf_data (stds, nstds, GP_GRAPHS(gp,i)) + call sf_graph1 (gp, stds, nstds, GP_GRAPHS(gp,i)) + case 'e': + call sf_egraph (gp, wextn, extn, nextn, ecv) + case 'c': + call sf_cgraph (gp, stds, nstds, cv) + case 'i', 'l': + if (GP_GRAPHS(gp,i) == 'i') + GP_LOG(gp) = NO + else + GP_LOG(gp) = YES + image = image + 1 + iferr (call sf_image (gp, image, stds, nstds, cv, wextn, extn, + nextn, ecv)) + call erract (EA_WARN) + } + } + + call gseti (gio, G_WCS, 0) + call gtext (gio, 0.5, 1., Memc[id], "h=c,v=t,f=b") + call gtext (gio, 0.5, 0.97, GP_TITLE(gp), "h=c,v=t,f=b") + + call sfree (sp) +end + + +# SF_GRAPH1 -- Make graph of sensitivity or residual sensitivity vs wavelength. + +procedure sf_graph1 (gp, stds, nstds, graph) + +pointer gp # Graphics structure +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +char graph # Graph type + +int i, j, n, mark, mdel, cdel, color +real szmark, szmdel, ymin, ymax, y1, y2 +pointer x, y, w, gio + +begin + gio = GP_GIO(gp) + + # Autoscale the included data in Y and set wavelength range. + j = 0 + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + j = j + 1 + n = STD_NWAVES(stds[i]) + x = STD_X(stds[i]) + y = STD_Y(stds[i]) + if (j == 1) + call alimr (Memr[y], n, ymin, ymax) + else { + call alimr (Memr[y], n, y1, y2) + ymin = min (ymin, y1) + ymax = max (ymax, y2) + } + } + y1 = 0.05 * (ymax - ymin) + ymin = ymin - y1 + ymax = ymax + y1 + + # Draw axes and title based on type of graph. + switch (graph) { + case 'a': + call gswind (gio, GP_AIRMASS(gp,1), GP_AIRMASS(gp,2), ymin, ymax) + call glabax (gio, "Sensitivity Residuals vs Airmass", "", "") + case 's': + call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax) + call glabax (gio, "Sensitivity vs Wavelength", "", "") + case 'r': + call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax) + call glabax (gio, "Sensitivity Residuals vs Wavelength", "", "") + } + + # Plot the data with appropriate mark types and sizes. + mdel = GP_MDEL(gp) + cdel = GP_CDEL(gp) + szmdel = GP_SZMDEL(gp) + szmark = GP_SZMARK(gp) + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + if (i != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + n = STD_NWAVES(stds[i]) + x = STD_X(stds[i]) - 1 + y = STD_Y(stds[i]) - 1 + w = STD_WTS(stds[i]) - 1 + do j = 1, n { + if (Memr[w+j] == 0.) { + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x+j], Memr[y+j], mdel, szmdel, szmdel) + } else { + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x+j], Memr[y+j], mark, szmark, szmark) + } + } + } +end + + +# SF_EGRAPH -- Graph extinction curves with and without residual extinction +# correction. + +procedure sf_egraph (gp, wextn, extn, nextn, ecv) + +pointer gp # Graphics structure +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 i, j +real xmin, xmax, dx, x, cveval() +pointer sp, ext, extnew, gio + +begin + call smark (sp) + call salloc (ext, NCURVE, TY_REAL) + + # Interpolate extinction table to a grid of wavelengths within + # the range of the sensitivity data. + + gio = GP_GIO(gp) + xmin = GP_WSTART(gp) + xmax = GP_WEND(gp) + dx = (xmax - xmin) / (NCURVE - 1) + x = xmin + do i = 0, NCURVE-1 { + call intrp (1, wextn, extn, nextn, x, Memr[ext+i], j) + x = x + dx + } + call gswind (gio, xmin, xmax, INDEF, INDEF) + call gascale (gio, Memr[ext], NCURVE, 2) + + # If there is a residual extinction curve determine a new extinction + # curve. + + if (ecv != NULL) { + call salloc (extnew, NCURVE, TY_REAL) + call amovr (Memr[ext], Memr[extnew], NCURVE) + x = xmin + do i = 0, NCURVE-1 { + Memr[extnew+i] = Memr[extnew+i] + cveval (ecv, x) + x = x + dx + } + call grscale (gio, Memr[extnew], NCURVE, 2) + } + + # Draw the axes and title and extinction curves. + call glabax (gio, "Extinction vs Wavelength", "", "") + call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp)) + call gvline (gio, Memr[ext], NCURVE, xmin, xmax) + if (ecv != NULL) { + call gseti (gio, G_PLTYPE, 2) + call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp)+1) + call gvline (gio, Memr[extnew], NCURVE, xmin, xmax) + call gseti (gio, G_PLTYPE, 1) + } + + call sfree (sp) +end + + +# SF_FITGRAPH -- Overplot the fitted sensitivity curve. + +procedure sf_fitgraph (gp, cv) + +pointer gp # Graphics structure +pointer cv # Sensitivity function curve + +int i, j +real x1, x2, y1, y2, cveval() +pointer gio + +begin + gio = GP_GIO(gp) + call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp)) + + # Only plot on sensitivity curve graph types. + for (j = 1; GP_GRAPHS(gp,j) != EOS; j = j + 1) { + if (GP_GRAPHS(gp,j) != 's') + next + call gseti (gio, G_WCS, j) + call ggwind (gio, x1, x2, y1, y2) + x2 = (x2 - x1) / NCURVE + y1 = cveval (cv, x1) + call gamove (gio, x1, y1) + do i = 1, NCURVE { + x1 = x1 + x2 + y1 = cveval (cv, x1) + call gadraw (gio, x1, y1) + } + } +end 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 diff --git a/noao/onedspec/sensfunc/sfmarks.x b/noao/onedspec/sensfunc/sfmarks.x new file mode 100644 index 00000000..39d85af6 --- /dev/null +++ b/noao/onedspec/sensfunc/sfmarks.x @@ -0,0 +1,46 @@ +include <gset.h> +include "sensfunc.h" + +define GMTYPES "|point|box|plus|cross|diamond|hline|vline|hebar|vebar|circle|" + + +# SF_MARKS -- Decode user mark types into GIO mark types. The input string +# consists of two whitespace separated mark types. + +procedure sf_marks (gp, marks) + +pointer gp +char marks[ARB] + +int i, nscan(), strdic() +pointer sp, str + +int gmtypes[10] +data gmtypes /GM_POINT,GM_BOX,GM_PLUS,GM_CROSS,GM_DIAMOND,GM_HLINE,GM_VLINE, + GM_HEBAR,GM_VEBAR,GM_CIRCLE/ + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call sscan (marks) + call gargwrd (Memc[str], SZ_LINE) + if (nscan() == 1) { + i = strdic (Memc[str], Memc[str], SZ_LINE, GMTYPES) + if (i != 0) + GP_MARK(gp) = gmtypes[i] + } + call gargwrd (Memc[str], SZ_LINE) + if (nscan() == 2) { + i = strdic (Memc[str], Memc[str], SZ_LINE, GMTYPES) + if (i != 0) + GP_MDEL(gp) = gmtypes[i] + } + call gargwrd (Memc[str], SZ_LINE) + if (nscan() == 3) { + i = strdic (Memc[str], Memc[str], SZ_LINE, GMTYPES) + if (i != 0) + GP_MADD(gp) = gmtypes[i] + } + + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sfmove.x b/noao/onedspec/sensfunc/sfmove.x new file mode 100644 index 00000000..4451938e --- /dev/null +++ b/noao/onedspec/sensfunc/sfmove.x @@ -0,0 +1,166 @@ +include <gset.h> +include "sensfunc.h" + + +# SF_MOVE -- Move point, star, or wavelength. + +procedure sf_move (gp, stds, nstds, key, istd, ipt, shift) + +pointer gp # GIO pointer +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +int key # Delete point, star, or wavelength +int istd # Index of standard star +int ipt # Index of point +real shift + +int i, j, n, wcs, mark, mdel, cdel, color, stridx() +real wave, szmark, szmdel +pointer x, y, z, w, gio + +begin + gio = GP_GIO(gp) + mdel = GP_MDEL(gp) + cdel = GP_CDEL(gp) + szmdel = GP_SZMDEL(gp) + szmark = GP_SZMARK(gp) + + # Move points in each displayed graph. + for (wcs = 1; GP_GRAPHS(gp,wcs) != EOS; wcs = wcs + 1) { + if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0) + next + + call gseti (gio, G_WCS, wcs) + call sf_data (stds, nstds, GP_GRAPHS(gp,wcs)) + switch (key) { + case 'p': + if (istd != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + x = STD_X(stds[istd])+ipt-1 + y = STD_Y(stds[istd],1)+ipt-1 + w = STD_WTS(stds[istd])+ipt-1 + if (Memr[w] != 0.) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x], Memr[y]+shift, mark, szmark, + szmark) + } else { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x], Memr[y]+shift, mdel, szmdel, + szmdel) + } + case 's': + if (istd != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + n = STD_NWAVES(stds[istd]) + x = STD_X(stds[istd]) + y = STD_Y(stds[istd]) + w = STD_WTS(stds[istd]) + do i = 1, n { + if (Memr[w] != 0.) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x], Memr[y]+shift, mark, szmark, + szmark) + } else { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x], Memr[y]+shift, mdel, szmdel, + szmdel) + } + x = x + 1 + y = y + 1 + w = w + 1 + } + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + if (i != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + n = STD_NWAVES(stds[i]) + x = STD_X(stds[i]) + y = STD_Y(stds[i]) + z = STD_WAVES(stds[i]) + w = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[z] == wave) { + if (Memr[w] != 0.) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mark, szmark, + szmark) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x], Memr[y]+shift, mark, + szmark, szmark) + } else { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, + szmdel) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, cdel) + call gmark (gio, Memr[x], Memr[y]+shift, mdel, + szmdel, szmdel) + } + } + x = x + 1 + y = y + 1 + z = z + 1 + w = w + 1 + } + } + } + } + + # Now add the shift to the data. + switch (key) { + case 'p': + y = STD_SENS(stds[istd])+ipt-1 + Memr[y] = Memr[y] + shift + case 's': + n = STD_NWAVES(stds[istd]) + y = STD_SENS(stds[istd]) + call aaddkr (Memr[y], shift, Memr[y], n) + STD_SHIFT(stds[istd]) = STD_SHIFT(stds[istd]) + shift + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + z = STD_WAVES(stds[i]) + y = STD_SENS(stds[i]) + do j = 1, n { + if (Memr[z] == wave) + Memr[y] = Memr[y] + shift + w = w + 1 + y = y + 1 + } + } + } +end diff --git a/noao/onedspec/sensfunc/sfnearest.x b/noao/onedspec/sensfunc/sfnearest.x new file mode 100644 index 00000000..540faad2 --- /dev/null +++ b/noao/onedspec/sensfunc/sfnearest.x @@ -0,0 +1,69 @@ +include <gset.h> +include <mach.h> +include "sensfunc.h" + +# SF_NEAREST -- Find the nearest point to the cursor. Return the standard +# star index and the wavelength point index. The metric is in NDC. +# The cursor is moved to the nearest point selected. Return zero for +# the standard star index if valid point not found. + +procedure sf_nearest (gp, stds, nstds, wx, wy, wcs, type, istd, ipt) + +pointer gp # Graphics pointer +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +real wx, wy # Cursor position +int wcs # WCS +int type # Type of points (0=not del, 1=del, 2=both) +int istd # Index of standard star (returned) +int ipt # Index of point (returned) + +int i, j, n, stridx() +real x0, y0, x1, y1, x2, y2, r2, r2min +pointer x, y, w, gio + +begin + # Check for valid wc. + istd = 0 + if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0) + return + + # Transform world cursor coordinates to NDC. + gio = GP_GIO(gp) + call gctran (gio, wx, wy, wx, wy, wcs, 0) + + # Search for nearest point. + r2min = MAX_REAL + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + x = STD_X(stds[i]) - 1 + y = STD_Y(stds[i]) - 1 + w = STD_WTS(stds[i]) - 1 + do j = 1, n { + if (type == 0) { + if (Memr[w+j] == 0.) + next + } else if (type == 1) { + if (Memr[w+j] != 0.) + next + } + x1 = Memr[x+j] + y1 = Memr[y+j] + call gctran (gio, x1, y1, x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + istd = i + ipt = j + x2 = x1 + y2 = y1 + } + } + } + + # Move the cursor to the selected point. + call gseti (gio, G_WCS, wcs) + call gscur (gio, x2, y2) +end diff --git a/noao/onedspec/sensfunc/sfoutput.x b/noao/onedspec/sensfunc/sfoutput.x new file mode 100644 index 00000000..e21df280 --- /dev/null +++ b/noao/onedspec/sensfunc/sfoutput.x @@ -0,0 +1,114 @@ +include <mach.h> +include <imhdr.h> +include "sensfunc.h" + + +# SF_OUTPUT -- Write the sensitivity function image. + +procedure sf_output (stds, nstds, cv, output, ignoreaps) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +pointer cv # Sensitivity function curve +char output[SZ_FNAME] # Output root image name (must be SZ_FNAME) +bool ignoreaps # Ignore apertures? + +int i, ap, nw, scan(), nscan() +real w1, w2, dw, dw1, aplow[2], aphigh[2], cveval() +pointer sp, fname, std, im, mw, buf, immap(), mw_open(), impl1r() +errchk imaddi, imaddr + +define makeim_ 99 + +begin + # Return if no output root sensitivity imagename is specified. + if (output[1] == EOS) + return + + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + + # Determine wavelength range and reference standard star. + w1 = MAX_REAL + w2 = -MAX_REAL + dw = MAX_REAL + do i = 1, nstds-2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + std = stds[i] + dw1 = abs ((STD_WEND(std) - STD_WSTART(std)) / (STD_NPTS(std) - 1)) + w1 = min (w1, STD_WSTART(std), STD_WEND(std)) + w2 = max (w2, STD_WSTART(std), STD_WEND(std)) + dw = min (dw, dw1) + } + nw = (w2 - w1) / dw + 1.5 + + # Make output image name with aperture number appended. If the + # image exists allow the user to change root name. +makeim_ + if (ignoreaps) { + call strcpy (output, Memc[fname], SZ_FNAME) + } else { + call sprintf (Memc[fname], SZ_FNAME, "%s.%04d") + call pargstr (output) + call pargi (STD_BEAM(std)) + } + + iferr (im = immap (Memc[fname], NEW_IMAGE, 0)) { + call printf ("Cannot create %s -- Enter new name: ") + call pargstr (Memc[fname]) + call flush (STDOUT) + if (scan() != EOF) { + call gargwrd (Memc[fname], SZ_FNAME) + if (nscan() == 1) { + call strcpy (Memc[fname], output, SZ_FNAME) + goto makeim_ + } + } + call printf ("No sensitivity function created for aperture %2d\n") + call pargi (STD_BEAM(std)) + call flush (STDOUT) + return + } + + # Define the image header. + IM_NDIM(im) = 1 + IM_LEN(im,1) = nw + IM_PIXTYPE(im) = TY_REAL + if (ignoreaps) { + call sprintf (IM_TITLE(im), SZ_FNAME, + "Sensitivity function for all apertures") + } else { + call sprintf (IM_TITLE(im), SZ_FNAME, + "Sensitivity function for aperture %d") + call pargi (STD_BEAM(std)) + } + + mw = mw_open (NULL, 1) + call mw_newsystem (mw, "equispec", 1) + call mw_swtype (mw, 1, 1, "linear", "label=Wavelength units=Angstroms") + call smw_open (mw, NULL, im) + ap = STD_BEAM(std) + aplow[1] = INDEF + aphigh[1] = INDEF + aplow[2] = INDEF + aphigh[2] = INDEF + call smw_swattrs (mw, 1, 1, ap, STD_BEAM(std), 0, + double(w1), double(dw), nw, 0D0, aplow, aphigh, "") + call smw_saveim (mw, im) + call smw_close (mw) + + # Write sensitivity data. + buf = impl1r (im) + do i = 0, nw-1 + Memr[buf+i] = cveval (cv, w1 + i * dw) + + # Notify user. + call printf ("%s --> %s\n") + call pargstr (IM_TITLE(im)) + call pargstr (Memc[fname]) + call flush (STDOUT) + + call imunmap (im) + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sfreset.x b/noao/onedspec/sensfunc/sfreset.x new file mode 100644 index 00000000..fc4d974e --- /dev/null +++ b/noao/onedspec/sensfunc/sfreset.x @@ -0,0 +1,62 @@ +include "sensfunc.h" + +# SF_RESET -- Reset the standard star data to the original input. +# This is called cancel changes and start over. + +procedure sf_reset (stds, nstds, wextn, extn, nextn, ecv, shift) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +real wextn[nextn] # Extinction table wavelengths +real extn[nextn] # Extinction table values +int nextn # Number of extinction values +pointer ecv # Residual extinction curve +int shift # Shift flag + +int i, j, n, err +real exptime, airmass, ext +pointer waves, fluxes, dwaves, counts, sens, iwts, wts + +begin + # Reset the flags, sensitivity, and weight values. + shift = NO + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + STD_FLAG(stds[i]) = SF_INCLUDE + STD_SHIFT(stds[i]) = 0. + n = STD_NWAVES(stds[i]) + exptime = STD_EXPTIME(stds[i]) + airmass = STD_AIRMASS(stds[i]) + waves = STD_WAVES(stds[i]) + fluxes = STD_FLUXES(stds[i]) + dwaves = STD_DWAVES(stds[i]) + counts = STD_COUNTS(stds[i]) + sens = STD_SENS(stds[i]) + iwts = STD_IWTS(stds[i]) + wts = STD_WTS(stds[i]) + do j = 1, n { + call intrp (1, wextn, extn, nextn, Memr[waves], ext, err) + Memr[sens] = Memr[counts] / + (Memr[fluxes] * Memr[dwaves] * exptime) + Memr[sens] = 2.5 * log10 (Memr[sens]) + airmass * ext + Memr[wts] = Memr[iwts] + + waves = waves + 1 + fluxes = fluxes + 1 + dwaves = dwaves + 1 + counts = counts + 1 + sens = sens + 1 + iwts = iwts + 1 + wts = wts + 1 + } + } + + # Reset the added and composite stars. + STD_NWAVES(stds[nstds-1]) = 0 + STD_FLAG(stds[nstds-1]) = SF_DELETE + STD_FLAG(stds[nstds]) = SF_EXCLUDE + + # Clear the residual extinction curve. + call cvfree (ecv) +end diff --git a/noao/onedspec/sensfunc/sfrms.x b/noao/onedspec/sensfunc/sfrms.x new file mode 100644 index 00000000..72b8ea98 --- /dev/null +++ b/noao/onedspec/sensfunc/sfrms.x @@ -0,0 +1,43 @@ +include "sensfunc.h" + + +# SF_RMS -- Compute the RMS of the sensitivity function fit. + +procedure sf_rms (stds, nstds, rms, npts) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +real rms # RMS about fit (returned) +int npts # Number of points in fit (excluding zero wts.) + +int i, j, f, n +pointer x, y, w + +begin + npts = 0 + rms = 0. + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + x = STD_WAVES(stds[i]) + y = STD_SENS(stds[i]) + f = STD_FIT(stds[i]) + w = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[w] != 0.) { + rms = rms + (Memr[y] - Memr[f]) ** 2 + npts = npts + 1 + } + x = x + 1 + y = y + 1 + f = f + 1 + w = w + 1 + } + } + + if (npts > 1) + rms = sqrt (rms / (npts - 1)) + else + rms = INDEF +end diff --git a/noao/onedspec/sensfunc/sfsensfunc.x b/noao/onedspec/sensfunc/sfsensfunc.x new file mode 100644 index 00000000..ee2f1b2a --- /dev/null +++ b/noao/onedspec/sensfunc/sfsensfunc.x @@ -0,0 +1,255 @@ +include <error.h> +include <gset.h> +include <mach.h> +include "sensfunc.h" + +define KEY "noao$onedspec/sensfunc/sensfunc.key" +define PROMPT "sensfunc options" + + +# SF_SENSFUNC -- Interactive sensitivity function determination. + +procedure sf_sensfunc (gp, stds, nstds, wextn, extn, nextn, sensimage, logfile, + ecv, function, order, ignoreaps, interactive) + +pointer gp # Graphics structure +pointer stds[nstds] # Pointer to standard observations +int nstds # Number of standards +real wextn[nextn] # Extinction table wavelengths +real extn[nextn] # Extinction table values +int nextn # Number of extinction table values +char sensimage[ARB] # Output rootname +char logfile[ARB] # Statistics filename +pointer ecv # Residual extinction curve +char function[ARB] # Fitting function type +int order # Function order +bool ignoreaps # Ignore apertures? +int interactive # Interactive? + +char cmd[SZ_FNAME] +int wc, key, newgraph, newfit +real wx, wy + +int i, j, aperture, shift, npts, fd, open() +real xmin, xmax, rms, delta +pointer cv + +int clgcur(), scan(), nscan(), clgwrd() +errchk open + +define output_ 99 + +begin + # Initialize data and do the initial fit. + call sf_reset (stds, nstds, wextn, extn, nextn, ecv, shift) + + xmin = MAX_REAL + xmax = -MAX_REAL + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + aperture = STD_BEAM(stds[i]) + xmin = min (xmin, STD_WSTART(stds[i]), STD_WEND(stds[i])) + xmax = max (xmax, STD_WSTART(stds[i]), STD_WEND(stds[i])) + } + cv = NULL + call sf_fit (stds, nstds, cv, function, order, xmin, xmax) + call sf_rms (stds, nstds, rms, npts) + + # If not interactive go to the output. + if (interactive == 3) + goto output_ + if (interactive != 4) { + call printf ("Fit aperture %d interactively? ") + call pargi (aperture) + interactive = clgwrd ("answer", cmd, SZ_FNAME, "|no|yes|NO|YES") + switch (interactive) { + case 1: + goto output_ + case 3: + call sf_gfree (gp) + goto output_ + } + } + + # Initialize graphics structure parameters: airmass and wavelength + # limits and default images to plot. + + if (gp == NULL) + call sf_ginit (gp) + GP_AIRMASS(gp,1) = MAX_REAL + GP_AIRMASS(gp,2) = -MAX_REAL + j = 0 + do i = 1, nstds - 2 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + GP_AIRMASS(gp,1) = min (GP_AIRMASS(gp,1), STD_AIRMASS(stds[i])) + GP_AIRMASS(gp,2) = max (GP_AIRMASS(gp,2), STD_AIRMASS(stds[i])) + if (j < SF_NGRAPHS) { + j = j + 1 + call strcpy (STD_IMAGE(stds[i]), Memc[GP_IMAGES(gp,j)], + SZ_FNAME) + call strcpy (STD_SKY(stds[i]), Memc[GP_SKYS(gp,j)], SZ_FNAME) + } + } + delta = GP_AIRMASS(gp,2) - GP_AIRMASS(gp,1) + GP_AIRMASS(gp,1) = GP_AIRMASS(gp,1) - 0.05 * delta + GP_AIRMASS(gp,2) = GP_AIRMASS(gp,2) + 0.05 * delta + GP_WSTART(gp) = xmin + GP_WEND(gp) = xmax + call sf_title (gp, aperture, function, order, npts, rms) + + # Enter cursor loop by drawing the graphs. + key = 'r' + repeat { + switch (key) { + case '?': + call gpagefile (GP_GIO(gp), KEY, PROMPT) + case ':': + call sf_colon (cmd, gp, stds, nstds, cv, wextn, extn, nextn, + ecv, function, order, npts, rms, newfit, newgraph) + case 'a': + call sf_add (gp, stds, nstds, cv, wx, wy, wc) + case 'c': + call sf_composite (stds, nstds, cv) + newfit = YES + newgraph = YES + case 'd': + call sf_data (stds, nstds, GP_GRAPHS(gp,wc)) + call sf_nearest (gp, stds, nstds, wx, wy, wc, 0, i, j) + if (i > 0) { + call printf ( + "%s - Delete p(oint), s(tar), or w(avelength):") + call pargstr (STD_IMAGE(stds[i])) + if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF) + break + call printf ("\n") + call sf_delete (gp, stds, nstds, key, i, j) + } + case 'e': + call sf_extinct (gp, stds, nstds, cv, ecv, function, order) + newfit = YES + newgraph = YES + case 'f': + newfit = YES + case 'g': + newgraph = YES + newfit = YES + case 'i': + call sf_data (stds, nstds, GP_GRAPHS(gp,wc)) + call sf_nearest (gp, stds, nstds, wx, wy, wc, 2, i, j) + if (i > 0) { + call printf ( + "%s: airmass=%6.3f wavelen=%6.3f sens=%6.3f fit=%6.3f weight=%3f") + call pargstr (STD_IMAGE(stds[i])) + call pargr (STD_AIRMASS(stds[i])) + call pargr (Memr[STD_WAVES(stds[i])+j-1]) + call pargr (Memr[STD_SENS(stds[i])+j-1]) + call pargr (Memr[STD_FIT(stds[i])+j-1]) + call pargr (Memr[STD_WTS(stds[i])+j-1]) + } + case 'm': + call sf_data (stds, nstds, GP_GRAPHS(gp,wc)) + call sf_nearest (gp, stds, nstds, wx, wy, wc, 2, i, j) + if (i > 0) { + call printf ( + "%s - Move p(oint), s(tar), or w(avelength) to cursor:") + call pargstr (STD_IMAGE(stds[i])) + if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF) + break + call printf ("\n") + delta = wy - Memr[STD_Y(stds[i])+j-1] + call sf_move (gp, stds, nstds, key, i, j, delta) + } + case 'o': + call sf_reset (stds, nstds, wextn, extn, nextn, ecv, shift) + newfit = YES + newgraph = YES + case 'q': + break + case 'I': + call fatal (0, "Interrupt") + case 'r': + newgraph = YES + case 's': + call sf_shift (stds, nstds, shift) + newfit=YES + newgraph=YES + case 'u': + call sf_data (stds, nstds, GP_GRAPHS(gp,wc)) + call sf_nearest (gp, stds, nstds, wx, wy, wc, 1, i, j) + if (i > 0) { + call printf ( + "%s - Undelete p(oint), s(tar), or w(avelength):") + call pargstr (STD_IMAGE(stds[i])) + if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF) + break + call printf ("\n") + call sf_undelete (gp, stds, nstds, key, i, j) + } + case 'w': + call sf_data (stds, nstds, GP_GRAPHS(gp,wc)) + call sf_nearest (gp, stds, nstds, wx, wy, wc, 0, i, j) + if (i > 0) { + call printf ( + "%s - Reweight p(oint), s(tar), or w(avelength):") + call pargstr (STD_IMAGE(stds[i])) + if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF) + break + call printf ("New weight (%g):") + call pargr (Memr[STD_IWTS(stds[i])+j-1]) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (delta) + if (nscan() == 1) + call sf_weights (stds, nstds, key, i, j, delta) + } + call printf ("\n") + } + default: + call printf ("\007") + } + + # Do a new fit and recompute the RMS, and title string. + if (newfit == YES) { + call sf_fit (stds, nstds, cv, function, order, xmin, xmax) + call sf_rms (stds, nstds, rms, npts) + call sf_title (gp, aperture, function, order, npts, rms) + do i = 1, SF_NGRAPHS + if (GP_SHDR(gp,i) != NULL) + call shdr_close (GP_SHDR(gp,i)) + } + + # Draw new graphs. + if (newgraph == YES) { + call sf_graph (gp, stds, nstds, cv, wextn, extn, nextn, ecv) + newgraph = NO + newfit = YES + } + + # Overplot new fit. + if (newfit == YES) { + call sf_fitgraph (gp, cv) + newfit = NO + } + } until (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME) == EOF) + + # Close any open images. + do i = 1, SF_NGRAPHS + if (GP_SHDR(gp,i) != NULL) + call shdr_close (GP_SHDR(gp,i)) + +output_ + # Output the sensitivity function and logfile statistics. + call sf_output (stds, nstds, cv, sensimage, ignoreaps) + if (logfile[1] != EOS) { + iferr { + fd = open (logfile, APPEND, TEXT_FILE) + call sf_stats (fd, stds, nstds, function, order, npts, rms) + call sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn, ecv) + call close (fd) + } then + call erract (EA_WARN) + } + call cvfree (cv) +end diff --git a/noao/onedspec/sensfunc/sfshift.x b/noao/onedspec/sensfunc/sfshift.x new file mode 100644 index 00000000..07b204f3 --- /dev/null +++ b/noao/onedspec/sensfunc/sfshift.x @@ -0,0 +1,81 @@ +include "sensfunc.h" + + +# SF_SHIFT -- Shift or unshift all standard stars to have zero mean residual. + +procedure sf_shift (stds, nstds, flag) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +int flag # Shift flag + +pointer x, y, w, f +int i, j, n, nshift +real shift, shift1, minshift + +begin + # If flag is YES then unshift the data. + if (flag == YES) { + do i = 1, nstds { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + n = STD_NWAVES(stds[i]) + if (n == 0) + next + y = STD_SENS(stds[i]) + call asubkr (Memr[y], STD_SHIFT(stds[i]), Memr[y], n) + STD_SHIFT(stds[i]) = 0. + } + flag = NO + call printf ("Data unshifted") + return + } + + # Determine the shifts needed to make the mean residual zero. + # Also determine the minimum shift. + + minshift = 0. + do i = 1, nstds { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + n = STD_NWAVES(stds[i]) + if (n == 0) + next + x = STD_WAVES(stds[i]) + y = STD_SENS(stds[i]) + w = STD_WTS(stds[i]) + f = STD_FIT(stds[i]) + nshift = 0 + shift = 0. + shift1 = 0. + do j = 1, n { + shift1 = shift1 + Memr[f+j-1] - Memr[y+j-1] + if (Memr[w+j-1] > 0.) { + shift = shift + Memr[f+j-1] - Memr[y+j-1] + nshift = nshift + 1 + } + } + if (nshift > 0) { + shift = STD_SHIFT(stds[i]) + shift / nshift + if (shift < minshift) + minshift = shift + } else + shift = STD_SHIFT(stds[i]) + shift1 / n + STD_SHIFT(stds[i]) = shift + } + + # Adjust the shifts to be upwards. + do i = 1, nstds { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + n = STD_NWAVES(stds[i]) + if (n == 0) + next + y = STD_SENS(stds[i]) + shift = STD_SHIFT(stds[i]) - minshift + call aaddkr (Memr[y], shift, Memr[y], n) + STD_SHIFT(stds[i]) = shift + } + flag = YES + call printf ("Data shifted") +end diff --git a/noao/onedspec/sensfunc/sfstats.x b/noao/onedspec/sensfunc/sfstats.x new file mode 100644 index 00000000..a94691a4 --- /dev/null +++ b/noao/onedspec/sensfunc/sfstats.x @@ -0,0 +1,152 @@ +include "sensfunc.h" + + +# SF_STATS -- Print basic statistics about the stars and the fit. + +procedure sf_stats (fd, stds, nstds, function, order, npts, rms) + +int fd # Output file descriptor (may be STDOUT) +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +char function[ARB] # Fitted function +int order # Order of function +int npts # Number of points in fit +real rms # RMS of fit + +int i, j, n +real rms1, dev1, dev2, dev3 +pointer sp, str, wts + +begin + # Start with system ID. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call sysid (Memc[str], SZ_LINE) + + # Determine beam from first standard star not excluded. + for (i=1; (i<nstds) && (STD_FLAG(stds[i])==SF_EXCLUDE); i=i+1) + ; + call fprintf (fd, "%s\n") + call pargstr (Memc[str]) + call fprintf (fd, "Sensitivity function for aperture %d:\n") + call pargi (STD_BEAM(stds[i])) + call fprintf (fd, + "Fitting function is %s of order %d with %d points and RMS of %6.4f.\n") + call pargstr (function) + call pargi (order) + call pargi (npts) + call pargr (rms) + + call fprintf (fd, "%12s %7s %7s %7s %7s %7s %7s %7s\n") + call pargstr ("Image") + call pargstr ("Airmass") + call pargstr ("Points") + call pargstr ("Shift") + call pargstr ("RMS Fit") + call pargstr ("Dev 1") + call pargstr ("Dev 2") + call pargstr ("Dev 3") + + do i = 1, nstds { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + + n = 0 + wts = STD_WTS(stds[i]) - 1 + for (j=1; j<=STD_NWAVES(stds[i]); j=j+1) + if (Memr[wts+j] != 0.) + n = n + 1 + if ((i == nstds-1) && (n == 0)) + next + + call sf_devs (stds[i], rms1, dev1, dev2, dev3) + + call fprintf (fd, "%12s %7.3f %7d %7.4f %7.4f %7.4f %7.4f %7.4f") + call pargstr (STD_IMAGE(stds[i])) + call pargr (STD_AIRMASS(stds[i])) + call pargi (n) + call pargr (STD_SHIFT(stds[i])) + call pargr (rms1) + call pargr (dev1) + call pargr (dev2) + call pargr (dev3) + + if (n == 0) { + call fprintf (fd, "%s") + call pargstr (" <-- deleted") + } + call fprintf (fd, "\n") + } + + # Trailing spacer + call fprintf (fd, "\n") +end + + +# SF_DEVS - Compute rms and mean deviations from the fit. +# The deviations are computed in three segments. + +procedure sf_devs (std, rms, dev1, dev2, dev3) + +pointer std # Standard star data +real rms # RMS about fit +real dev1 # Average deviation in first third of data +real dev2 # Average deviation in second third of data +real dev3 # Average deviation in last third of data + +int i, ndev1, ndev2, ndev3, nrms, nbin, nwaves +real dev +pointer sens, fit, wts + +begin + # Get elements froms standard star structure. + nwaves = STD_NWAVES(std) + sens = STD_SENS(std) + fit = STD_FIT(std) + wts = STD_WTS(std) + + # Divide into thirds. + rms = 0. + ndev1 = 0 + dev1 = 0. + nbin = nwaves / 3 + for (i=1; i<= nbin; i=i+1) + if (Memr[wts+i-1] != 0.) { + dev = Memr[sens+i-1] - Memr[fit+i-1] + dev1 = dev1 + dev + rms = rms + dev ** 2 + ndev1 = ndev1 + 1 + } + if (ndev1 > 0) + dev1 = dev1 / ndev1 + + ndev2 = 0 + dev2 = 0. + nbin = 2 * nwaves / 3 + for (; i<=nbin; i=i+1) + if (Memr[wts+i-1] != 0.) { + dev = Memr[sens+i-1] - Memr[fit+i-1] + dev2 = dev2 + dev + rms = rms + dev ** 2 + ndev2 = ndev2 + 1 + } + if (ndev2 > 0) + dev2 = dev2 / ndev2 + + ndev3 = 0 + dev3 = 0. + nbin = nwaves + for (; i<=nbin; i=i+1) + if (Memr[wts+i-1] != 0.) { + dev = Memr[sens+i-1] - Memr[fit+i-1] + dev3 = dev3 + dev + rms = rms + dev ** 2 + ndev3 = ndev3 + 1 + } + if (ndev3 > 0) + dev3 = dev3 / ndev3 + + nrms = ndev1 + ndev2 + ndev3 + if (nrms > 0) + rms = sqrt (rms / nrms) +end diff --git a/noao/onedspec/sensfunc/sfstds.x b/noao/onedspec/sensfunc/sfstds.x new file mode 100644 index 00000000..07219729 --- /dev/null +++ b/noao/onedspec/sensfunc/sfstds.x @@ -0,0 +1,266 @@ +include "sensfunc.h" + + +# SF_STDS -- Get the standard observations for the specified apertures. +# If ignoring aperture set all apertures to 1. +# This routine knows the output of the task STANDARD. + +procedure sf_stds (standards, aps, ignoreaps, stds, nstds) + +char standards # Standard star data file +pointer aps # Aperture list +bool ignoreaps # Ignore apertures? +pointer stds # Pointer to standard observations (returned) +int nstds # Number of standard observations (returned) + +int i, j, fd, beam, npts, nwaves, nalloc +real exptime, airmass, wstart, wend +real wavelength, flux, dwave, count +pointer sp, image, title, std +pointer waves, fluxes, dwaves, counts, sens, fit, wts, iwts, x, y + +bool rng_elementi() +int open(), fscan(), nscan(), stridxs() +errchk open, malloc, realloc + +begin + call smark (sp) + call salloc (image, SZ_STDIMAGE, TY_CHAR) + call salloc (title, SZ_STDTITLE, TY_CHAR) + + # Open the standard observation data file. + fd = open (standards, READ_ONLY, TEXT_FILE) + + # Read the standard observations and create a structure for each one. + # The beginning of a new star is found by a line whose first word + # begins with the character '['. Otherwise the line is interpreted + # as a data line. All unrecognized formats are skipped. + + nwaves = 0 + nstds = 0 + while (fscan (fd) != EOF) { + call gargwrd (Memc[image], SZ_STDIMAGE) + if (Memc[image] == '[') { + call gargi (beam) + call gargi (npts) + call gargr (exptime) + call gargr (airmass) + call gargr (wstart) + call gargr (wend) + call gargstr (Memc[title], SZ_STDTITLE) + if (nscan() < 7) + next + if (!rng_elementi (aps, beam)) + next + if (IS_INDEF (exptime) || exptime <= 0.) { + call eprintf ( + "%s: Warning - exposure time missing or zero, using 1 second\n") + call pargstr (Memc[image]) + exptime = 1. + } + + # For the first one create the pointer to the array of + # structures. For the following stars increase the size + # of the pointer array and finish up the previous standard + # star. + + if (nstds == 0) { + nstds = nstds + 1 + call calloc (stds, nstds, TY_INT) + call calloc (std, LEN_STD, TY_STRUCT) + Memi[stds+nstds-1] = std + } else { + if (nwaves > 0) { + call realloc (waves, nwaves, TY_REAL) + call realloc (fluxes, nwaves, TY_REAL) + call realloc (dwaves, nwaves, TY_REAL) + call realloc (counts, nwaves, TY_REAL) + call realloc (wts, nwaves, TY_REAL) + call malloc (sens, nwaves, TY_REAL) + call malloc (fit, nwaves, TY_REAL) + call malloc (iwts, nwaves, TY_REAL) + call malloc (x, nwaves, TY_REAL) + call malloc (y, nwaves, TY_REAL) + call amovr (Memr[wts], Memr[iwts], nwaves) + STD_NWAVES(std) = nwaves + STD_WAVES(std) = waves + STD_FLUXES(std) = fluxes + STD_DWAVES(std) = dwaves + STD_COUNTS(std) = counts + STD_SENS(std) = sens + STD_FIT(std) = fit + STD_WTS(std) = wts + STD_IWTS(std) = iwts + STD_X(std) = x + STD_Y(std) = y + + nstds = nstds + 1 + call realloc (stds, nstds, TY_INT) + call calloc (std, LEN_STD, TY_STRUCT) + Memi[stds+nstds-1] = std + } + } + + # Start a new standard star. + std = Memi[stds+nstds-1] + if (ignoreaps) + STD_BEAM(std) = 1 + else + STD_BEAM(std) = beam + STD_NPTS(std) = npts + STD_EXPTIME(std) = exptime + STD_AIRMASS(std) = airmass + STD_WSTART(std) = wstart + STD_WEND(std) = wend + STD_SHIFT(std) = 0. + STD_NWAVES(std) = 0 + + # Decode the image and sky strings. + call strcpy (Memc[title], STD_TITLE(std), SZ_STDTITLE) + i = stridxs ("]", Memc[image]) + if (Memc[image+i] == ']') + i = i + 1 + Memc[image+i-1] = EOS + call strcpy (Memc[image+1], STD_IMAGE(std), SZ_STDIMAGE) + if (Memc[image+i] == '-') { + i = i + 2 + j = stridxs ("]", Memc[image+i]) + i + Memc[image+j-1] = EOS + call strcpy (Memc[image+i], STD_SKY(std), SZ_STDIMAGE) + } else + STD_SKY(std) = EOS + nwaves = 0 + + # Interprete the line as standard star wavelength point. + } else if (nstds > 0) { + call reset_scan() + call gargr (wavelength) + call gargr (flux) + call gargr (dwave) + call gargr (count) + if (nscan() < 3) + next + if (wavelength < min (wstart, wend) || + wavelength > max (wstart, wend) || + flux<=0. || dwave<=0. || count<=0.) + next + if (!rng_elementi (aps, beam)) + next + nwaves = nwaves + 1 + + # Allocate in blocks to minimize the number of reallocations. + if (nwaves == 1) { + nalloc = 100 + call malloc (waves, nalloc, TY_REAL) + call malloc (fluxes, nalloc, TY_REAL) + call malloc (dwaves, nalloc, TY_REAL) + call malloc (counts, nalloc, TY_REAL) + call malloc (wts, nalloc, TY_REAL) + } else if (nwaves > nalloc) { + nalloc = nalloc + 100 + call realloc (waves, nalloc, TY_REAL) + call realloc (fluxes, nalloc, TY_REAL) + call realloc (dwaves, nalloc, TY_REAL) + call realloc (counts, nalloc, TY_REAL) + call realloc (wts, nalloc, TY_REAL) + } + + # Record the data and compute the sensitivity. + Memr[waves+nwaves-1] = wavelength + Memr[fluxes+nwaves-1] = flux + Memr[dwaves+nwaves-1] = dwave + Memr[counts+nwaves-1] = count + Memr[wts+nwaves-1] = 1. + } + } + + # Finish up the last standard star and close the file. + if (nstds > 0) { + STD_NWAVES(std) = nwaves + if (nwaves > 0) { + call realloc (waves, nwaves, TY_REAL) + call realloc (fluxes, nwaves, TY_REAL) + call realloc (dwaves, nwaves, TY_REAL) + call realloc (counts, nwaves, TY_REAL) + call realloc (wts, nwaves, TY_REAL) + call malloc (sens, nwaves, TY_REAL) + call malloc (fit, nwaves, TY_REAL) + call malloc (iwts, nwaves, TY_REAL) + call malloc (x, nwaves, TY_REAL) + call malloc (y, nwaves, TY_REAL) + call amovr (Memr[wts], Memr[iwts], nwaves) + STD_WAVES(std) = waves + STD_FLUXES(std) = fluxes + STD_DWAVES(std) = dwaves + STD_COUNTS(std) = counts + STD_SENS(std) = sens + STD_FIT(std) = fit + STD_WTS(std) = wts + STD_IWTS(std) = iwts + STD_X(std) = x + STD_Y(std) = y + } + } + call close (fd) + call sfree (sp) + + # Add standard stars for any added points and composite points. + nstds = nstds + 2 + call realloc (stds, nstds, TY_INT) + call calloc (std, LEN_STD, TY_STRUCT) + Memi[stds+nstds-2] = std + call strcpy ("Added", STD_IMAGE(std), SZ_STDIMAGE) + STD_BEAM(std) = STD_BEAM(Memi[stds]) + STD_NPTS(std) = STD_NPTS(Memi[stds]) + STD_EXPTIME(std) = 1. + STD_AIRMASS(std) = 1. + STD_WSTART(std) = STD_WSTART(Memi[stds]) + STD_WEND(std) = STD_WEND(Memi[stds]) + STD_SHIFT(std) = 0. + STD_NWAVES(std) = 0 + call calloc (std, LEN_STD, TY_STRUCT) + Memi[stds+nstds-1] = std + call strcpy ("Composite", STD_IMAGE(std), SZ_STDIMAGE) + STD_BEAM(std) = STD_BEAM(Memi[stds]) + STD_NPTS(std) = STD_NPTS(Memi[stds]) + STD_EXPTIME(std) = 1. + STD_AIRMASS(std) = 1. + STD_WSTART(std) = STD_WSTART(Memi[stds]) + STD_WEND(std) = STD_WEND(Memi[stds]) + STD_SHIFT(std) = 0. + STD_NWAVES(std) = 0 +end + + +# SF_FREE -- Free the standard observations and aperture array. + +procedure sf_free (stds, nstds, apertures, napertures) + +pointer stds # Pointer to standard observations +int nstds # Number of standard observations +pointer apertures # Pointer to apertures array +int napertures # Number of apertures + +int i +pointer std + +begin + do i = 1, nstds { + std = Memi[stds+i-1] + if (STD_NWAVES(std) > 0) { + call mfree (STD_WAVES(std), TY_REAL) + call mfree (STD_FLUXES(std), TY_REAL) + call mfree (STD_DWAVES(std), TY_REAL) + call mfree (STD_COUNTS(std), TY_REAL) + call mfree (STD_SENS(std), TY_REAL) + call mfree (STD_FIT(std), TY_REAL) + call mfree (STD_WTS(std), TY_REAL) + call mfree (STD_IWTS(std), TY_REAL) + call mfree (STD_X(std), TY_REAL) + call mfree (STD_Y(std), TY_REAL) + } + call mfree (std, TY_STRUCT) + } + call mfree (stds, TY_INT) + call mfree (apertures, TY_INT) +end diff --git a/noao/onedspec/sensfunc/sftitle.x b/noao/onedspec/sensfunc/sftitle.x new file mode 100644 index 00000000..50cece9a --- /dev/null +++ b/noao/onedspec/sensfunc/sftitle.x @@ -0,0 +1,23 @@ +include "sensfunc.h" + + +# SF_TITLE -- Make title string for graphs. + +procedure sf_title (gp, aperture, function, order, npts, rms) + +pointer gp +int aperture +char function[ARB] +int order +int npts +real rms + +begin + call sprintf (GP_TITLE(gp), GP_SZTITLE, + "Aperture=%d Function=%s Order=%d Points=%d RMS=%6.4f") + call pargi (aperture) + call pargstr (function) + call pargi (order) + call pargi (npts) + call pargr (rms) +end diff --git a/noao/onedspec/sensfunc/sfundelete.x b/noao/onedspec/sensfunc/sfundelete.x new file mode 100644 index 00000000..25161cc3 --- /dev/null +++ b/noao/onedspec/sensfunc/sfundelete.x @@ -0,0 +1,141 @@ +include <gset.h> +include "sensfunc.h" + + +# SF_UNDELETE -- Unelete point, star, or wavelength. + +procedure sf_undelete (gp, stds, nstds, key, istd, ipt) + +pointer gp # GIO pointer +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +int key # Delete point, star, or wavelength +int istd # Index of standard star +int ipt # Index of point + +int i, j, n, wcs, mark, mdel, color, stridx() +real wave, szmark, szmdel +pointer x, y, z, w, w1, gio + +begin + gio = GP_GIO(gp) + mdel = GP_MDEL(gp) + szmdel = GP_SZMDEL(gp) + szmark = GP_SZMARK(gp) + + # Undelete points from each displayed graph. + for (wcs = 1; GP_GRAPHS(gp,wcs) != EOS; wcs = wcs + 1) { + if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0) + next + + call gseti (gio, G_WCS, wcs) + call gseti (gio, G_PMLTYPE, 0) + call sf_data (stds, nstds, GP_GRAPHS(gp,wcs)) + switch (key) { + case 'p': + if (istd != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + x = STD_X(stds[istd])+ipt-1 + y = STD_Y(stds[istd])+ipt-1 + w = STD_WTS(stds[istd])+ipt-1 + w1 = STD_IWTS(stds[istd])+ipt-1 + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x], Memr[y], mark, szmark , szmark) + case 's': + if (istd != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + n = STD_NWAVES(stds[istd]) + x = STD_X(stds[istd]) + y = STD_Y(stds[istd]) + w = STD_WTS(stds[istd]) + do j = 1, n { + if (Memr[w] == 0.) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark) + } + x = x + 1 + y = y + 1 + w = w + 1 + } + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + if (i != nstds-1) { + mark = GP_MARK(gp) + color = GP_CMARK(gp) + } else { + mark = GP_MADD(gp) + color = GP_CADD(gp) + } + n = STD_NWAVES(stds[i]) + x = STD_X(stds[i]) + y = STD_Y(stds[i]) + z = STD_WAVES(stds[i]) + w = STD_WTS(stds[i]) + do j = 1, n { + if ((Memr[z] == wave) && (Memr[w] == 0.)) { + call gseti (gio, G_PMLTYPE, 0) + call gmark (gio, Memr[x], Memr[y], mdel, szmdel, + szmdel) + call gseti (gio, G_PMLTYPE, 1) + call gseti (gio, G_PLCOLOR, color) + call gmark (gio, Memr[x], Memr[y], mark, szmark, + szmark) + } + x = x + 1 + y = y + 1 + z = z + 1 + w = w + 1 + } + } + } + } + + # Now actually undelete the points by resetting the weights. + switch (key) { + case 'p': + w = STD_WTS(stds[istd])+ipt-1 + w1 = STD_IWTS(stds[istd])+ipt-1 + Memr[w] = Memr[w1] + case 's': + n = STD_NWAVES(stds[istd]) + w = STD_WTS(stds[istd]) + w1 = STD_IWTS(stds[istd]) + call amovr (Memr[w1], Memr[w], n) + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + z = STD_WAVES(stds[i]) + w = STD_WTS(stds[i]) + w1 = STD_IWTS(stds[i]) + do j = 1, n { + if (Memr[z] == wave) + Memr[w] = Memr[w1] + z = z + 1 + w = w + 1 + w1 = w1 + 1 + } + } + } +end diff --git a/noao/onedspec/sensfunc/sfvstats.x b/noao/onedspec/sensfunc/sfvstats.x new file mode 100644 index 00000000..add49da7 --- /dev/null +++ b/noao/onedspec/sensfunc/sfvstats.x @@ -0,0 +1,104 @@ +include "sensfunc.h" + +# SF_VSTATS -- Verbose statistics output. + +procedure sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn, ecv) + +int fd # Output file descriptor (may be STDOUT) +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +pointer cv # Sensitivity function curve +real wextn[nextn] # Extinction table wavelength +real extn[nextn] # Extinction table values +int nextn # Number of extinction table values +pointer ecv # Residual extinction curve + +int i, j, n, nwaves +real w, fit, ext, dext, cveval() +double sum, sum2, s +pointer sp, waves, sens, xp, yp, zp + +begin + nwaves = 0 + do i = 1, nstds-1 + if (STD_FLAG(stds[i]) != SF_EXCLUDE) + nwaves = nwaves + STD_NWAVES(stds[i]) + + call smark (sp) + call salloc (waves, nwaves, TY_REAL) + call salloc (sens, nwaves, TY_REAL) + + nwaves = 0 + do i = 1, nstds-1 { + if (STD_FLAG(stds[i]) == SF_EXCLUDE) + next + n = STD_NWAVES(stds[i]) + xp = STD_WAVES(stds[i]) + yp = STD_SENS(stds[i]) + zp = STD_WTS(stds[i]) + do j = 1, n { + if (Memr[zp] != 0.) { + Memr[waves+nwaves] = Memr[xp] + Memr[sens+nwaves] = Memr[yp] + nwaves = nwaves + 1 + } + xp = xp + 1 + yp = yp + 1 + zp = zp + 1 + } + } + call xt_sort2 (Memr[waves], Memr[sens], nwaves) + + call fprintf (fd, "%8s %7s %7s %7s %7s %5s %7s %7s\n") + call pargstr ("Lambda") + call pargstr ("Fit") + call pargstr ("Avg") + call pargstr ("Resid") + call pargstr ("SD Avg") + call pargstr ("N") + call pargstr ("Ext") + call pargstr ("Dext") + + dext = 0. + n = 0 + sum = 0. + sum2 = 0. + do i = 0, nwaves-1 { + w = Memr[waves+i] + s = Memr[sens+i] + n = n + 1 + sum = sum + s + sum2 = sum2 + s * s + + if ((i < nwaves-1) && (w == Memr[waves+i+1])) + next + + sum = sum / n + sum2 = sum2 / n - sum * sum + if (sum2 > 0) + sum2 = sqrt (sum2 / n) + else + sum2 = 0. + fit = cveval (cv, w) + call intrp (1, wextn, extn, nextn, w, ext, j) + if (ecv != NULL) + dext = cveval (ecv, w) + call fprintf (fd, "%8.2f %7.3f %7.3f %7.4f %7.4f %5d %7.4f %7.4f\n") + call pargr (w) + call pargr (fit) + call pargd (sum) + call pargd (sum - fit) + call pargd (sum2) + call pargi (n) + call pargr (ext) + call pargr (dext) + n = 0 + sum = 0. + sum2 = 0. + } + + # Trailing spacer + call fprintf (fd, "\n") + + call sfree (sp) +end diff --git a/noao/onedspec/sensfunc/sfweights.x b/noao/onedspec/sensfunc/sfweights.x new file mode 100644 index 00000000..2ce24b1a --- /dev/null +++ b/noao/onedspec/sensfunc/sfweights.x @@ -0,0 +1,51 @@ +include "sensfunc.h" + + +# SF_WEIGHTS -- Change weights for point, star, or wavelength. +# The original input weight is permanently lost. + +procedure sf_weights (stds, nstds, key, istd, ipt, weight) + +pointer stds[nstds] # Standard star data +int nstds # Number of standard stars +int key # Delete point, star, or wavelength +int istd # Index of standard star +int ipt # Index of point +real weight # New weight + +int i, j, n +real wave +pointer z, w, iw + +begin + switch (key) { + case 'p': + Memr[STD_WTS(stds[istd])+ipt-1] = weight + Memr[STD_IWTS(stds[istd])+ipt-1] = weight + case 's': + n = STD_NWAVES(stds[istd]) + w = STD_WTS(stds[istd]) + iw = STD_IWTS(stds[istd]) + call amovkr (weight, Memr[w], n) + call amovkr (weight, Memr[iw], n) + case 'w': + wave = Memr[STD_WAVES(stds[istd])+ipt-1] + do i = 1, nstds { + if (STD_FLAG(stds[i]) != SF_INCLUDE) + next + n = STD_NWAVES(stds[i]) + z = STD_WAVES(stds[i]) + w = STD_WTS(stds[i]) + iw = STD_IWTS(stds[i]) + do j = 1, n { + if (Memr[z] == wave) { + Memr[w] = weight + Memr[iw] = weight + } + w = w + 1 + iw = iw + 1 + z = z + 1 + } + } + } +end diff --git a/noao/onedspec/sensfunc/t_sensfunc.x b/noao/onedspec/sensfunc/t_sensfunc.x new file mode 100644 index 00000000..82f18678 --- /dev/null +++ b/noao/onedspec/sensfunc/t_sensfunc.x @@ -0,0 +1,99 @@ +include "sensfunc.h" + + +# T_SENSFUNC -- Determine sensitivities and residual extinctions. +# The input is a file of standard star produced by the task STANDARD. +# The input data is read into an array of structures, one per standard +# star. The stars common to the aperture to be fit are flagged +# and then the data is passed to the main routine SF_SENSFUNC. +# This routine determines a sensitivity curve for the aperture which +# is output by the procedure as well as some optional statistical +# information. It returns an optional residual extinction curve +# for each aperture. The residual extinctions curves are finally combined +# and output as a revised extinction table. + +procedure t_sensfunc () + +pointer standards # Input standard star data filename +pointer sensitivity # Output root sensitivity function imagename +pointer aps # Aperture list +bool ignoreaps # Ignore apertures? +pointer logfile # Output log for statistics +pointer function # Sensitivity function type +int order # Order of sensitivity function +int interactive # Interactive? + +int i, j, aperture, nstds, napertures, nextn, clgeti() +pointer sp, str, stds, apertures, wextn, extn, ecvs, gp +bool clgetb() +pointer rng_open() +errchk sf_sensfunc + +begin + call smark (sp) + call salloc (standards, SZ_FNAME, TY_CHAR) + call salloc (sensitivity, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (logfile, SZ_FNAME, TY_CHAR) + call salloc (function, SZ_FNAME, TY_CHAR) + + # CL parameter input. + call clgstr ("standards", Memc[standards], SZ_FNAME) + call clgstr ("sensitivity", Memc[sensitivity], SZ_FNAME) + call clgstr ("apertures", Memc[str], SZ_LINE) + ignoreaps = clgetb ("ignoreaps") + call clgstr ("logfile", Memc[logfile], SZ_FNAME) + call clgstr ("function", Memc[function], SZ_FNAME) + order = clgeti ("order") + if (clgetb ("interactive")) + interactive = 2 + else + interactive = 3 + + # Decode aperture list. + iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF)) + call error (0, "Bad aperture list") + + # Get the standard star data, the aperture array, and the + # extinction table, and allocate and initialize an array of + # residual extinction curves for each aperture. + + call sf_stds (Memc[standards], aps, ignoreaps, stds, nstds) + if (nstds == 0) { + call sfree (sp) + return + } + call sf_apertures (Memi[stds], nstds, apertures, napertures) + call ext_load (wextn, extn, nextn) + call salloc (ecvs, napertures, TY_INT) + call amovki (NULL, Memi[ecvs], napertures) + + # For each aperture flag standard stars to be used and call sf_sensfunc. + gp = NULL + do j = 1, napertures { + aperture = Memi[apertures+j-1] + do i = 1, nstds - 2 + if (STD_BEAM(Memi[stds+i-1]) == aperture) + STD_FLAG(Memi[stds+i-1]) = SF_INCLUDE + else + STD_FLAG(Memi[stds+i-1]) = SF_EXCLUDE + call sf_sensfunc (gp, Memi[stds], nstds, Memr[wextn], Memr[extn], + nextn, Memc[sensitivity], Memc[logfile], Memi[ecvs+j-1], + Memc[function], order, ignoreaps, interactive) + } + call sf_gfree (gp) + + # Output a revised extinction table by combining the residual + # extinction curves for the apertures. The table name is obtained + # by this proceudre. + + call sf_eout (Memr[wextn], Memr[extn], nextn, Memi[ecvs], napertures) + + # Finish up. + call sf_free (stds, nstds, apertures, napertures) + call ext_free (wextn, extn) + do j = 1, napertures + call cvfree (Memi[ecvs+j-1]) + call rng_close (aps) + call sfree (sp) +end |