diff options
Diffstat (limited to 'noao/onedspec/t_sbands.x')
| -rw-r--r-- | noao/onedspec/t_sbands.x | 585 |
1 files changed, 585 insertions, 0 deletions
diff --git a/noao/onedspec/t_sbands.x b/noao/onedspec/t_sbands.x new file mode 100644 index 00000000..d3c3d1e2 --- /dev/null +++ b/noao/onedspec/t_sbands.x @@ -0,0 +1,585 @@ +include <error.h> +include <smw.h> + +# Band structure +define LEN_BAND 9 # length of structure +define BAND_ID Memi[$1] # ptr to band id string +define BAND_FILTER Memi[$1+1] # ptr to filter string +define BAND_WC Memd[P2D($1+2)] # center wavelength +define BAND_DW Memd[P2D($1+4)] # wavelength width +define BAND_FN Memi[$1+6] # no. of filter points +define BAND_FW Memi[$1+7] # ptr to filter wavelengths +define BAND_FR Memi[$1+8] # ptr to filter responses + +# Multiple bands for indices and equivalent widths. +define NBANDS 3 # maximum number of bands +define BAND1 1 +define BAND2 2 +define BAND3 3 +define BAND Memi[$1+($2-1)*NBANDS+($3-1)] + + +# T_SBANDS -- Compute band fluxes, indices, and equivalent widths. +# A list of bandpasses is supplied in a text file, and all of them are applied +# to each spectrum in the list. The output is written to an output file +# in multicolumn format. + +procedure t_sbands () + +pointer inlist # Input list of spectra +pointer output # Output file name +pointer fbands # Band file name +pointer apertures # Aperture list string +bool norm # Normalize bands by response? +bool mag # Output magnitudes instead of fluxes? +double magzero # Magnitude zeropoint for magnitude output +bool verbose # Verbose header? + +int i, nbands, nsubbands, nimages, fd +pointer bands, aps, im, smw, sh +pointer sp, input + +int open(), imtgetim() +bool clgetb(), rng_elementi() +double clgetd() +pointer imtopenp(), immap(), smw_openim(), rng_open() + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (fbands, SZ_FNAME, TY_CHAR) + call salloc (apertures, SZ_LINE, TY_CHAR) + + # Get task parameters. + inlist = imtopenp ("input") + call clgstr ("output", Memc[output], SZ_FNAME) + call clgstr ("bands", Memc[fbands], SZ_FNAME) + call clgstr ("apertures", Memc[apertures], SZ_LINE) + norm = clgetb ("normalize") + mag = clgetb ("mag") + magzero = clgetd ("magzero") + verbose = clgetb ("verbose") + + # Read bands from the band file. + fd = open (Memc[fbands], READ_ONLY, TEXT_FILE) + call sb_bands (fd, bands, nbands, nsubbands) + call close (fd) + + # Open the aperture list. + iferr (aps = rng_open (Memc[apertures], INDEF, INDEF, INDEF)) + call error (1, "Bad aperture list") + + # Loop over the input spectra. + fd = 0 + nimages = 0 + while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) { + nimages = nimages + 1 + + # Open the input image and get the WCS. + iferr (im = immap (Memc[input], READ_ONLY, 0)) { + call erract (EA_WARN) + next + } + iferr (smw = smw_openim (im)) { + call imunmap (im) + call erract (EA_WARN) + next + } + + # Open output file and write a verbose header if desired. + # It is delayed until now to avoid output if an error occurs + # such as image not found. + + if (nimages == 1) { + fd = open (Memc[output], APPEND, TEXT_FILE) + if (verbose) + call sb_header (fd, norm, mag, magzero, + Memc[fbands], bands, nbands, nsubbands) + } + + # Measure selected apertures. + do i = 1, SMW_NSPEC(smw) { + call shdr_open (im, smw, i, 1, INDEFI, SHHDR, sh) + if (!rng_elementi (aps, AP(sh))) + next + call shdr_open (im, smw, i, 1, INDEFI, SHDATA, sh) + + call sb_proc (fd, sh, bands, nbands, norm, mag, magzero) + } + + # Finish with image. + call shdr_close (sh) + call smw_close (smw) + call imunmap (im) + } + + # Finish up. + call sb_free (bands, nbands) + if (fd != 0) + call close (fd) + call imtclose (inlist) + call sfree (sp) +end + + +# SB_BANDS - Read bands from the band file and put them into an array +# of band pointers. + +procedure sb_bands (fd, bands, nbands, nsubbands) + +int fd #I Bandpass file descriptor +pointer bands #O Bandpass table descriptor +int nbands #O Number of bandpasses +int nsubbands #O Number of individual bands + +bool bandok +int ip +double center, width +pointer sp, line, id, filter + +int getline(), ctowrd(), ctod() + +begin + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + call salloc (id, SZ_FNAME, TY_CHAR) + call salloc (filter, SZ_FNAME, TY_CHAR) + + # Read the bands. If the first band is not seen + # skip the line. Check for 1, 2, or 3 bandpasses. + # Can't use fscan() because fscan() will be called later to + # read any filter file. + + bands = NULL + nbands = 0 + nsubbands = 0 + while (getline (fd, Memc[line]) != EOF) { + ip = 1 + bandok = (ctowrd (Memc[line], ip, Memc[id], SZ_FNAME) > 0) + bandok = (bandok && ctod (Memc[line], ip, center) > 0) + bandok = (bandok && ctod (Memc[line], ip, width) > 0) + bandok = (bandok && ctowrd (Memc[line],ip,Memc[filter],SZ_FNAME)>0) + if (!bandok || Memc[id] == '#') + next + + # Allocate and reallocate the array of band pointers. + if (nbands == 0) + call malloc (bands, 10 * NBANDS, TY_POINTER) + else if (mod (nbands, 10) == 0) + call realloc (bands, (nbands + 10) * NBANDS, TY_POINTER) + nbands = nbands + 1 + + call sb_alloc (BAND(bands,nbands,BAND1), + Memc[id], Memc[filter], center, width) + nsubbands = nsubbands + 1 + + bandok = (ctowrd (Memc[line], ip, Memc[id], SZ_FNAME) > 0) + bandok = (bandok && ctod (Memc[line], ip, center) > 0) + bandok = (bandok && ctod (Memc[line], ip, width) > 0) + bandok = (bandok && ctowrd (Memc[line],ip,Memc[filter],SZ_FNAME)>0) + if (bandok) { + call sb_alloc (BAND(bands,nbands,BAND2), + Memc[id], Memc[filter], center, width) + nsubbands = nsubbands + 1 + } else + BAND(bands,nbands,BAND2) = NULL + + bandok = (ctowrd (Memc[line], ip, Memc[id], SZ_FNAME) > 0) + bandok = (bandok && ctod (Memc[line], ip, center) > 0) + bandok = (bandok && ctod (Memc[line], ip, width) > 0) + bandok = (bandok && ctowrd (Memc[line],ip,Memc[filter],SZ_FNAME)>0) + if (bandok) { + call sb_alloc (BAND(bands,nbands,BAND3), + Memc[id], Memc[filter], center, width) + nsubbands = nsubbands + 1 + } else + BAND(bands,nbands,BAND3) = NULL + } + + call sfree (sp) +end + + +# SB_ALLOC -- Allocate a band structure. + +procedure sb_alloc (band, id, filter, center, width) + +pointer band #O Band pointer +char id[ARB] #I Band id +char filter[ARB] #I Band filter +double center #I Band wavelength +double width #I Band width + +int fn, fd, strlen(), open(), fscan(), nscan() +double w, r +pointer fw, fr +bool streq() +errchk open() + +begin + call calloc (band, LEN_BAND, TY_STRUCT) + call malloc (BAND_ID(band), strlen(id), TY_CHAR) + call malloc (BAND_FILTER(band), strlen(filter), TY_CHAR) + call strcpy (id, Memc[BAND_ID(band)], ARB) + call strcpy (filter, Memc[BAND_FILTER(band)], ARB) + BAND_WC(band) = center + BAND_DW(band) = width + BAND_FN(band) = 0 + BAND_FW(band) = NULL + BAND_FR(band) = NULL + + if (streq (filter, "none")) + return + + # Read the filter file. + fd = open (filter, READ_ONLY, TEXT_FILE) + fn = 0 + while (fscan (fd) != EOF) { + call gargd (w) + call gargd (r) + if (nscan() != 2) + next + if (fn == 0) { + call malloc (fw, 100, TY_DOUBLE) + call malloc (fr, 100, TY_DOUBLE) + } else if (mod (fn, 100) == 0) { + call realloc (fw, fn+100, TY_DOUBLE) + call realloc (fr, fn+100, TY_DOUBLE) + } + Memd[fw+fn] = w + Memd[fr+fn] = r + fn = fn + 1 + } + call close (fd) + + BAND_FN(band) = fn + BAND_FW(band) = fw + BAND_FR(band) = fr +end + + +# SB_FREE -- Free band structures. + +procedure sb_free (bands, nbands) + +pointer bands #I bands descriptor +int nbands #I number of bands + +int i, j +pointer band + +begin + do i = 1, nbands { + do j = 1, NBANDS { + band = BAND(bands,i,j) + if (band != NULL) { + call mfree (BAND_ID(band), TY_CHAR) + call mfree (BAND_FILTER(band), TY_CHAR) + call mfree (BAND_FW(band), TY_DOUBLE) + call mfree (BAND_FR(band), TY_DOUBLE) + call mfree (band, TY_STRUCT) + } + } + } + call mfree (bands, TY_POINTER) +end + + +# SB_HEADER -- Print output header. + +procedure sb_header (fd, norm, mag, magzero, fbands, bands, nbands, nsubbands) + +pointer fd #I Output file descriptor +bool norm #I Normalization flag +bool mag #I Magnitude flag +double magzero #I Magnitude zeropoint +char fbands[ARB] #I Band file +pointer bands #I Pointer to array of bands +int nbands #I Number of bands +int nsubbands #I Number of subbands + +int i, j +pointer sp, str, band + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Output a banner and task parameters. + call sysid (Memc[str], SZ_LINE) + call fprintf (fd, "\n# SBANDS: %s\n# ") + call pargstr (Memc[str]) + if (fbands[1] != EOS) { + call fprintf (fd, " bands = %s,") + call pargstr (fbands) + } + call fprintf (fd, " norm = %b, mag = %b") + call pargb (norm) + call pargb (mag) + if (mag) { + call fprintf (fd, ", magzero = %.2f") + call pargd (magzero) + call strcpy ("mag", Memc[str], SZ_LINE) + } else + call strcpy ("flux", Memc[str], SZ_LINE) + + # Output the bands. + call fprintf (fd, "\n# %10s %10s %10s %10s\n") + call pargstr ("band") + call pargstr ("filter") + call pargstr ("wavelength") + call pargstr ("width") + do i = 1, nbands { + do j = 1, NBANDS { + band = BAND(bands,i,j) + if (band == NULL) + next + call fprintf (fd, "# %10s %10s %10g %10g\n") + call pargstr (Memc[BAND_ID(band)]) + call pargstr (Memc[BAND_FILTER(band)]) + call pargd (BAND_WC(band)) + call pargd (BAND_DW(band)) + } + } + + # Output column headings. + call fprintf (fd, + "#\n# %24s %7.7s %11.11s") + call pargstr ("spectrum") + call pargstr ("band") + call pargstr (Memc[str]) + if (nsubbands > nbands) { + call fprintf (fd, " %7.7s %11.11s %9.9s %9.9s") + call pargstr ("band") + call pargstr (Memc[str]) + call pargstr ("index") + call pargstr ("eqwidth") + } + call fprintf (fd, "\n") + + call sfree (sp) +end + + +# SB_PROC -- Measure the band fluxes and possibly a band index and eq. width. + +procedure sb_proc (fd, sh, bands, nbands, norm, mag, magzero) + +int fd #I Output file descriptor +pointer sh #I Spectrum descriptor +pointer bands #I Bandpass table pointer +int nbands #I Number of bandpasses +bool norm #I Normalize? +bool mag #I Magnitude output? +double magzero #I Magnitude zero point + +int i +double flux, contval, index, eqwidth +double flux1, norm1, flux2, norm2, flux3, norm3, a, b +pointer sp, imname, band1, band2, band3 + +begin + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + + call sprintf (Memc[imname], SZ_FNAME, "%s%s(%d)") + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargi (AP(sh)) + + # Loop over all bandpasses + do i = 1, nbands { + # Measure primary band flux, normalize, and print result. + band1 = BAND(bands,i,BAND1) + call sb_flux (sh, band1, flux1, norm1) + if (IS_INDEFD(flux1)) + next + + if (norm) { + flux1 = flux1 / norm1 + norm1 = 1 + } + if (mag && flux1 > 0.) + flux = magzero - 2.5 * log10 (flux1) + else + flux = flux1 + + call fprintf (fd, "%26s %7.7s %11.6g") + call pargstr (Memc[imname]) + call pargstr (Memc[BAND_ID(band1)]) + call pargd (flux) + + # Measure the alternate band fluxes and compute and output + # the band index and equivalent width. + + band2 = BAND(bands,i,BAND2) + band3 = BAND(bands,i,BAND3) + call sb_flux (sh, band2, flux2, norm2) + call sb_flux (sh, band3, flux3, norm3) + if (IS_INDEFD(flux2) && IS_INDEFD(flux3)) { + call fprintf (fd, "\n") + next + } + + if (norm) { + if (!IS_INDEFD(flux2)) { + flux2 = flux2 / norm2 + norm2 = 1 + } + if (!IS_INDEFD(flux3)) { + flux3 = flux3 / norm3 + norm3 = 1 + } + } + + contval = INDEFD + index = INDEFD + eqwidth = INDEFD + if (!IS_INDEFD(flux2) && !IS_INDEFD(flux3)) { + # Interpolate to the center of the primary band. + a = (flux2 / norm2 - flux3 / norm3) / + (BAND_WC(band2) - BAND_WC(band3)) + b = flux2 / norm2 - a * BAND_WC(band2) + contval = (a * BAND_WC(band1) + b) * norm1 + call fprintf (fd, " %7.7s") + call pargstr ("cont") + } else if (!IS_INDEFD(flux2)) { + contval = flux2 + call fprintf (fd, " %7.7s") + call pargstr (Memc[BAND_ID(band2)]) + } else if (!IS_INDEFD(flux3)) { + contval = flux3 + call fprintf (fd, " %7.7s") + call pargstr (Memc[BAND_ID(band3)]) + } + + if (mag && contval > 0.) + flux = magzero - 2.5 * log10 (contval) + else + flux = contval + call fprintf (fd, " %11.6g") + call pargd (flux) + + if (flux1 > 0. && contval > 0.) { + index = flux1 / contval + eqwidth = (1 - index) * BAND_DW(band1) + } + if (mag) { + if (!IS_INDEFD(contval) && contval > 0.) + contval = magzero - 2.5 * log10 (contval) + if (!IS_INDEFD(index)) + index = -2.5 * log10 (index) + } + + call fprintf (fd, " %9.6g %9.6g\n") + call pargd (index) + call pargd (eqwidth) + } + + # Flush output and finish up. + call flush (fd) + call sfree (sp) +end + + +# SB_FLUX - Compute the flux and total response in a given band. +# Return INDEF if the band is outside of the spectrum. + +procedure sb_flux (sh, band, flux, norm) + +pointer sh #I spectrum descriptor +pointer band #I band descriptor +double flux #O flux +double norm #O normalization + +int i, i1, i2 +double a, b, w1, w2, x1, x2, wt +pointer x, y +double sb_filter(), shdr_wl() + +begin + # Return if no band is defined. + flux = INDEFD + norm = 1 + if (band == NULL) + return + + # Check end points. + a = BAND_WC(band) - BAND_DW(band) / 2. + b = BAND_WC(band) + BAND_DW(band) / 2. + w1 = min (a, b) + w2 = max (a, b) + a = shdr_wl (sh, w1) + b = shdr_wl (sh, w2) + x1 = min (a, b) + x2 = max (a, b) + i1 = nint (x1) + i2 = nint (x2) + if (x1 == x2 || i1 < 1 || i2 > SN(sh)) + return + + x = SX(sh) + i1 - 1 + y = SY(sh) + i1 - 1 + + if (i1 == i2) { + wt = sb_filter (double(Memr[x]), band) * (x2 - x1) + flux = wt * Memr[y] + norm = wt + } else { + wt = sb_filter (double(Memr[x]), band) * (i1 + 0.5 - x1) + flux = wt * Memr[y] + norm = wt + x = x + 1 + y = y + 1 + for (i = i1 + 1; i <= i2 - 1; i = i + 1) { + wt = sb_filter (double(Memr[x]), band) + flux = flux + wt * Memr[y] + norm = norm + wt + x = x + 1 + y = y + 1 + } + wt = sb_filter (double(Memr[x]), band) * (x2 - i2 + 0.5) + flux = flux + wt * Memr[y] + norm = norm + wt + } +end + + +# SB_FILTER -- Given a filter array interpolate to the specified wavelength. + +double procedure sb_filter (w, band) + +double w # Wavelength desired +pointer band # Band pointer + +int i, n +double x1, x2 +pointer x, y + +begin + n = BAND_FN(band) + if (n == 0) + return (1.) + + x = BAND_FW(band) + y = BAND_FR(band) + x1 = Memd[x] + x2 = Memd[x+n-1] + + if (w <= x1) + return (Memd[y]) + else if (w >= x2) + return (Memd[y+n-1]) + + if ((w - x1) < (x2 - w)) + for (i = 1; w > Memd[x+i]; i=i+1) + ; + else + for (i = n - 1; w < Memd[x+i-1]; i=i-1) + ; + + x1 = Memd[x+i-1] + x2 = Memd[x+i] + return ((w - x1) / (x2 - x1) * (Memd[y+i] - Memd[y+i-1]) + Memd[y+i-1]) +end |