diff options
Diffstat (limited to 'noao/onedspec/t_standard.x')
| -rw-r--r-- | noao/onedspec/t_standard.x | 835 |
1 files changed, 835 insertions, 0 deletions
diff --git a/noao/onedspec/t_standard.x b/noao/onedspec/t_standard.x new file mode 100644 index 00000000..9a596150 --- /dev/null +++ b/noao/onedspec/t_standard.x @@ -0,0 +1,835 @@ +include <error.h> +include <imhdr.h> +include <imset.h> +include <gset.h> +include <mach.h> +include <pkg/gtools.h> +include <smw.h> + +define KEY "noao$onedspec/standard.key" +define PROMPT "STANDARD options" + +define VLIGHT 2.997925e18 # Velocity of light in Angstroms/sec +define EXT_LOOKUP 1 # Interp entry ID for extinction table +define MAG_LOOKUP 2 # Interp entry ID for magnitude table + +define STD_LEN 13 # Length of standard structure +define STD_AP Memi[$1] # Aperture number +define STD_TYPE Memi[$1+1] # Spectrum type +define STD_SH Memi[$1+2] # Pointer to spectrum parameters +define STD_IFLAG Memi[$1+3] # Interactive flag +define STD_NWAVES Memi[$1+4] # Number of calibration points +define STD_WAVES Memi[$1+5] # Pointer to standard star wavelengths +define STD_DWAVES Memi[$1+6] # Pointer to standard star bandpasses +define STD_MAGS Memi[$1+7] # Pointer to standard star magnitudes +define STD_FLUXES Memi[$1+8] # Pointer to standard star fluxes +define CAL_NWAVES Memi[$1+9] # Number of calibration points +define CAL_WAVES Memi[$1+10] # Pointer to calibration wavelengths +define CAL_DWAVES Memi[$1+11] # Pointer to calibration bandpasses +define CAL_MAGS Memi[$1+12] # Pointer to calibration magnitudes + +# Object flags +define NONE -1 # No object flag +define SKY 0 # Sky +define OBJ 1 # Object + +# Interactive flags +define ANSWERS "|no|yes|N|Y|NO|YES|NO!|YES!|" +define NO1 1 # No for a single spectrum +define YES1 2 # Yes for a single spectrum +define N2 3 # No for all spectra of the same aperture +define Y2 4 # Yes for all spectra of the same aperture +define NO2 5 # No for all spectra of the same aperture +define YES2 6 # Yes for all spectra of the same aperture +define NO3 7 # No for all spectra +define YES3 8 # Yes for all spectra + +# T_STANDARD -- Read standard star spectrum and compare with tabulated +# fluxes for given star to ascertain the system sensitivity +# across the spectrum. The user may optionally define +# new and arbitrary bandpasses +# +# The sensitivity function is stored in tabular form in a file +# containing the wavelength, sensitivity factor, and counts per +# bandpass at each required position along the spectrum. +# The file is appended to for each new measurement from either +# same or different standard stars. + +procedure t_standard() + +int list # List of input spectra +pointer output # Output standard file +pointer observatory # Observatory +pointer aps # Aperture list +real bandwidth # Width of bandpass +real bandsep # Separation of bandpass +bool bswitch # Beam switch? +bool samestar # Same star in all apertures? +int interactive # Interactive bandpass definition + +bool newobs, obshead +int i, j, line, enwaves, nstds +real wave, dwave, latitude +pointer sp, units, errstr, str, image, ewaves, emags +pointer im, mw, un, unang, sh, obj, sky, std, stds, obs, gp + +int imtgetim(), errget() +real clgetr(), obsgetr() +bool clgetb(), rng_elementi(), streq() +pointer imtopenp(), rng_open(), immap(), smw_openim(), un_open() +errchk immap, smw_openim, shdr_open, std_calib, get_airm, ext_load, obsimopen +errchk un_open, std_gcalib + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (observatory, SZ_FNAME, TY_CHAR) + call salloc (units, SZ_FNAME, TY_CHAR) + call salloc (errstr, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Get task parameters. + list = imtopenp ("input") + call clgstr ("records", Memc[image], SZ_FNAME) + call odr_openp (list, Memc[image]) + call clgstr ("output", Memc[output], SZ_FNAME) + call clgstr ("apertures", Memc[image], SZ_FNAME) + bandwidth = clgetr ("bandwidth") + bandsep = clgetr ("bandsep") + bswitch = clgetb ("beam_switch") + if (bswitch) + samestar = true + else + samestar = clgetb ("samestar") + if (clgetb ("interact")) + interactive = YES1 + else + interactive = NO3 + + # Expand the aperture list. + iferr (aps = rng_open (Memc[image], INDEF, INDEF, INDEF)) + call error (0, "Bad aperture list") + + call ext_load (ewaves, emags, enwaves) + + un = NULL + sh = NULL + obj = NULL + sky = NULL + obs = NULL + gp = NULL + nstds = 0 + while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) { + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + mw = smw_openim (im) + call shdr_open (im, mw, 1, 1, INDEFI, SHHDR, sh) + + if (DC(sh) == DCNO) { + call eprintf ("%s: No dispersion function\n") + call pargstr (Memc[image]) + call smw_close (MW(sh)) + call imunmap (IM(sh)) + next + } + + # Work in units of first spectrum. + if (un == NULL) { + call strcpy (UNITS(sh), Memc[units], SZ_FNAME) + un = un_open (Memc[units]) + unang = un_open ("Angstroms") + call un_ctranr (unang, un, Memr[ewaves], Memr[ewaves], enwaves) + } + if (IS_INDEF (IT(sh))) { + call printf ("%s: ") + call pargstr (Memc[image]) + call flush (STDOUT) + IT(sh) = clgetr ("exptime") + call imunmap (IM(sh)) + ifnoerr (im = immap (Memc[image], READ_WRITE, 0)) { + IM(sh) = im + call imseti (IM(sh), IM_WHEADER, YES) + call imaddr (IM(sh), "exptime", IT(sh)) + } else { + im = immap (Memc[image], READ_ONLY, 0) + IM(sh) = im + } + } + + do line = 1, SMW_NSPEC(mw) { + call shdr_open (im, mw, line, 1, INDEFI, SHDATA, sh) + if (!rng_elementi (aps, AP(sh))) + next + call shdr_units (sh, Memc[units]) + + if (!bswitch || OFLAG(sh) == OBJ) { + call printf ("%s%s(%d): %s\n") + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargi (AP(sh)) + call pargstr (TITLE(sh)) + call flush (STDOUT) + } + + if (IS_INDEF (AM(sh))) { + call clgstr ("observatory", Memc[observatory], SZ_FNAME) + call obsimopen (obs, im, Memc[observatory], NO, newobs, + obshead) + if (newobs) + call obslog (obs, "STANDARD", "latitude", STDOUT) + latitude = obsgetr (obs, "latitude") + iferr (call get_airm (RA(sh), DEC(sh), HA(sh), ST(sh), + latitude, AM(sh))) { + call printf ("%s: ") + call pargstr (Memc[image]) + call flush (STDOUT) + AM(sh) = clgetr ("airmass") + call imunmap (IM(sh)) + ifnoerr (im = immap (Memc[image], READ_WRITE, 0)) { + IM(sh) = im + call imseti (IM(sh), IM_WHEADER, YES) + call imaddr (IM(sh), "airmass", AM(sh)) + } else { + im = immap (Memc[image], READ_ONLY, 0) + IM(sh) = im + } + } + } + + for (i=0; i<nstds; i=i+1) { + std = Memi[stds+i] + if (STD_AP(std) == AP(sh)) + break + } + + # Allocate space for this aperture if not already done. + if (i >= nstds) { + if (nstds == 0) + call malloc (stds, 10, TY_INT) + else if (mod (nstds, 10) == 0) + call realloc (stds, nstds+10, TY_INT) + call salloc (std, STD_LEN, TY_STRUCT) + Memi[stds+i] = std + nstds = nstds + 1 + + STD_AP(std) = AP(sh) + STD_TYPE(std) = NONE + STD_SH(std) = NULL + STD_IFLAG(std) = interactive + STD_NWAVES(std) = 0 + + if (!samestar || i == 0) { + # Read calibration data + Memc[str] = EOS + repeat { + iferr (call std_gcalib (std, un)) { + j = errget (Memc[errstr], SZ_FNAME) + if (streq (Memc[errstr], Memc[str])) + call erract (EA_ERROR) + call strcpy (Memc[errstr], Memc[str], SZ_LINE) + call erract (EA_WARN) + next + } + break + } + } else { + CAL_NWAVES(std) = CAL_NWAVES(Memi[stds]) + CAL_WAVES(std) = CAL_WAVES(Memi[stds]) + CAL_DWAVES(std) = CAL_DWAVES(Memi[stds]) + CAL_MAGS(std) = CAL_MAGS(Memi[stds]) + } + + if (IS_INDEF (bandwidth)) { + do j = 1, CAL_NWAVES(std) { + wave = Memr[CAL_WAVES(std)+j-1] + dwave = Memr[CAL_DWAVES(std)+j-1] + call std_addband (std, wave, dwave, 0.) + } + } else { + wave = W0(sh) + bandwidth / 2 + dwave = W0(sh) + (SN(sh)-1) * WP(sh) - bandwidth / 2 + while (wave <= dwave) { + call std_addband (std, wave, bandwidth, 0.) + wave = wave + bandsep + } + } + } + + # The copying of SHDR structures and associated MWCS only + # occurs with beam switched data. + + if (bswitch) { + switch (STD_TYPE(std)) { + case NONE: + STD_TYPE(std) = OFLAG(sh) + call shdr_copy (sh, STD_SH(std), YES) + next + case SKY: + obj = sh + sky = STD_SH(std) + if (OFLAG(sh) == SKY) { + call eprintf ("%s[%d]: Object spectrum not found\n") + call pargstr (IMNAME(sky)) + call pargi (AP(sky)) + + call smw_close (MW(sky)) + call shdr_copy (sh, STD_SH(std), YES) + next + } + case OBJ: + obj = STD_SH(std) + sky = sh + if (OFLAG(sh) == OBJ) { + obj = STD_SH(std) + call eprintf ("%s[%d]: Sky spectrum not found\n") + call pargstr (IMNAME(obj)) + call pargi (AP(obj)) + + call smw_close (MW(obj)) + call shdr_copy (sh, STD_SH(std), YES) + next + } + } + } else { + obj = sh + sky = NULL + } + + # Generate a calibration table + call std_calib (obj, sky, std, gp, Memr[ewaves], Memr[emags], + enwaves) + call std_output (obj, sky, std, Memc[output]) + + if (interactive == YES1) { + if (STD_IFLAG(std) == NO3 || STD_IFLAG(std) == YES3) { + interactive = STD_IFLAG(std) + do i = 0, nstds-1 + STD_IFLAG(Memi[stds+i]) = interactive + } + if (interactive == NO3 && gp != NULL) + call gclose (gp) + } + + if (bswitch) { + call smw_close (MW(STD_SH(std))) + STD_TYPE(std) = NONE + } + } + + call smw_close (MW(sh)) + call imunmap (IM(sh)) + } + + if (un != NULL) { + call un_close (un) + call un_close (unang) + } + if (obs != NULL) + call obsclose (obs) + if (gp != NULL) + call gclose (gp) + do i = 0, nstds-1 { + std = Memi[stds+i] + obj = STD_SH(std) + switch (STD_TYPE(std)) { + case SKY: + call eprintf ("%s[%d]: Object spectrum not found\n") + call pargstr (IMNAME(obj)) + call pargi (AP(obj)) + case OBJ: + call eprintf ("%s[%d]: Sky spectrum not found\n") + call pargstr (IMNAME(obj)) + call pargi (AP(obj)) + } + if (obj != NULL) + call shdr_close (obj) + if (!samestar || i == 0) { + call mfree (CAL_WAVES(std), TY_REAL) + call mfree (CAL_DWAVES(std), TY_REAL) + call mfree (CAL_MAGS(std), TY_REAL) + } + call mfree (STD_WAVES(std), TY_REAL) + call mfree (STD_DWAVES(std), TY_REAL) + call mfree (STD_MAGS(std), TY_REAL) + call mfree (STD_FLUXES(std), TY_REAL) + } + call mfree (stds, TY_INT) + call mfree (ewaves, TY_REAL) + call mfree (emags, TY_REAL) + call shdr_close (sh) + call rng_close (aps) + call imtclose (list) + call sfree (sp) +end + + +# STD_CALIB -- Compute standard star calibrations + +procedure std_calib (obj, sky, std, gp, ewaves, emags, enwaves) + +pointer obj # Object pointer +pointer sky # Sky pointer +pointer std # Standard pointer +pointer gp # Graphics pointer +real ewaves[enwaves] # Extinction wavelengths +real emags[enwaves] # Extinction magnitudes +int enwaves # Extinction points + +int i, j, n, nwaves, wcs, key, newgraph +real wave, dwave, flux, wx1, wx2, wy +pointer sp, cmd, gt, waves, dwaves, fluxes, x, y + +real std_flux() +double shdr_wl() +int clgcur(), strdic(), clgwrd() +pointer gopen(), gt_init() +errchk gopen, std_output + +define beep_ 99 + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Sky subtract + if (sky != NULL) { + call shdr_rebin (sky, obj) + call asubr (Memr[SY(obj)], Memr[SY(sky)], Memr[SY(obj)], SN(obj)) + } + + # Remove extinction correction + if (EC(obj) == ECYES) { + x = SX(obj) + y = SY(obj) + n = SN(obj) + do i = 1, n { + call intrp (EXT_LOOKUP, ewaves, emags, enwaves, + Memr[x], flux, j) + Memr[y] = Memr[y] * 10.0 ** (-0.4 * flux * AM(obj)) + x = x + 1 + y = y + 1 + } + } + + nwaves = STD_NWAVES(std) + waves = STD_WAVES(std) + dwaves = STD_DWAVES(std) + fluxes = STD_FLUXES(std) + do i = 0, nwaves-1 { + wave = Memr[waves+i] + dwave = Memr[dwaves+i] + Memr[fluxes+i] = std_flux (obj, wave, dwave, ewaves, emags, enwaves) + } + + # Plot spectrum if user wants to see whats happening + if (STD_IFLAG(std) == NO1 || STD_IFLAG(std) == YES1) { + call printf ("%s[%d]: Edit bandpasses? ") + call pargstr (IMNAME(obj)) + call pargi (AP(obj)) + STD_IFLAG(std) = clgwrd ("answer", Memc[cmd], SZ_FNAME, ANSWERS) + } + + i = STD_IFLAG(std) + if (i==YES1||i==Y2||i==YES2||i==YES3) { + if (gp == NULL) { + call clgstr ("graphics", Memc[cmd], SZ_FNAME) + gp = gopen (Memc[cmd], NEW_FILE, STDGRAPH) + } + gt = gt_init() + call gt_sets (gt, GTTITLE, TITLE(obj)) + call gt_sets (gt, GTPARAMS, IMNAME(obj)) + call gt_sets (gt, GTXLABEL, LABEL(obj)) + call gt_sets (gt, GTXUNITS, UNITS(obj)) + call gt_sets (gt, GTYLABEL, "instrumental flux") + call gt_sets (gt, GTTYPE, "line") + + key = 'r' + repeat { + switch (key) { + case '?': + call gpagefile (gp, KEY, PROMPT) + case ':': + if (Memc[cmd] == '/') + call gt_colon (Memc[cmd], gp, gt, newgraph) + else { + switch (strdic (Memc[cmd],Memc[cmd],SZ_LINE,"|show|")) { + case 1: + call mktemp ("std", Memc[cmd], SZ_LINE) + call std_output (obj, sky, std, Memc[cmd]) + call gpagefile (gp, Memc[cmd], "standard star data") + call delete (Memc[cmd]) + default: + goto beep_ + } + } + case 'a': + call printf ("a again:\n") + i = clgcur ("cursor", wx2, wy, wcs, key, Memc[cmd], SZ_LINE) + call printf ("\n") + if (wx1 == wx2) { + call printf ("\07Two cursor positions required") + goto beep_ + } + + # Create artificial standard wavelength and bandpass + wave = (wx1 + wx2) / 2.0 + dwave = wx2 - wx1 + flux = std_flux (obj, wave, dwave, ewaves, emags, enwaves) + call std_addband (std, wave, dwave, flux) + flux = flux / abs (shdr_wl (obj, double(wx1)) - + shdr_wl (obj, double (wx2))) + call gmark (gp, wave, flux, GM_BOX, -dwave, 3.) + + nwaves = STD_NWAVES(std) + waves = STD_WAVES(std) + dwaves = STD_DWAVES(std) + fluxes = STD_FLUXES(std) + case 'd': + dwave = MAX_REAL + do i = 0, nwaves-1 { + wave = Memr[waves+i] + if (abs (wx1 - wave) < dwave) { + dwave = abs (wx1 - wave) + j = i + } + } + wave = Memr[waves+j] + dwave = Memr[dwaves+j] + flux = Memr[fluxes+j] + flux = flux / abs (shdr_wl (obj, double(wave-dwave/2)) - + shdr_wl (obj, double (wave+dwave/2))) + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, wave, flux, GM_BOX, -dwave, 3.) + call gseti (gp, G_PMLTYPE, 1) + call gscur (gp, wave, flux) + call std_delband (std, j) + case 'q': + break + case 'I': + call fatal (0, "Interrupt") + case 'r': + newgraph = YES + case 'w': + call gt_window (gt, gp, "cursor", newgraph) + default: # Invalid keystroke +beep_ call printf ("\007") + } + + if (newgraph == YES) { + call std_graph (obj, std, gp, gt, YES) + newgraph = NO + } + } until (clgcur ("cursor",wx1,wy,wcs,key,Memc[cmd],SZ_LINE) == EOF) + call gt_free (gt) + } + + call sfree (sp) +end + + +# STD_OUTPUT -- Output standard star data. +# For now we do this in Angstroms. + +procedure std_output (obj, sky, std, output) + +pointer obj # Object pointer +pointer sky # Sky pointer +pointer std # Standard pointer +char output[ARB] # Output file name + +int i, fd, open() +real wave, dwave, mag, flux, fnuzero, flambda, clgetr() +pointer unang, un_open() +errchk open, un_open, un_ctranr, std_units + +begin + fd = open (output, APPEND, TEXT_FILE) + call fprintf (fd, "[%s]") + call pargstr (IMNAME(obj)) + if (sky != NULL) { + call fprintf (fd, "-[%s]") + call pargstr (IMNAME(sky)) + } + + unang = un_open ("Angstroms") + call un_ctranr (UN(obj), unang, W0(obj), wave, 1) + call un_ctranr (UN(obj), unang, W0(obj)+(SN(obj)-1)*WP(obj), dwave, 1) + call fprintf (fd, " %d %d %.2f %5.3f %9.3f %9.3f %s\n") + call pargi (AP(obj)) + call pargi (SN(obj)) + call pargr (IT(obj)) + call pargr (AM(obj)) + #call pargr (W0(obj)) + #call pargr (W0(obj) + (SN(obj)-1) * WP(obj)) + call pargr (wave) + call pargr (dwave) + call pargstr (TITLE(obj)) + + fnuzero = clgetr ("fnuzero") + do i = 0, STD_NWAVES(std)-1 { + wave = Memr[STD_WAVES(std)+i] + dwave = Memr[STD_DWAVES(std)+i] + mag = Memr[STD_MAGS(std)+i] + flux = Memr[STD_FLUXES(std)+i] + if (flux == 0.) + next + + call std_units (UN(obj), unang, wave, dwave, 1) + flambda = fnuzero * 10. ** (-0.4 * mag) * VLIGHT / wave**2 + call fprintf (fd, "%8.2f %12.5g %8.3f %12.5g\n") + call pargr (wave) + call pargr (flambda) + call pargr (dwave) + call pargr (flux) + } + call close (fd) + + call un_close (unang) +end + + +# STD_FLUX -- Add up the flux in a given bandpass centered on a given +# wavelength. The bandpass must be entirely within the data. +# A correction for differential extinction across the bandpass is made +# by applying the extinction correction and then removing the correction +# at the bandpass center + +real procedure std_flux (sh, wave, dwave, ewaves, emags, enwaves) + +pointer sh # Spectrum +real wave # Bandpass wavelength +real dwave # Bandpass width +real ewaves[enwaves] # Extinction wavelengths +real emags[enwaves] # Extinction magnitudes +int enwaves # Extinction points + +real flux # Bandpass flux + +int i, i1, i2, ierr +real a, e, ec, x1, x2 +double w1, w2, w3, w4, shdr_lw(), shdr_wl() +pointer x, y + +begin + # Determine bandpass limits in pixel and return if out of bounds. + w1 = wave - dwave / 2. + w2 = wave + dwave / 2. + w3 = shdr_lw (sh, 0.5D0) + w4 = shdr_lw (sh, double (SN(sh)+0.5)) + if (w1 < min (w3, w4) || w2 > max (w3, w4)) + return (0.) + + a = shdr_wl (sh, w1) + x2 = shdr_wl (sh, w2) + x1 = min (a, x2) + x2 = max (a, x2) + i1 = nint (x1) + i2 = nint (x2 - 0.00001) + if (x1 == x2 || i1 < 1 || i2 > SN(sh)) + return (0.) + + a = AM(sh) + x = SX(sh) + i1 - 1 + y = SY(sh) + i1 - 1 + + call intrp (EXT_LOOKUP, ewaves, emags, enwaves, wave, ec, ierr) + call intrp (EXT_LOOKUP, ewaves, emags, enwaves, Memr[x], e, ierr) + + if (i1 == i2) { + flux = (x2-x1) * Memr[y] * 10.0 ** (0.4 * a * (e - ec)) + return (flux) + } + + flux = (i1+0.5-x1) * Memr[y] * 10.0 ** (0.4 * a * (e - ec)) + x = x + 1 + y = y + 1 + + for (i=i1+1; i<=i2-1; i=i+1) { + call intrp (EXT_LOOKUP, ewaves, emags, enwaves, Memr[x], e, ierr) + flux = flux + Memr[y] * 10.0 ** (0.4 * a * (e - ec)) + x = x + 1 + y = y + 1 + } + + call intrp (EXT_LOOKUP, ewaves, emags, enwaves, Memr[x], e, ierr) + flux = flux + (x2-i2+0.5) * Memr[y] * 10.0 ** (0.4 * a * (e - ec)) + + return (flux) +end + + +# STD_ADDBAND -- Add a standard bandpass + +procedure std_addband (std, wave, dwave, flux) + +pointer std # Pointer to standard star data +real wave # Wavelength to be added +real dwave # Bandpass to be added +real flux # Flux to be added + +int i, nwaves +real mag +pointer waves, dwaves, mags, fluxes + +begin + nwaves = STD_NWAVES(std) + if (nwaves == 0) { + call malloc (STD_WAVES(std), 10, TY_REAL) + call malloc (STD_DWAVES(std), 10, TY_REAL) + call malloc (STD_MAGS(std), 10, TY_REAL) + call malloc (STD_FLUXES(std), 10, TY_REAL) + } else if (mod (nwaves, 10) == 0) { + call realloc (STD_WAVES(std), nwaves+10, TY_REAL) + call realloc (STD_DWAVES(std), nwaves+10, TY_REAL) + call realloc (STD_MAGS(std), nwaves+10, TY_REAL) + call realloc (STD_FLUXES(std), nwaves+10, TY_REAL) + } + + call intrp (MAG_LOOKUP, Memr[CAL_WAVES(std)], Memr[CAL_MAGS(std)], + CAL_NWAVES(std), wave, mag, i) + + waves = STD_WAVES(std) + dwaves = STD_DWAVES(std) + mags = STD_MAGS(std) + fluxes = STD_FLUXES(std) + for (i=nwaves; (i>0)&&(Memr[waves+i-1]>wave); i=i-1) { + Memr[waves+i] = Memr[waves+i-1] + Memr[dwaves+i] = Memr[dwaves+i-1] + Memr[mags+i] = Memr[mags+i-1] + Memr[fluxes+i] = Memr[fluxes+i-1] + } + Memr[waves+i] = wave + Memr[dwaves+i] = dwave + Memr[mags+i] = mag + Memr[fluxes+i] = flux + STD_NWAVES(std) = nwaves + 1 +end + + +# STD_DELBAND -- Delete a bandpass + +procedure std_delband (std, band) + +pointer std # Pointer to standard star data +int band # Band to be deleted + +int i, nwaves +pointer waves, dwaves, mags, fluxes + +begin + nwaves = STD_NWAVES(std) + waves = STD_WAVES(std) + dwaves = STD_DWAVES(std) + mags = STD_MAGS(std) + fluxes = STD_FLUXES(std) + for (i=band+1; i<nwaves; i=i+1) { + Memr[waves+i-1] = Memr[waves+i] + Memr[dwaves+i-1] = Memr[dwaves+i] + Memr[mags+i-1] = Memr[mags+i] + Memr[fluxes+i-1] = Memr[fluxes+i] + } + nwaves = nwaves - 1 + + STD_NWAVES(std) = nwaves + if (nwaves == 0) { + call mfree (STD_WAVES(std), TY_REAL) + call mfree (STD_DWAVES(std), TY_REAL) + call mfree (STD_MAGS(std), TY_REAL) + call mfree (STD_FLUXES(std), TY_REAL) + } +end + + +# STD_GRAPH -- Graph the spectrum and standard star calibration points. + +procedure std_graph (sh, std, gp, gt, clear) + +pointer sh # Spectrum pointer +pointer std # Standard star data +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +int clear # Clear flag + +int i +real dw, wave, dwave, flux +double shdr_wl() + +begin + if (clear == YES) { + call gclear (gp) + call gascale (gp, Memr[SX(sh)], SN(sh), 1) + call gascale (gp, Memr[SY(sh)], SN(sh), 2) + call gt_swind (gp, gt) + call gt_labax (gp, gt) + } + + call gt_plot (gp, gt, Memr[SX(sh)], Memr[SY(sh)], SN(sh)) + + do i = 0, STD_NWAVES(std)-1 { + wave = Memr[STD_WAVES(std)+i] + dwave = Memr[STD_DWAVES(std)+i] + flux = Memr[STD_FLUXES(std)+i] + if (flux == 0.) + next + dw = abs (shdr_wl (sh, double(wave-dwave/2)) - + shdr_wl (sh, double (wave+dwave/2))) + flux = flux / dw + call gmark (gp, wave, flux, GM_BOX, -dwave, 3.) + } +end + + +# STD_GCALIB -- Get calibration data in desired units. + +procedure std_gcalib (std, un) + +pointer std #I Standard pointer +pointer un #I Desired units pointer + +pointer unang, un_open() +errchk getcalib, std_units + +begin + call getcalib (CAL_WAVES(std), CAL_DWAVES(std), CAL_MAGS(std), + CAL_NWAVES(std)) + + # Cnvert to desired units. + unang = un_open ("Angstroms") + call std_units (unang, un, + Memr[CAL_WAVES(std)], Memr[CAL_DWAVES(std)], CAL_NWAVES(std)) + call un_close (unang) +end + + +# STD_UNITS -- Convert bandpass information to different units. + +procedure std_units (unin, unout, center, width, n) + +pointer unin #I Input units +pointer unout #I Output units +real center[ARB] #U Bandpass centers +real width[ARB] #U Bandpass widths +int n #I Number of bandpasses + +int i +real x1, x2 +bool un_compare() +errchk un_ctranr + + +begin + if (un_compare (unin, unout)) + return + + do i = 1, n { + x1 = center[i] - width[i] / 2 + x2 = center[i] + width[i] / 2 + call un_ctranr (unin, unout, x1, x1, 1) + call un_ctranr (unin, unout, x2, x2, 1) + center[i] = (x1 + x2) / 2. + width[i] = abs (x1 - x2) + } +end |