diff options
Diffstat (limited to 'noao/twodspec/apextract/apextract.x')
| -rw-r--r-- | noao/twodspec/apextract/apextract.x | 1834 |
1 files changed, 1834 insertions, 0 deletions
diff --git a/noao/twodspec/apextract/apextract.x b/noao/twodspec/apextract/apextract.x new file mode 100644 index 00000000..176ab251 --- /dev/null +++ b/noao/twodspec/apextract/apextract.x @@ -0,0 +1,1834 @@ +include <error.h> +include <imhdr.h> +include <mach.h> +include <math/iminterp.h> +include <pkg/gtools.h> +include "apertures.h" + +# Background fitting types +define BACKGROUND "|none|average|median|minimum|fit|" +define B_NONE 1 +define B_AVERAGE 2 +define B_MEDIAN 3 +define B_MINIMUM 4 +define B_FIT 5 + +# Weight types +define WEIGHTS "|none|variance|" +define W_NONE 1 +define W_VARIANCE 2 + +# Profile fitting algorithms +define P_FIT "|fit1d|fit2d|" +define P_FIT1D 1 +define P_FIT2D 2 + +# Output formats +define FORMATS "|onedspec|multispec|echelle|strip|normalize|flatten\ + |ratio|difference|fit|noise|" +define ONEDSPEC 1 # Individual 1D spectra +define MULTISPEC 2 # Multiple spectra +define ECHELLE 3 # Echelle spectra +define STRIP 4 # Strip spectra +define NORM 5 # Normalized spectra +define FLAT 6 # Flat spectra +define RATIO 7 # Ratio of data to model +define DIFF 8 # Difference of data and model +define FIT 9 # Model +define NOISE 10 # Noise calculation + + +# AP_EXTRACT -- Extract spectra by a weighted sum across the apertures. +# +# This routine does clobber checks on the output images, manages the I/O +# from the input image in as big of pieces as possible, and loops through +# each aperture calling routines to determine the sky, do any fitting and +# extraction, and output the spectra. +# The extraction may be either a simple, unweighted extraction +# which is very fast or a weighted extraction using CCD noise +# parameters. The weights require dividing out the basic spectrum and +# smoothing the 2D spectral profile. The general approach of variance +# weighting is described by K. Horne (PASP V98, P609, 1986). The +# smoothing has two algorithms, fitting columns or lines parallel to the +# dispersion axis for nearly aligned spectra or fitting a 2D function +# using a method given by T. Marsh (PASP V101, P1032, 1989). The profile +# may also be used to reject cosmic rays by iteration. +# +# The extractions require enough memory to get at least one aperture plus +# background (if needed) into memory. If possible the region containing +# all the apertures is read into memory. The target maximum amount of +# memory is set by the maxmimum size returned by BEGMEM and the +# appropriate working set size is requested. The optimal size can be +# tuned through BEGMEM, which references a machine dependent include +# file, if needed. The algorithm should work well (minimize I/O as well +# as paging) in all cases but very large image formats with highly tilted +# spectra (where aperture extraction along the image axes is not really +# appropriate). These memory requirements were chosen to minimize image +# I/O and because the variance weighted algorithms need to make multiple +# passes through the image. In principle simple, unweighted extractions +# with no sky smoothing can be done sequentially but this was not done in +# order to use nearly the same code for both weighted and unweighted +# cases. +# +# If using variance weighting and a profile image is given then it is used +# to determine the profile which is then applied to the target image +# during the final extraction. If the same profile image is used multiple +# times it would be more efficient to store the profile but then issues +# of consistency arise. For now this possible feature is not implemented. + +procedure ap_extract (input, output, format, profiles, aps, naps) + +char input[SZ_FNAME] # Input image +char output[SZ_FNAME] # Output image (optional root name) +char format[SZ_LINE] # Output format +char profiles[SZ_FNAME] # Profile filename (optional) +pointer aps[ARB] # Apertures +int naps # Number of apertures + +# CL parameters +int fmt # Output format +int bkg # Background type +int weights # Extraction weights +int pfit # Profile fitting algorithm +bool clean # Reject cosmic rays? +real gain # Photon/DN gain +real rdnoise # Read out noise +int nsubaps # Number of subapertures +int interptype # Edge interpolation type + +int i, j, k, napsex, aaxis, baxis, namax, na, nb, na1, interpbuf +int amin, amax, bmin, bmax +int new_size, old_size, max_size, best_size +real cmin, cmax, xmin, xmax, shift +pointer sp, str, bkgstr, wtstr, cleanstr, apsex +pointer a, b, c, astart, spec, specsky, specsig, raw, profile +pointer a1, a2, b1, b2, c1, c2, im, pim, ap, cv, ic, dbuf, pbuf, sbuf, svar, ptr +pointer asi + +bool clgetb(), apgetb(), strne() +int apgeti(), apgwrd(), begmem(), ap_check() +real apgimr(), ap_cveval(), ic_getr() +pointer ap_immap(), imgs2r(), imgl2r() +errchk salloc, malloc, ap_immap, imgs2r, imgl2r, asiinit +errchk ap_check, ap_skyeval, ap_profile, ap_variance, ap_output, apgimr + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + napsex = 0 + do i = 1, naps + if (AP_SELECT(aps[i]) == YES) + napsex = napsex + 1 + if (napsex == 0) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - No apertures defined for %s") + call pargstr (input) + call ap_log (Memc[str], YES, NO, YES) + call sfree (sp) + return + } + + call salloc (bkgstr, SZ_FNAME, TY_CHAR) + call salloc (wtstr, SZ_FNAME, TY_CHAR) + call salloc (cleanstr, SZ_FNAME, TY_CHAR) + call salloc (apsex, napsex, TY_POINTER) + + # Select apertures to extract and fix possible limit error. + napsex = 0 + do i = 1, naps { + if (AP_LOW(aps[i],1) > AP_HIGH(aps[i],1)) { + xmax = AP_LOW(aps[i],1) + AP_LOW(aps[i],1) = AP_HIGH(aps[i],1) + AP_HIGH(aps[i],1) = xmax + } + if (AP_LOW(aps[i],2) > AP_HIGH(aps[i],2)) { + xmax = AP_LOW(aps[i],2) + AP_LOW(aps[i],2) = AP_HIGH(aps[i],2) + AP_HIGH(aps[i],2) = xmax + } + if (AP_SELECT(aps[i]) == NO) + next + Memi[apsex+napsex] = aps[i] + napsex = napsex + 1 + } + + # Get CL parameters + bkg = apgwrd ("background", Memc[bkgstr], SZ_FNAME, BACKGROUND) + pfit = apgwrd ("pfit", Memc[str], SZ_LINE, P_FIT) + clean = apgetb ("clean") + if (clean) + call strcpy ("yes", Memc[cleanstr], SZ_FNAME) + else + call strcpy ("no", Memc[cleanstr], SZ_FNAME) + nsubaps = apgeti ("nsubaps") + interptype = II_LINEAR + + # Do clobber checking. Return if output exists and not clobbering. + call apgstr ("ansclobber", Memc[str], SZ_LINE) + call appstr ("ansclobber1", Memc[str]) + fmt = ap_check (input, output, format, Memi[apsex], napsex, nsubaps) + if (fmt == 0) { + call sfree (sp) + return + } + + # Force weights depending on format or cleaning. + switch (fmt) { + case FLAT, RATIO, DIFF, FIT, NOISE: + weights = W_VARIANCE + default: + if (clean) { + call strcpy ("variance", Memc[wtstr], SZ_FNAME) + weights = W_VARIANCE + } else + weights = apgwrd ("weights", Memc[wtstr], SZ_FNAME, WEIGHTS) + } + + if (clgetb ("verbose")) { + call printf ("Extracting apertures ...\n") + call flush (STDOUT) + } + + # Open input image and profile image if given. Set axis parameters + # where 'a' is the aperture axis across the dispersion and 'b' is + # along the dispersion. + + im = ap_immap (input, aaxis, baxis) + namax = IM_LEN(im, aaxis) + nb = IM_LEN(im, baxis) + + pim = NULL + if (strne(profiles,input) && weights==W_VARIANCE && profiles[1]!=EOS) { + pim = ap_immap (profiles, i, j) + if (i!=aaxis||j!=baxis||IM_LEN(pim,i)!=namax||IM_LEN(pim,j)!=nb) { + call imunmap (pim) + call imunmap (im) + call sfree (sp) + call error (1, + "Input image and profile image are not compatible") + } + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Using profile image %s for %s") + call pargstr (profiles) + call pargstr (input) + call ap_log (Memc[str], YES, YES, NO) + } + + # Determine limits of apertures for use in defining memory requirements + # and I/O. + + call salloc (a, 2 * napsex, TY_INT) + call salloc (b, 2 * napsex, TY_INT) + call salloc (c, 2 * napsex, TY_REAL) + a1 = a - 1 + a2 = a1 + napsex + b1 = b - 1 + b2 = b1 + napsex + c1 = c - 1 + c2 = c1 + napsex + + # Initialize image interpolator for edge pixel weighting. + switch (interptype) { + case II_LINEAR: + interpbuf = 2 + case II_POLY3: + interpbuf = 3 + case II_SINC: + interpbuf = 16 + default: + interpbuf = 0 + } + if (interptype > 0) + call asiinit (asi, interptype) + else + asi = NULL + + na1 = 0 + do i = 1, napsex { + ap = Memi[apsex+i-1] + cv = AP_CV(ap) + ic = AP_IC(ap) + + # Dispersion axis limits + bmin = min (nb, max (1, nint (AP_CEN(ap,baxis)+AP_LOW(ap,baxis)))) + bmax = max (1, min (nb, nint (AP_CEN(ap,baxis)+AP_HIGH(ap,baxis)))) + + # Aperture axis shifts + if (cv != NULL) { + cmin = MAX_REAL + cmax = -MAX_REAL + do j = bmin, bmax { + shift = ap_cveval (cv, real (j)) + cmin = min (cmin, shift) + cmax = max (cmax, shift) + } + } else { + cmin = 0. + cmax = 0. + } + + # Background region limits. + xmin = AP_LOW(ap,aaxis) + xmax = AP_HIGH(ap,aaxis) + if (weights == W_VARIANCE) { + xmin = xmin - 2 + xmax = xmax + 2 + } + xmin = xmin - interpbuf + xmax = xmax + interpbuf + if (bkg != B_NONE && AP_IC(ap) != NULL) { + xmin = min (xmin, ic_getr (ic, "xmin")) + xmax = max (xmax, ic_getr (ic, "xmax")) + } + + Memi[a1+i] = min (namax, max (1, nint (AP_CEN(ap,aaxis)+xmin+cmin))) + Memi[a2+i] = max (1, min (namax, nint (AP_CEN(ap,aaxis)+xmax+cmax))) + Memi[b1+i] = bmin + Memi[b2+i] = bmax + Memr[c1+i] = cmin + Memr[c2+i] = cmax + } + call alimi (Memi[a], 2*napsex, amin, amax) + call alimi (Memi[b], 2*napsex, bmin, bmax) + + # The maximum size of the image in memory is 80% of the maximum + # working set size returned by begmem or 40% if a profile image + # is used. Later I/O may exceed this since at least one + # aperture + background is needed in memory. + + new_size = begmem (0, old_size, max_size) + namax = (amax - amin + 1) + nb = (bmax - bmin + 1) + if (pim == NULL) + namax = min (namax, int (0.8 * max_size / SZ_REAL / nb)) + else + namax = min (namax, int (0.8 * max_size / SZ_REAL / nb / 2)) + best_size = 1.2 * namax * nb * SZ_REAL + new_size = begmem (best_size, old_size, max_size) + + # Allocate auxilary memory. Some memory is only dependent on the + # number of dispersion points and subapertures and is the same for + # all apertures. Other memory, such as the sky and profile depend on + # the aperture widths and tilts which may vary. The input data is + # expected to have the aperture axis along the first dimension. If + # the image is in this orientation then the IMIO buffer is used. + # Otherwise sequential I/O is used and transposed into the allocated + # memory. + + iferr { + call salloc (astart, nb, TY_INT) + call salloc (spec, nsubaps * nb, TY_REAL) + if (weights == W_VARIANCE) { + call salloc (raw, nsubaps * nb, TY_REAL) + call salloc (specsig, nsubaps * nb, TY_REAL) + } else { + raw = NULL + specsig = NULL + } + profile = NULL + if (aaxis == 2) { + call calloc (dbuf, namax * nb, TY_REAL) + if (pim != NULL) + call calloc (pbuf, namax * nb, TY_REAL) + } + + # For variance weighting the computations are done in photon units. + if (weights == W_VARIANCE) { + gain = apgimr ("gain", im) + rdnoise = apgimr ("readnoise", im) + } else { + gain = 1 + rdnoise = 0 + } + + # Loop through each aperture doing the extractions. + amax = 0 + do i = 1, napsex { + ap = Memi[apsex+i-1] + + # Check if a new input data buffer is needed. As many apertures + # as possible are read at once within the given memory limits + # though at least one aperture must be read. Do a transpose if + # needed. + + if (Memi[a1+i] < amin || Memi[a2+i] > amax) { + amin = Memi[a1+i] + amax = Memi[a2+i] + do j = i, napsex { + amin = min (amin, Memi[a1+j]) + amax = max (amax, Memi[a2+j]) + na = amax - amin + 1 + if (na > namax) + break + } + + if (aaxis == 1) { + if (fmt == DIFF) { + call mfree (dbuf, TY_REAL) + call malloc (dbuf, na*nb, TY_REAL) + call amovr (Memr[imgs2r(im,amin,amax,bmin,bmax)], + Memr[dbuf], na*nb) + } else + dbuf = imgs2r (im, amin, amax, bmin, bmax) + } else { + if (na > namax) { + call mfree (dbuf, TY_REAL) + namax = na + call calloc (dbuf, namax * nb, TY_REAL) + } + do j = amin, amax { + sbuf = imgl2r (im, j) + sbuf = sbuf + bmin - 1 + ptr = dbuf + j - amin + do k = bmin, bmax { + Memr[ptr] = Memr[sbuf] + sbuf = sbuf + 1 + ptr = ptr + na + } + } + } + if (pim != NULL) { + if (aaxis == 1) + pbuf = imgs2r (pim, amin, amax, bmin, bmax) + else { + if (na > namax) { + call mfree (pbuf, TY_REAL) + namax = na + call calloc (pbuf, namax * nb, TY_REAL) + } + do j = amin, amax { + sbuf = imgl2r (pim, j) + sbuf = sbuf + bmin - 1 + ptr = pbuf + j - amin + do k = bmin, bmax { + Memr[ptr] = Memr[sbuf] + sbuf = sbuf + 1 + ptr = ptr + na + } + } + } + } + if (weights == W_VARIANCE && gain != 1.) { + j = na * nb + call amulkr (Memr[dbuf], gain, Memr[dbuf], j) + if (pim != NULL) + call amulkr (Memr[pbuf], gain, Memr[pbuf], j) + } + } + + # To minimize memory a variable integer offset is used to + # accomodate the aperture tilts. The offsets are stored in + # the astart array and the width of any one line determined. + # If a stored profile is used it is read and it is ASSUMED to + # be valid for the input aperture with the same ID. If no + # stored profile is found the profile fitting algorithm + # parameter determines whether to fit 1D function along the + # image axes (in which case all the profile offsets are the + # same) or if the Marsh algorithm for tilted spectra is + # used. In the latter the offsets can be adjusted to mimize + # memory and a buffer of two pixels around the aperture is + # required by the algorithm. + + if (weights == W_NONE) { + xmin = AP_CEN(ap,aaxis) + AP_LOW(ap,aaxis) + xmax = AP_CEN(ap,aaxis) + AP_HIGH(ap,aaxis) + xmin = xmin - interpbuf + xmax = xmax + interpbuf + na1 = nint (xmax) - nint (xmin) + 1 + cv = AP_CV(ap) + do j = bmin, bmax { + shift = ap_cveval (cv, real (j)) + Memi[astart+j-bmin] = nint (xmin + shift) + } + } else { + if (pfit == P_FIT1D) { + xmin = AP_CEN(ap,aaxis) + AP_LOW(ap,aaxis) + Memr[c1+i] + xmax = AP_CEN(ap,aaxis) + AP_HIGH(ap,aaxis) + Memr[c2+i] + xmin = xmin - interpbuf + xmax = xmax + interpbuf + na1 = nint (xmax) - nint (xmin) + 1 + call amovki (nint (xmin), Memi[astart], nb) + } else if (pfit == P_FIT2D) { + xmin = AP_CEN(ap,aaxis) + AP_LOW(ap,aaxis) - 2 + xmax = AP_CEN(ap,aaxis) + AP_HIGH(ap,aaxis) + 2 + xmin = xmin - interpbuf + xmax = xmax + interpbuf + na1 = nint (xmax) - nint (xmin) + 1 + cv = AP_CV(ap) + do j = bmin, bmax { + shift = ap_cveval (cv, real (j)) + Memi[astart+j-bmin] = nint (xmin + shift) + } + } + } + + # Do the sky or background determination if needed. An array + # of the same size as the 2D aperture is returned as well as + # a single estimate of the variance in the sky value at each + # line based on the fit. If a profile image is used then the + # sky is for the profile image and the object sky is + # determined later in order to reuse the sky buffers. + + if (bkg != B_NONE && AP_IC(ap) != NULL) { + call malloc (sbuf, na1 * nb, TY_REAL) + call malloc (svar, nb, TY_REAL) + call malloc (specsky, nsubaps * nb, TY_REAL) + if (pim == NULL) + call ap_skyeval (im, ap, dbuf, na, nb, amin, 1, + Memr[sbuf], Memr[svar], Memr[specsky], na1, nb, + Memi[astart], 1, nsubaps, rdnoise) + else + call ap_skyeval (pim, ap, pbuf, na, nb, amin, 1, + Memr[sbuf], Memr[svar], Memr[specsky], na1, nb, + Memi[astart], 1, nsubaps, rdnoise) + } else { + sbuf = NULL + svar = NULL + specsky = NULL + } + + # Use a quick sum for unweighted extraction. For weighed + # extractions we use either a previously determined profile + # or call the profile routine. If desired the profile is + # stored for later use. Then the variance weighted + # extraction routine is called. + + if (weights == W_NONE) + call ap_sum (ap, dbuf, na, nb, amin, 1, sbuf, na1, nb, + Memi[astart], 1, Memr[spec], nsubaps, asi) + else { + call malloc (profile, na1 * nb, TY_REAL) + if (pim == NULL) + call ap_profile (im, ap, dbuf, na, nb, amin, 1, sbuf, + svar, Memr[profile], na1, nb, Memi[astart], 1, + asi) + else { + call ap_profile (pim, ap, pbuf, na, nb, amin, 1, sbuf, + svar, Memr[profile], na1, nb, Memi[astart], 1, + asi) + if (sbuf != NULL) + call ap_skyeval (im, ap, dbuf, na, nb, amin, 1, + Memr[sbuf], Memr[svar], Memr[specsky], na1, nb, + Memi[astart], 1, nsubaps, rdnoise) + } + + call ap_variance (im, ap, dbuf, na, nb, amin, 1, sbuf, svar, + Memr[profile], na1, nb, Memi[astart], 1, Memr[spec], + Memr[raw], Memr[specsig], nsubaps, asi) + } + + # Output the extracted spectrum. The extras of sky, sigma, + # and unweighted spectrum may also be stored. If the extra + # information is not available the pointers will be NULL. + + if (weights == W_VARIANCE && gain != 1.) { + call adivkr (Memr[spec], gain, Memr[spec], nb) + if (raw != NULL) + call adivkr (Memr[raw], gain, Memr[raw], nb) + if (specsky != NULL) + call adivkr (Memr[specsky], gain, Memr[specsky], nb) + if (specsig != NULL) + call adivkr (Memr[specsig], gain, Memr[specsig], nb) + call amulkr (Memr[profile], gain, Memr[profile], nb*na1) + } + + call ap_output (input, output, format, Memc[bkgstr], + Memc[wtstr], Memc[cleanstr], gain, im, Memi[apsex], napsex, + i, nsubaps, spec, raw, specsky, specsig, dbuf, na, nb, amin, + 1, sbuf, profile, na1, nb, Memi[astart], 1) + + call mfree (profile, TY_REAL) + call mfree (sbuf, TY_REAL) + call mfree (svar, TY_REAL) + call mfree (specsky, TY_REAL) + } + + # Finish up and restore the working set size. + if (asi != NULL) + call asifree (asi) + if (pim != NULL) { + if (aaxis == 2) + call mfree (pbuf, TY_REAL) + call imunmap (pim) + } + if (aaxis == 2) + call mfree (dbuf, TY_REAL) + call imunmap (im) + call fixmem (old_size) + call sfree (sp) + + } then { + call mfree (profile, TY_REAL) + call mfree (sbuf, TY_REAL) + call mfree (svar, TY_REAL) + call mfree (specsky, TY_REAL) + + if (asi != NULL) + call asifree (asi) + if (pim != NULL) { + if (aaxis == 2) + call mfree (pbuf, TY_REAL) + call imunmap (pim) + } + if (aaxis == 2) + call mfree (dbuf, TY_REAL) + call imunmap (im) + call fixmem (old_size) + call sfree (sp) + + call erract (EA_ERROR) + } +end + + +# AP_CHECK -- Check if output spectra exist. If the user allows clobbering, +# delete the spectra. Return the format. + +int procedure ap_check (input, output, format, aps, naps, nsubaps) + +char input[ARB] # Input image name +char output[ARB] # Output root name +char format[ARB] # Output format +pointer aps[naps] # Apertures +int naps # Number of apertures +int nsubaps # Number of subapertures + +int i, j, fmt +pointer sp, name, name1, input1, ksection, ans + +int strdic(), imaccess(), stridxs() +bool streq(), ap_answer() + +begin + call smark (sp) + call salloc (name, SZ_LINE, TY_CHAR) + call salloc (name1, SZ_LINE, TY_CHAR) + call salloc (input1, SZ_LINE, TY_CHAR) + call salloc (ksection, SZ_LINE, TY_CHAR) + call salloc (ans, SZ_LINE, TY_CHAR) + + fmt = strdic (format, format, SZ_LINE, FORMATS) + call imgimage (input, Memc[input1], SZ_LINE) + + switch (fmt) { + case MULTISPEC, NORM, FLAT, RATIO, DIFF, FIT: + i = stridxs ("[", Memc[input1]) + if (i > 0) { + call strcpy (Memc[input1+i-1], Memc[ksection], SZ_LINE) + Memc[input1+i-1] = EOS + } else + Memc[ksection] = EOS + if (output[1] == EOS) + call strcpy (Memc[input1], Memc[name], SZ_LINE) + else + call strcpy (output, Memc[name], SZ_LINE) + + switch (fmt) { + case MULTISPEC: + if (streq (Memc[input1], Memc[name])) { + call strcat (".ms", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case NORM: + if (streq (Memc[input1], Memc[name])) { + call strcat (".norm", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case FLAT: + if (streq (Memc[input1], Memc[name])) { + call strcat (".flat", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case RATIO: + if (streq (Memc[input1], Memc[name])) { + call strcat (".ratio", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case DIFF: + if (streq (Memc[input1], Memc[name])) { + call strcat (".diff", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case FIT: + if (streq (Memc[input1], Memc[name])) { + call strcat (".fit", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + } + if (imaccess (Memc[name], 0) == YES) { + call sprintf (Memc[ans], SZ_LINE, + "Clobber existing output image %s?") + call pargstr (Memc[name]) + if (ap_answer ("ansclobber1", Memc[ans])) + call imdelete (Memc[name]) + else { + call sprintf (Memc[ans], SZ_LINE, + "EXTRACT - Output spectrum %s already exists") + call pargstr (Memc[name]) + call ap_log (Memc[ans], YES, NO, YES) + fmt = 0 + } + } + case ECHELLE: + if (output[1] == EOS) + call strcpy (Memc[input1], Memc[name], SZ_LINE) + else + call strcpy (output, Memc[name], SZ_LINE) + + do i = 1, nsubaps { + if (nsubaps == 1) + call strcpy (Memc[name], Memc[name1], SZ_LINE) + else { + call sprintf (Memc[name1], SZ_LINE, "%s%0*d") + call pargstr (Memc[name]) + call pargi (int(log10(real(nsubaps)))+1) + call pargi (i) + } + if (streq (Memc[input1], Memc[name])) { + call strcat (".ec", Memc[name1], SZ_LINE) + call strcat (Memc[ksection], Memc[name1], SZ_LINE) + } + + if (imaccess (Memc[name1], 0) == YES) { + call sprintf (Memc[ans], SZ_LINE, + "Clobber existing output image %s?") + call pargstr (Memc[name1]) + if (ap_answer ("ansclobber1", Memc[ans])) + call imdelete (Memc[name1]) + else { + call sprintf (Memc[ans], SZ_LINE, + "EXTRACT - Output spectrum %s already exists") + call pargstr (Memc[name1]) + call ap_log (Memc[ans], YES, NO, YES) + fmt = 0 + } + } + } + case ONEDSPEC, STRIP: + do i = 1, naps { + do j = 1, nsubaps { + call sprintf (Memc[name], SZ_LINE, "%s.%0*d") + if (output[1] == EOS) + call pargstr (Memc[input1]) + else + call pargstr (output) + call pargi (int(log10(real(nsubaps)))+4) + call pargi (AP_ID(aps[i])+(j-1)*1000) + if (imaccess (Memc[name], 0) == YES) { + call sprintf (Memc[ans], SZ_LINE, + "Clobber existing output image %s?") + call pargstr (Memc[name]) + if (ap_answer ("ansclobber1", Memc[ans])) + call imdelete (Memc[name]) + else { + call sprintf (Memc[ans], SZ_LINE, + "EXTRACT - Output spectrum %s already exists") + call pargstr (Memc[name]) + call ap_log (Memc[ans], YES, NO, YES) + fmt = 0 + } + } + } + } + case NOISE: + ; + default: + call sfree (sp) + call error (1, "EXTRACT - Unknown output format") + } + + call sfree (sp) + return (fmt) +end + + +# AP_OUTPUT -- Review the extracted spectra and write them to an image. +# This routine determines the output format and whether to also output sky +# unweighted, and sigma spectra. The appropriate header keywords have +# to be added. + +procedure ap_output (image, output, format, bkg, wt, clean, gain, in, aps, + naps, iap, nsubaps, spec, raw sky, sig, dbuf, nc, nl, c1, l1, sbuf, + profile, nx, ny, xs, ys) + +char image[ARB] # Input image name +char output[ARB] # Output root name +char format[ARB] # Output format +char bkg[ARB] # Background type +char wt[ARB] # Weight type +char clean[ARB] # Clean? +real gain # Gain +pointer in # Input IMIO pointer +pointer aps[naps] # Apertures +int naps # Number of apertures +int iap # Aperture +int nsubaps # Number of subapertures +pointer spec # Output spectrum +pointer raw # Output raw spectrum +pointer sky # Output sky +pointer sig # Output sigma +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +pointer sbuf # Sky values (NULL if none) +pointer profile # Profile (NULL if none) +int nx, ny # Size of sky and profile array +int xs[ny], ys # Origin of sky and profile array + +int fmt # Output format +bool extras # Include raw spectrum, sky, and sigma + +real low, high, step +int i, k, l, m, apid, apaxis, dispaxis +pointer sp, str, str1, name, name1, input, ksection +pointer ap, out, outsave, gt, apmw, buf +pointer sum2, sum4, nsum + +real clgetr() +int scan(), strdic(), imaccf(), stridxs() +bool streq(), ap_answer(), apgetb() +pointer immap(), imgl2r(), impl2r(), impl3r() +pointer gt_init(), apmw_open() +errchk immap, impl2r, impl3r, imps2r, ap_strip, ap_pstrip, apmw_open +errchk ap_fitspec, ap_lnorm, ap_cnorm, ap_lflat, ap_cflat + +begin + # Allocate string and file name arrays. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (name, SZ_LINE, TY_CHAR) + call salloc (name1, SZ_LINE, TY_CHAR) + call salloc (input, SZ_LINE, TY_CHAR) + call salloc (ksection, SZ_LINE, TY_CHAR) + + fmt = strdic (format, format, SZ_LINE, FORMATS) + extras = apgetb ("extras") + + ap = aps[iap] + apaxis = AP_AXIS(ap) + dispaxis = mod (apaxis, 2) + 1 + + # Set output name. + call imgimage (image, Memc[input], SZ_LINE) + i = stridxs ("[", Memc[input]) + if (i > 0) { + call strcpy (Memc[input+i-1], Memc[ksection], SZ_LINE) + Memc[input+i-1] = EOS + i = stridxs ("]", Memc[ksection]) + call strcpy (",append]", Memc[ksection+i-1], SZ_LINE) + } else + Memc[ksection] = EOS + if (output[1] == EOS) + call strcpy (Memc[input], Memc[name], SZ_LINE) + else + call strcpy (output, Memc[name], SZ_LINE) + + switch (fmt) { + case ECHELLE: + ; + case MULTISPEC: + if (streq (Memc[input], Memc[name])) { + call strcat (".ms", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case NORM: + if (streq (Memc[input], Memc[name])) { + call strcat (".norm", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case FLAT: + if (streq (Memc[input], Memc[name])) { + call strcat (".flat", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case RATIO: + if (streq (Memc[input], Memc[name])) { + call strcat (".ratio", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case DIFF: + if (streq (Memc[input], Memc[name])) { + call strcat (".diff", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case FIT: + if (streq (Memc[input], Memc[name])) { + call strcat (".fit", Memc[name], SZ_LINE) + call strcat (Memc[ksection], Memc[name], SZ_LINE) + } + case NOISE: + Memc[name] = EOS + } + + + # Set the review graph title. + call sprintf (Memc[str], SZ_LINE, "%s: %s - Aperture %s") + call pargstr (image) + call pargstr (IM_TITLE(in)) + call pargi (AP_ID(ap)) + + gt = gt_init () + call gt_sets (gt, GTTITLE, Memc[str]) + + # Query the user whether to review the extraction. + call sprintf (Memc[str], SZ_LINE, + "Review extracted spectrum for aperture %d from %s?") + call pargi (AP_ID(ap)) + call pargstr (image) + + # If reviewing graph the spectrum, do a cursor loop, and allow + # the user to skip the output or define a new output image. + if (ap_answer ("ansreview1", Memc[str])) { + call ap_graph1 (gt, Memr[spec], ny, nsubaps) + + if (fmt == ONEDSPEC && nsubaps == 1) { + call printf ( + "Output image name [use # to skip output] (%s): ") + call pargstr (Memc[name]) + call flush (STDOUT) + if (scan() != EOF) { + call gargwrd (Memc[str], SZ_LINE) + if (Memc[str] == '#') { + call gt_free (gt) + call sfree (sp) + return + } + if (Memc[str] != EOS) + call strcpy (Memc[str], Memc[name], SZ_LINE) + } + } + } + + # Output the image. + switch (fmt) { + case MULTISPEC: + if (iap == 1) { + out = immap (Memc[name], NEW_COPY, in) + + IM_PIXTYPE(out) = TY_REAL + IM_NDIM(out) = 1 + IM_LEN(out, 1) = ny + IM_LEN(out, 2) = nsubaps * naps + IM_LEN(out, 3) = 1 + if (extras) { + if (sky != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + if (raw != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + if (sig != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + } + if (IM_LEN(out, 2) > 1) + IM_NDIM(out) = 2 + if (IM_LEN(out, 3) > 1) + IM_NDIM(out) = 3 + + apmw = apmw_open (in, out, dispaxis, nsubaps*naps, ny) + + # Write BAND IDs. + k = 1 + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "spectrum - background %s, weights %s, clean %s") + call pargstr (bkg) + call pargstr (wt) + call pargstr (clean) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + if (extras) { + if (raw != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "raw - background %s, weights none, clean no") + call pargstr (bkg) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + } + if (sky != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "background - background %s") + call pargstr (bkg) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + } + if (sig != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "sigma - background %s, weights %s, clean %s") + call pargstr (bkg) + call pargstr (wt) + call pargstr (clean) + call imastr (out, Memc[str1], Memc[str]) + } + } + + do k = 1, naps { + low = AP_CEN(aps[k],apaxis) + AP_LOW(aps[k],apaxis) + high = AP_CEN(aps[k],apaxis) + AP_HIGH(aps[k],apaxis) + step = (high - low) / nsubaps + low = low - step + do l = 1, nsubaps { + apid = AP_ID(aps[k]) + (l - 1) * 1000 + low = low + step + high = low + step + call apmw_setap (apmw, (k-1)*nsubaps+l, + apid, AP_BEAM(aps[k]), low, high) + } + } + do k = 1, naps { + if (AP_TITLE(aps[k]) != NULL) { + do l = 1, nsubaps { + call sprintf (Memc[str], SZ_LINE, "APID%d") + call pargi ((k-1)*nsubaps+l) + call imastr (out, Memc[str], + Memc[AP_TITLE(aps[k])]) + } + } + } + } + + do l = 1, nsubaps { + k = (iap - 1) * nsubaps + l + buf = impl2r (out, k) + call amovr (Memr[spec+(l-1)*ny], Memr[buf], ny) + if (extras) { + m = 2 + if (raw != NULL) { + buf = impl3r (out, k, m) + call amovr (Memr[raw+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + if (sky != NULL) { + buf = impl3r (out, k, m) + call amovr (Memr[sky+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + if (sig != NULL) { + buf = impl3r (out, k, m) + call amovr (Memr[sig+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + } + } + if (iap == naps) { + call apmw_saveim (apmw, out, fmt) + call apmw_close (apmw) + call imunmap (out) + } + + if (Memc[name] != EOS) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s") + call pargi (AP_ID(ap)) + call pargstr (image) + call pargstr (Memc[name]) + call ap_log (Memc[str], YES, YES, NO) + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + } + + case ECHELLE: + do l = 1, nsubaps { + if (nsubaps == 1) + call strcpy (Memc[name], Memc[name1], SZ_LINE) + else { + call sprintf (Memc[name1], SZ_LINE, "%s%0*d") + call pargstr (Memc[name]) + call pargi (int(log10(real(nsubaps)))+1) + call pargi (l) + } + if (streq (Memc[input], Memc[name])) { + call strcat (".ec", Memc[name1], SZ_LINE) + call strcat (Memc[ksection], Memc[name1], SZ_LINE) + } + + if (iap == 1) { + out = immap (Memc[name1], NEW_COPY, in) + + IM_PIXTYPE(out) = TY_REAL + IM_NDIM(out) = 1 + IM_LEN(out, 1) = ny + IM_LEN(out, 2) = naps + IM_LEN(out, 3) = 1 + if (extras) { + if (sky != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + if (raw != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + if (sig != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + } + if (IM_LEN(out, 2) > 1) + IM_NDIM(out) = 2 + if (IM_LEN(out, 3) > 1) + IM_NDIM(out) = 3 + + apmw = apmw_open (in, out, dispaxis, naps, ny) + + # Write BAND IDs. + k = 1 + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "spectrum - background %s, weights %s, clean %s") + call pargstr (bkg) + call pargstr (wt) + call pargstr (clean) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + if (extras) { + if (raw != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "raw - background %s, weights none, clean no") + call pargstr (bkg) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + } + if (sky != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "background - background %s") + call pargstr (bkg) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + } + if (sig != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "sigma - background %s, weights %s, clean %s") + call pargstr (bkg) + call pargstr (wt) + call pargstr (clean) + call imastr (out, Memc[str1], Memc[str]) + } + } + + # Write keyword to allow matching by subaperture. + if (nsubaps > 1) + call imaddi (out, "SUBAP", l) + + do k = 1, naps { + low = AP_CEN(aps[k],apaxis) + AP_LOW(aps[k],apaxis) + high = AP_CEN(aps[k],apaxis) + AP_HIGH(aps[k],apaxis) + step = (high - low) / nsubaps + call apmw_setap (apmw, k, AP_ID(aps[k]), + AP_BEAM(aps[k]), low+(l-1)*step, low+l*step) + } + do k = 1, naps { + if (AP_TITLE(aps[k]) != NULL) { + call sprintf (Memc[str], SZ_LINE, "APID%d") + call pargi (k) + call imastr (out, Memc[str], + Memc[AP_TITLE(aps[k])]) + } + } + } else { + if (l == 1) + out = outsave + else + out = immap (Memc[name1], READ_WRITE, 0) + } + + k = iap + buf = impl2r (out, k) + call amovr (Memr[spec+(l-1)*ny], Memr[buf], ny) + if (extras) { + m = 2 + if (raw != NULL) { + buf = impl3r (out, k, m) + call amovr (Memr[raw+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + if (sky != NULL) { + buf = impl3r (out, k, m) + call amovr (Memr[sky+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + if (sig != NULL) { + buf = impl3r (out, k, m) + call amovr (Memr[sig+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + } + + if (iap == 1) { + call apmw_saveim (apmw, out, fmt) + call apmw_close (apmw) + } + if (l != 1 || iap == naps) + call imunmap (out) + if (l == 1) + outsave = out + + if (nsubaps == 1) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s") + call pargi (AP_ID(ap)) + call pargstr (image) + call pargstr (Memc[name1]) + } else { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d-%d from %s --> %s") + call pargi (AP_ID(ap)) + call pargi (l) + call pargstr (image) + call pargstr (Memc[name1]) + } + call ap_log (Memc[str], YES, YES, NO) + } + + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + + case ONEDSPEC: + do l = 1, nsubaps { + apid = AP_ID(ap) + (l - 1) * 1000 + low = AP_CEN(ap,apaxis) + AP_LOW(ap,apaxis) + high = AP_CEN(ap,apaxis) + AP_HIGH(ap,apaxis) + step = (high - low) / nsubaps + low = low + (l - 1) * step + high = low + step + + call sprintf (Memc[str], SZ_LINE, "%s.%0*d") + call pargstr (Memc[name]) + call pargi (int(log10(real(nsubaps)))+4) + call pargi (apid) + out = immap (Memc[str], NEW_COPY, in) + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s.%0*d") + call pargi (apid) + call pargstr (image) + call pargstr (Memc[name]) + call pargi (int(log10(real(nsubaps)))+4) + call pargi (apid) + call ap_log (Memc[str], YES, YES, NO) + + apmw = apmw_open (in, out, dispaxis, 1, ny) + call apmw_setap (apmw, 1, apid, AP_BEAM(ap), low, high) + if (AP_TITLE(ap) != NULL) + call imastr (out, "APID1", Memc[AP_TITLE(ap)]) + + IM_PIXTYPE(out) = TY_REAL + IM_NDIM(out) = 1 + IM_LEN(out, 1) = ny + IM_LEN(out, 2) = 1 + IM_LEN(out, 3) = 1 + if (extras) { + if (sky != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + if (raw != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + if (sig != NULL) + IM_LEN(out, 3) = IM_LEN(out, 3) + 1 + } + if (IM_LEN(out, 2) > 1) + IM_NDIM(out) = 2 + if (IM_LEN(out, 3) > 1) + IM_NDIM(out) = 3 + + # Write BAND IDs. + k = 1 + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "spectrum: background %s, weights %s, clean %s") + call pargstr (bkg) + call pargstr (wt) + call pargstr (clean) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + if (extras) { + if (raw != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "spectrum: background %s, weights none, clean no") + call pargstr (bkg) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + } + if (sky != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "background: background %s") + call pargstr (bkg) + call imastr (out, Memc[str1], Memc[str]) + k = k + 1 + } + if (sig != NULL) { + call sprintf (Memc[str1], SZ_LINE, "BANDID%d") + call pargi (k) + call sprintf (Memc[str], SZ_LINE, + "sigma - background %s, weights %s, clean %s") + call pargstr (bkg) + call pargstr (wt) + call pargstr (clean) + call imastr (out, Memc[str1], Memc[str]) + } + } + + buf = impl2r (out, 1) + call amovr (Memr[spec+(l-1)*ny], Memr[buf], ny) + if (extras) { + m = 2 + if (raw != NULL) { + buf = impl3r (out, 1, m) + call amovr (Memr[raw+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + if (sky != NULL) { + buf = impl3r (out, 1, m) + call amovr (Memr[sky+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + if (sig != NULL) { + buf = impl3r (out, 1, m) + call amovr (Memr[sig+(l-1)*ny], Memr[buf], ny) + m = m + 1 + } + } + + call apmw_saveim (apmw, out, fmt) + call apmw_close (apmw) + call imunmap (out) + + } + + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + + case STRIP: + do l = 1, nsubaps { + apid = AP_ID(ap) + (l - 1) * 1000 + low = AP_CEN(ap,apaxis) + AP_LOW(ap,apaxis) + high = AP_CEN(ap,apaxis) + AP_HIGH(ap,apaxis) + step = (high - low) / nsubaps + low = low + (l - 1) * step + high = low + step + + call sprintf (Memc[str], SZ_LINE, "%s.%0*d") + call pargstr (Memc[name]) + call pargi (int(log10(real(nsubaps)))+4) + call pargi (apid) + out = immap (Memc[str], NEW_COPY, in) + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s.%0*d") + call pargi (apid) + call pargstr (image) + call pargstr (Memc[name]) + call pargi (int(log10(real(nsubaps)))+4) + call pargi (apid) + call ap_log (Memc[str], YES, YES, NO) + + apmw = apmw_open (in, out, dispaxis, 1, ny) + call apmw_setap (apmw, 1, apid, AP_BEAM(ap), low, high) + call sprintf (Memc[str], SZ_LINE, "%s - Aperture %d") + call pargstr (IM_TITLE(out)) + call pargi (AP_ID(ap)) + call strcpy (Memc[str], IM_TITLE(out), SZ_IMTITLE) + if (AP_TITLE(ap) != NULL) + call imastr (out, "APID1", Memc[AP_TITLE(ap)]) + + IM_PIXTYPE(out) = TY_REAL + IM_NDIM(out) = 2 + IM_LEN(out, 1) = ny + IM_LEN(out, 2) = high - low + 1 + + if (profile == NULL) + call ap_strip (ap, low, high, out, dbuf, nc, nl, c1, l1, + sbuf, nx, ny, xs, ys) + else + call ap_pstrip (ap, low, high, out, gain, Memr[spec], + Memr[profile], nx, ny, xs, ys) + + call apmw_saveim (apmw, out, fmt) + call apmw_close (apmw) + call imunmap (out) + } + + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + + case NORM, FLAT: + if (iap == 1) { + out = immap (Memc[name], NEW_COPY, in) + IM_PIXTYPE(out) = TY_REAL + if (imaccf (out, "CCDMEAN") == YES) + call imdelf (out, "CCDMEAN") + call ap_fitspec (ap, in, Memr[spec], ny) + k = YES + } else { + call ap_fitspec (ap, in, Memr[spec], ny) + k = NO + } + if (apaxis == 1) { + if (fmt == NORM) + call ap_lnorm (ap, out, gain, dbuf, nc, nl, c1, l1, + Memr[spec], ny, ys, k) + else + call ap_lflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec], + sbuf, Memr[profile], nx, ny, xs, ys, k) + } else { + if (fmt == NORM) + call ap_cnorm (ap, out, gain, dbuf, nc, nl, c1, l1, + Memr[spec], ny, ys, k) + else + call ap_cflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec], + sbuf, Memr[profile], nx, ny, xs, ys, k) + } + if (iap == naps) + call imunmap (out) + + if (Memc[name] != EOS) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s") + call pargi (AP_ID(ap)) + call pargstr (image) + call pargstr (Memc[name]) + call ap_log (Memc[str], YES, YES, NO) + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + } + + case RATIO, FIT: + if (iap == 1) { + out = immap (Memc[name], NEW_COPY, in) + IM_PIXTYPE(out) = TY_REAL + k = YES + } else + k = NO + if (apaxis == 1) { + switch (fmt) { + case RATIO: + call ap_lflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec], + sbuf, Memr[profile], nx, ny, xs, ys, k) + case FIT: + call ap_lfit (ap, out, gain, Memr[spec], Memr[profile], + nx, ny, xs, ys, k) + } + } else { + switch (fmt) { + case RATIO: + call ap_cflat (ap, out, dbuf, nc, nl, c1, l1, Memr[spec], + sbuf, Memr[profile], nx, ny, xs, ys, k) + case FIT: + call ap_cfit (ap, out, gain, Memr[spec], Memr[profile], + nx, ny, xs, ys, k) + } + } + if (iap == naps) + call imunmap (out) + + if (Memc[name] != EOS) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s") + call pargi (AP_ID(ap)) + call pargstr (image) + call pargstr (Memc[name]) + call ap_log (Memc[str], YES, YES, NO) + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + } + + case DIFF: + if (iap == 1) { + out = immap (Memc[name], NEW_COPY, in) + IM_PIXTYPE(out) = TY_REAL + do k = 1, IM_LEN(in,2) { + buf = impl2r (out, k) + call amovr (Memr[imgl2r(in,k)], Memr[buf], IM_LEN(out,1)) + } + k = NO + } else + k = NO + if (apaxis == 1) + call ap_ldiff (ap, out, gain, dbuf, nc, nl, c1, l1, Memr[spec], + Memr[profile], nx, ny, xs, ys, k) + else + call ap_cdiff (ap, out, gain, dbuf, nc, nl, c1, l1, Memr[spec], + Memr[profile], nx, ny, xs, ys, k) + if (iap == naps) + call imunmap (out) + + if (Memc[name] != EOS) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s") + call pargi (AP_ID(ap)) + call pargstr (image) + call pargstr (Memc[name]) + call ap_log (Memc[str], YES, YES, NO) + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + } + + case NOISE: + if (iap == 1) { + low = clgetr ("dmin") + high = clgetr ("dmax") + l = clgetr ("nbins") + if (high < low) { + step = low; low = high; high = step + } + step = (high - low) / l + call malloc (sum2, l, TY_REAL) + call malloc (sum4, l, TY_REAL) + call malloc (nsum, l, TY_INT) + call aclrr (Memr[sum2], l) + call aclrr (Memr[sum4], l) + call aclri (Memi[nsum], l) + } + call ap_noise (ap, gain, dbuf, nc, nl, c1, l1, sbuf, Memr[spec], + Memr[profile], nx, ny, xs, ys, Memr[sum2], Memr[sum4], + Memi[nsum], l, low, high) + if (iap == naps) { + do k = 0, l-1 { + m = Memi[nsum+k] + if (m > 10) { + Memr[sum2+k] = sqrt (Memr[sum2+k] / (m - 1)) + step = max (0., Memr[sum4+k] / m - Memr[sum2+k]**2) + Memr[sum4+k] = sqrt (sqrt (step / m)) + } else { + Memr[sum2+k] = 0. + Memr[sum4+k] = 0. + } + } + call ap_nplot (image, in, Memr[sum2], Memr[sum4], l, + low, high) + call mfree (sum2, TY_REAL) + call mfree (sum4, TY_REAL) + call mfree (nsum, TY_INT) + } + + if (Memc[name] != EOS) { + call sprintf (Memc[str], SZ_LINE, + "EXTRACT - Aperture %d from %s --> %s") + call pargi (AP_ID(ap)) + call pargstr (image) + call pargstr (Memc[name]) + call ap_log (Memc[str], YES, YES, NO) + call ap_plot1 (gt, Memr[spec], ny, nsubaps) + } + } + + call gt_free (gt) + call sfree (sp) +end + + +# AP_SUM -- Simple, unweighted aperture sum. + +procedure ap_sum (ap, dbuf, nc, nl, c1, l1, sbuf, nx, ny, xs, ys, spec, + nsubaps, asi) + +pointer ap # Aperture structure +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +pointer sbuf # Sky values (NULL if none) +int nx, ny # Size of profile array +int xs[ny], ys # Origin of sky array +real spec[ny, nsubaps] # Spectrum +int nsubaps # Number of subapertures +pointer asi # Interpolator for edge pixel weighting + +int i, ix, iy, ix1, ix2 +real low, high, step, x1, x2, wt1, wt2, s, sval, skyval +real ap_cveval() +pointer cv, data, sky +errchk asifit + +begin + i = AP_AXIS(ap) + low = AP_CEN(ap,i) + AP_LOW(ap,i) + high = AP_CEN(ap,i) + AP_HIGH(ap,i) + step = (high - low) / nsubaps + cv = AP_CV(ap) + do iy = 1, ny { + s = ap_cveval (cv, real (iy + ys - 1)) - c1 + 1 + call ap_asifit (dbuf+(iy+ys-1-l1)*nc, nc, xs[iy]-c1+1, + low+s, high+s, data, asi) +# data = dbuf + (iy + ys - 1 - l1) * nc + xs[iy] - c1 - 1 +# if (asi != NULL) +# call asifit (asi, Memr[data], nc-xs[iy]+c1) + do i = 1, nsubaps { + x1 = max (0.5, low + (i - 1) * step + s) + c1 - xs[iy] + x2 = min (nc + 0.49, low + i * step + s) + c1 - xs[iy] + if (x2 <= x1) { + spec[iy,i] = 0. + next + } + ix1 = nint (x1) + ix2 = nint (x2) + + # Compute end pixel weights. Remember asi is offset by 1. + call ap_edge (asi, x1+1, x2+1, wt1, wt2) + + # Sum pixels. + sval = wt1 * Memr[data+ix1] + wt2 * Memr[data+ix2] + do ix = ix1+1, ix2-1 + sval = sval + Memr[data+ix] + + # Subtract sky if desired. + if (sbuf != NULL) { + sky = sbuf + (iy - 1) * nx - 1 + skyval = wt1 * Memr[sky+ix1] + wt2 * Memr[sky+ix2] + do ix = ix1+1, ix2-1 + skyval = skyval + Memr[sky+ix] + sval = sval - skyval + } + + # Save extracted pixel value. + spec[iy,i] = sval + } + } +end + + +# AP_EDGE -- Compute edge weights. + +procedure ap_edge (asi, x1, x2, wt1, wt2) + +pointer asi #I Image interpolator pointer +real x1, x2 #I Aperture edges +real wt1, wt2 #I Weights + +int ix1, ix2 +real a, b +real asieval(), asigrl() + +begin + # Edge pixel centers. + ix1 = nint (x1) + ix2 = nint (x2) + + # Default weights are fractions of pixel. + if (ix1 == ix2) { + wt1 = (x2 - x1) + wt2 = 0 + } else { + wt1 = (ix1 - x1 + 0.5) + wt2 = (x2 - ix2 + 0.5) + } + + # If there is an interpolator compute fraction of integral. + # We require that data and integrals be positive. + if (asi != NULL) { + if (asieval (asi, real(ix1)) > 0) { + b = asigrl (asi, ix1-0.5, ix1+0.5) + if (b > 0) { + if (ix1 == ix2) + a = asigrl (asi, x1, x2) + else + a = asigrl (asi, x1, ix1+0.5) + if (a > 0 && a < b) + wt1 = a / b + } + } + if (ix1 != ix2 && asieval (asi, real(ix2)) > 0) { + b = asigrl (asi, ix2-0.5, ix2+0.5) + if (b > 0) { + a = asigrl (asi, ix2-0.5, x2) + if (a > 0 && a < b) + wt2 = a / b + } + } + } +end + + +# AP_STRIP -- Simple, unweighted aperture strip. +# Interpolate so that the lower edge of the aperture is the first pixel. + +procedure ap_strip (ap, aplow, aphigh, out, dbuf, nc, nl, c1, l1, sbuf, nx, ny, + xs, ys) + +pointer ap # Aperture structure +real aplow, aphigh # Aperture limits +pointer out # Output IMIO pointer +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +pointer sbuf # Sky values (NULL if none) +int nx, ny # Size of profile array +int xs[ny], ys # Origin of sky array + +int i, na, iy, ix1, ix2, nasi +real low, high, s, x, ap_cveval(), asieval() +pointer obuf, cv, asi, data, sky, ptr, imps2r() + +begin + i = AP_AXIS(ap) + low = aplow - c1 + 1 + high = aphigh - c1 + 1 + cv = AP_CV(ap) + call asiinit (asi, II_LINEAR) + + na = IM_LEN(out,2) + obuf = imps2r (out, 1, ny, 1, na) + call aclrr (Memr[obuf], na * ny) + + do iy = 1, ny { + i = iy + ys - 1 + s = ap_cveval (cv, real (i)) + ix1 = max (1, nint (low + s) - 1) + ix2 = min (nc, nint (high + s) + 1) + nasi = ix2 - ix1 + 1 + if (nasi < 3) + next + data = dbuf + (i - l1) * nc + ix1 - 1 + iferr (call asifit (asi, Memr[data], nasi)) + next + + x = low + s - ix1 + 1 + ptr = obuf + iy - 1 + if (sbuf == NULL) { + do i = 1, na { + if (x >= 1 && x <= nasi) + Memr[ptr] = asieval (asi, x) + x = x + 1. + ptr = ptr + ny + } + } else { + sky = sbuf + (iy - 1) * nx + nint (low + s) - xs[iy] + c1 - 2 + do i = 1, na { + if (x >= 1 && x <= nasi) + Memr[ptr] = asieval (asi, x) - Memr[sky+i] + x = x + 1. + ptr = ptr + ny + } + } + } + + call asifree (asi) +end + + +# AP_PSTRIP -- Profile based strip. +# Interpolate the profile spectrum so that the lower aperture edge is the +# first pixel. + +procedure ap_pstrip (ap, aplow, aphigh, out, gain, spec, profile, nx, ny, + xs, ys) + +pointer ap # Aperture structure +real aplow, aphigh # Aperture limits +pointer out # Output IMIO pointer +real gain # Gain +real spec[ny] # Spectrum +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Origin of profile array + +int na, ix, iy +real low, high, s, x, ap_cveval(), asieval() +pointer sp, cv, asi, data, impl2r() + +begin + call smark (sp) + call salloc (data, nx, TY_REAL) + + ix = AP_AXIS(ap) + low = aplow + high = aphigh + cv = AP_CV(ap) + na = IM_LEN(out,2) + call asiinit (asi, II_LINEAR) + + do iy = 1, ny { + s = spec[iy] / gain + do ix = 1, nx + Memr[data+ix-1] = s * profile[iy,ix] + call asifit (asi, Memr[data], nx) + s = ap_cveval (cv, real (iy+ys-1)) - xs[iy] + 1 + x = low + s + do ix = 1, na { + profile[iy,ix] = asieval (asi, x) + x = x + 1 + } + } + + do ix = 1, na + call amovr (profile[1,ix], Memr[impl2r(out,ix)], ny) + + call asifree (asi) +end + + +# AP_ASIFIT -- Return interpolation pointer and data pointer. +# +# The main reason for this routine is to shift the origin of the data by +# one pixel so that the interpolator may be called to evaluate across +# the extent of the first and last pixels. This means the calling program +# will reference asi fit between 1.5 and N+1.5. It also means the returned +# data pointer may start before the first point but will never be +# dereferenced outside of the data range. + +procedure ap_asifit (dbuf, nc, xs, low, high, data, asi) + +pointer dbuf #I Data buffer pointer +int nc #I Size of data buffer +int xs #I Start of aperture array (in dbuf coords) +real low #I Low aperture edge (in dbuf coords) +real high #I High aperture edge (in dbuf coords) +pointer data #O Data pointer +pointer asi #I ASI pointer + +int i, ix1, ix2, n +real x1, x2 +pointer fit + +begin + # Check for in bounds data. + x1 = max (0.5, low) + x2 = min (nc + 0.49, high) + if (x1 >= x2) + return + + # Set data pointer relative to the aperture start with an offset for + # one indexing; i.e. pixel i is referenced as Memr[data+i]. The + # aperture start may put this outside the data buffer but we expect + # routines using the pointer to never index outside of the buffer. + + data = (dbuf + xs - 1) - 1 + + # If not using an interpolator we are done. + + if (asi == NULL) + return + + # If the aperture, with one extra pixel on each end for integration + # across the end pixel, is within the data buffer then fit an + # interpolator directly. Otherwise we need to use a temporary + # padded buffer. The origin of the fitted buffer is relative + # to the data pointer. Note that this means that evaluating the + # fit requires the aperture start coordinates to be incremented + # by 1. + + ix1 = 0 + ix2 = nint (x2) + 1 - (xs - 1) + n = ix2 + ix1 + 1 + if (data + ix1 >= dbuf && data + ix2 <= dbuf + nc - 1) { + call asifit (asi, Memr[data+ix1], n) + return + } + + # One or the other end point is out of bounds so to avoid potential + # NAN and segmentation errors use an internal array to pad. + + call malloc (fit, n, TY_REAL) + do i = 0, n-1 { + if (data + i < dbuf) + Memr[fit+i] = Memr[dbuf] + else if (data + i > dbuf + nc - 1) + Memr[fit+i] = Memr[dbuf+nc-1] + else + Memr[fit+i] = Memr[data+i] + } + call asifit (asi, Memr[fit], n) + call mfree (fit, TY_REAL) +end |