aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/smw/shdr.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/onedspec/smw/shdr.x')
-rw-r--r--noao/onedspec/smw/shdr.x1269
1 files changed, 1269 insertions, 0 deletions
diff --git a/noao/onedspec/smw/shdr.x b/noao/onedspec/smw/shdr.x
new file mode 100644
index 00000000..bdcc4b95
--- /dev/null
+++ b/noao/onedspec/smw/shdr.x
@@ -0,0 +1,1269 @@
+include <error.h>
+include <imhdr.h>
+include <imset.h>
+include <smw.h>
+include <units.h>
+include <funits.h>
+include <math/iminterp.h>
+
+
+# SHDR_OPEN -- Open the SHDR spectrum header structure.
+# SHDR_TYPE -- Determine spectrum type.
+# SHDR_GTYPE -- Get the selected spectrum type.
+# SHDR_CLOSE -- Close and free the SHDR structure.
+# SHDR_COPY -- Make a copy of an SHDR structure.
+# SHDR_SYSTEM -- Set or change the WCS system.
+# SHDR_UNITS -- Set or change user units.
+# SHDR_LW -- Logical to world coordinate transformation.
+# SHDR_WL -- World to logical coordinate transformation.
+# SHDR_REBIN -- Rebin spectrum to dispersion of reference spectrum.
+# SHDR_LINEAR -- Rebin spectrum to linear dispersion.
+# SHDR_EXTRACT -- Extract a specific wavelength region.
+# SHDR_GI -- Load an integer value from the header.
+# SHDR_GR -- Load a real value from the header.
+
+
+# SHDR_OPEN -- Open SHDR spectrum header structure.
+#
+# This routine sets header information, WCS transformations, and extracts the
+# spectrum from EQUISPEC, MULTISPEC, and NDSPEC format images. The spectrum
+# from a 2D/3D format is specified by a logical line and band number.
+# Optionally an EQUISPEC or MULTISPEC spectrum may be selected by it's
+# aperture number. The access modes are header only or header and data.
+# Special checks are made to avoid repeated setting of the header and WCS
+# information common to all spectra in an image provided the previously set
+# structure is input. Note that the logical to world and world to logical
+# transformations require that the MWCS pointer not be closed.
+
+procedure shdr_open (im, smw, index1, index2, ap, mode, sh)
+
+pointer im # IMIO pointer
+pointer smw # SMW pointer
+int index1 # Image index desired
+int index2 # Image index desired
+int ap # Aperture number desired
+int mode # Access mode
+pointer sh # SHDR pointer
+
+int i, j, k, l, n, np, np1, np2, aplow[2], aphigh[2], strncmp()
+real smw_c1tranr(), asumr()
+double dval, shdr_lw()
+bool newim, streq()
+pointer sp, key, str, coeff, mw, ct, buf
+pointer smw_sctran(), imgs3r(), un_open(), fun_open()
+errchk smw_sctran, imgstr, imgeti, imgetr, smw_gwattrs
+errchk un_open, fun_open, fun_ctranr, imgs3r, shdr_gtype
+
+define data_ 90
+
+begin
+ call smark (sp)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Allocate basic structure or check if the same spectrum is requested.
+ if (sh == NULL) {
+ call calloc (sh, LEN_SHDR, TY_STRUCT)
+ call calloc (SID(sh,1), LEN_SHDRS, TY_CHAR)
+ newim = true
+ } else {
+ call imstats (im, IM_IMAGENAME, Memc[str], SZ_LINE)
+ newim = !streq (Memc[str], IMNAME(sh))
+ if (!newim) {
+ if (LINDEX(sh,1)==index1 && LINDEX(sh,2)==index2) {
+ if (IS_INDEFI(ap) || AP(sh)==ap) {
+ np1 = NP1(sh)
+ np2 = NP2(sh)
+ np = np2 - np1 + 1
+ if (CTLW(sh) == NULL || CTWL(sh) == NULL)
+ goto data_
+ if (mode == SHHDR) {
+ do i = 1, SH_NTYPES
+ call mfree (SPEC(sh,i), TY_REAL)
+ } else {
+ switch (SMW_FORMAT(smw)) {
+ case SMW_ND:
+ if (mode == SHDATA && SPEC(sh,mode) == NULL)
+ goto data_
+ case SMW_ES, SMW_MS:
+ if (SPEC(sh,mode) == NULL)
+ goto data_
+ }
+ }
+ call sfree (sp)
+ return
+ }
+ }
+ }
+ }
+
+ # Set parameters common to an entire image.
+ if (newim) {
+ call imstats (im, IM_IMAGENAME, IMNAME(sh), LEN_SHDRS)
+ IM(sh) = im
+ MW(sh) = smw
+
+ # Get standard parameters.
+ call shdr_gi (im, "OFLAG", OBJECT, OFLAG(sh))
+ call shdr_gr (im, "EXPOSURE", INDEF, IT(sh))
+ call shdr_gr (im, "ITIME", IT(sh), IT(sh))
+ call shdr_gr (im, "EXPTIME", IT(sh), IT(sh))
+ call shdr_gr (im, "RA", INDEF, RA(sh))
+ call shdr_gr (im, "DEC", INDEF, DEC(sh))
+ call shdr_gr (im, "UT", INDEF, UT(sh))
+ call shdr_gr (im, "ST", INDEF, ST(sh))
+ call shdr_gr (im, "HA", INDEF, HA(sh))
+ call shdr_gr (im, "AIRMASS", INDEF, AM(sh))
+ call shdr_gi (im, "DC-FLAG", DCNO, DC(sh))
+ call shdr_gi (im, "EX-FLAG", ECNO, EC(sh))
+ call shdr_gi (im, "CA-FLAG", FCNO, FC(sh))
+ iferr (call imgstr (im, "DEREDDEN", RC(sh), LEN_SHDRS))
+ RC(sh) = EOS
+
+ # Flag bad airmass value; i.e. 0.
+ if (!IS_INDEF (AM(sh)) && AM(sh) < 1.)
+ AM(sh) = INDEF
+
+ # Set the SMW information.
+ if (SMW_FORMAT(smw) == SMW_MS)
+ i = 3B
+ else
+ i = 2 ** (SMW_PAXIS(smw,1) - 1)
+ CTLW1(sh) = smw_sctran (smw, "logical", "world", i)
+ CTWL1(sh) = smw_sctran (smw, "world", "logical", i)
+
+ # Set the units.
+ mw = SMW_MW(smw,0)
+ i = SMW_PAXIS(smw,1)
+ iferr (call mw_gwattrs (mw, i, "label", LABEL(sh),LEN_SHDRS))
+ call strcpy ("", LABEL(sh), LEN_SHDRS)
+ if (streq (LABEL(sh), "equispec") || streq (LABEL(sh), "multispe"))
+ call strcpy ("", LABEL(sh), LEN_SHDRS)
+ iferr (call mw_gwattrs (mw, i, "units", UNITS(sh),LEN_SHDRS)) {
+ call sprintf (Memc[key], SZ_FNAME, "cunit%d")
+ call pargi (i)
+ iferr (call imgstr (im, Memc[key], UNITS(sh), LEN_SHDRS)) {
+ call strlwr (LABEL(sh))
+ if (LABEL(sh) == EOS)
+ call strcpy ("", UNITS(sh), LEN_SHDRS)
+ else if (streq (LABEL(sh), "lambda"))
+ call strcpy ("angstroms", UNITS(sh), LEN_SHDRS)
+ else if (streq (LABEL(sh), "freq"))
+ call strcpy ("hertz", UNITS(sh), LEN_SHDRS)
+ else if (strncmp (LABEL(sh), "velo", 4) == 0)
+ call strcpy ("m/s", UNITS(sh), LEN_SHDRS)
+ else if (streq (LABEL(sh), "waveleng"))
+ call strcpy ("angstroms", UNITS(sh), LEN_SHDRS)
+ else
+ call strcpy ("", UNITS(sh), LEN_SHDRS)
+ }
+ if (strncmp (LABEL(sh), "velo", 4) == 0)
+ call strcat (" 21 centimeters", UNITS(sh), LEN_SHDRS)
+ }
+ if (UNITS(sh) == EOS && DC(sh) != DCNO)
+ call strcpy ("Angstroms", UNITS(sh), LEN_SHDRS)
+ MWUN(sh) = un_open (UNITS(sh))
+ call un_copy (MWUN(sh), UN(sh))
+
+ iferr (call imgstr (im, "bunit", Memc[str], SZ_LINE))
+ call strcpy ("", Memc[str], SZ_LINE)
+ FUNIM(sh) = fun_open (Memc[str])
+ if (FUN_CLASS(FUNIM(sh)) != FUN_UNKNOWN)
+ FC(sh) = FCYES
+
+ call fun_copy (FUNIM(sh), FUN(sh))
+ call strcpy (FUN_LABEL(FUN(sh)), FLABEL(sh), LEN_SHDRS)
+ call strcpy (FUN_UNITS(FUN(sh)), FUNITS(sh), LEN_SHDRS)
+ }
+
+ # Set WCS parameters for spectrum type.
+ switch (SMW_FORMAT(smw)) {
+ case SMW_ND:
+ # Set physical and logical indices.
+ if (!IS_INDEFI (ap)) {
+ i = max (1, min (SMW_NSPEC(smw), ap))
+ j = 1
+ } else {
+ i = max (1, index1)
+ j = max (1, index2)
+ }
+ call smw_mw (smw, i, j, mw, k, l)
+
+ LINDEX(sh,1) = max (1, min (SMW_LLEN(smw,2), k))
+ LINDEX(sh,2) = max (1, min (SMW_LLEN(smw,3), l))
+ PINDEX(sh,1) = LINDEX(sh,1)
+ PINDEX(sh,2) = LINDEX(sh,2)
+ APINDEX(sh) = LINDEX(sh,1)
+
+ # Set aperture information. Note the use of the logical index.
+ np1 = 1
+ call smw_gwattrs (smw, i, j, AP(sh), BEAM(sh), DC(sh),
+ dval, dval, np2, dval, APLOW(sh,1), APHIGH(sh,1), coeff)
+
+ call smw_gapid (smw, i, j, TITLE(sh), LEN_SHDRS)
+ Memc[SID(sh,1)] = EOS
+
+ switch (SMW_LDIM(smw)) {
+ case 1:
+ IMSEC(sh) = EOS
+ case 2:
+ if (APLOW(sh,1) == APHIGH(sh,1)) {
+ if (SMW_LAXIS(smw,1) == 1)
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d]")
+ else
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,*]")
+ call pargi (nint (APLOW(sh,1)))
+ } else {
+ if (SMW_LAXIS(smw,1) == 1)
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d:%d]")
+ else
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,*]")
+ call pargi (nint (APLOW(sh,1)))
+ call pargi (nint (APHIGH(sh,1)))
+ }
+ case 3:
+ if (APLOW(sh,1)==APHIGH(sh,1) && APLOW(sh,2)==APHIGH(sh,2)) {
+ switch (SMW_LAXIS(smw,1)) {
+ case 1:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d,%d]")
+ case 2:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,*,%d]")
+ case 3:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,%d,*]")
+ }
+ call pargi (nint (APLOW(sh,1)))
+ call pargi (nint (APLOW(sh,2)))
+ } else if (APLOW(sh,1) == APHIGH(sh,1)) {
+ switch (SMW_LAXIS(smw,1)) {
+ case 1:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d,%d:%d]")
+ case 2:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,*,%d:%d]")
+ case 3:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d,%d:%d,*]")
+ }
+ call pargi (nint (APLOW(sh,1)))
+ call pargi (nint (APLOW(sh,2)))
+ call pargi (nint (APHIGH(sh,2)))
+ } else if (APLOW(sh,2) == APHIGH(sh,2)) {
+ switch (SMW_LAXIS(smw,1)) {
+ case 1:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d:%d,%d]")
+ case 2:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,*,%d]")
+ case 3:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,%d,*]")
+ }
+ call pargi (nint (APLOW(sh,1)))
+ call pargi (nint (APHIGH(sh,1)))
+ call pargi (nint (APLOW(sh,2)))
+ } else {
+ switch (SMW_LAXIS(smw,1)) {
+ case 1:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d:%d,%d:%d]")
+ case 2:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,*,%d:%d]")
+ case 3:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[%d:%d,%d:%d,*]")
+ }
+ call pargi (nint (APLOW(sh,1)))
+ call pargi (nint (APHIGH(sh,1)))
+ call pargi (nint (APLOW(sh,2)))
+ call pargi (nint (APHIGH(sh,2)))
+ }
+ }
+
+ case SMW_ES, SMW_MS:
+ # Set the image and aperture indices.
+ if (SMW_PAXIS(smw,2) != 3) {
+ PINDEX(sh,1) = max (1, min (SMW_LLEN(smw,2), index1))
+ PINDEX(sh,2) = max (1, min (SMW_LLEN(smw,3), index2))
+ LINDEX(sh,1) = PINDEX(sh,1)
+ LINDEX(sh,2) = PINDEX(sh,2)
+ APINDEX(sh) = LINDEX(sh,1)
+ } else {
+ PINDEX(sh,1) = 1
+ PINDEX(sh,2) = max (1, min (SMW_LLEN(smw,2), index2))
+ LINDEX(sh,1) = PINDEX(sh,2)
+ LINDEX(sh,2) = 1
+ APINDEX(sh) = 1
+ }
+
+ # If an aperture is specified first try and find it.
+ # If it is not specified or found then use the physical index.
+
+ coeff = NULL
+ AP(sh) = 0
+ if (!IS_INDEFI(ap)) {
+ do i = 1, SMW_NSPEC(smw) {
+ call smw_gwattrs (smw, i, 1, AP(sh), BEAM(sh), DC(sh),
+ dval, dval, np2, dval, APLOW(sh,1), APHIGH(sh,1), coeff)
+ if (AP(sh) == ap && SMW_PAXIS(smw,2) != 3) {
+ PINDEX(sh,1) = i
+ LINDEX(sh,1) = i
+ APINDEX(sh) = i
+ break
+ }
+ }
+ }
+ if (AP(sh) != ap)
+ call smw_gwattrs (smw, APINDEX(sh), 1, AP(sh), BEAM(sh), DC(sh),
+ dval, dval, np2, dval, APLOW(sh,1), APHIGH(sh,1), coeff)
+ call mfree (coeff, TY_CHAR)
+
+ np1 = 1
+ if (SMW_PDIM(smw) > 1) {
+ ct = smw_sctran (smw, "logical", "physical", 2)
+ PINDEX(sh,1) = nint (smw_c1tranr (ct, real(PINDEX(sh,1))))
+ call smw_ctfree (ct)
+ }
+ if (SMW_PDIM(smw) > 2) {
+ ct = smw_sctran (smw, "logical", "physical", 4)
+ PINDEX(sh,2) = nint (smw_c1tranr (ct, real(PINDEX(sh,2))))
+ call smw_ctfree (ct)
+ }
+
+ call smw_gapid (smw, APINDEX(sh), 1, TITLE(sh), LEN_SHDRS)
+ call shdr_type (sh, 1, PINDEX(sh,2))
+
+ switch (SMW_LDIM(smw)) {
+ case 1:
+ IMSEC(sh) = EOS
+ case 2:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d]")
+ call pargi (APINDEX(sh))
+ case 3:
+ call sprintf (IMSEC(sh), LEN_SHDRS, "[*,%d,%d]")
+ call pargi (APINDEX(sh))
+ call pargi (LINDEX(sh,2))
+ }
+ }
+
+ # Set NP1 and NP2 in logical coordinates.
+ i = 2 ** (SMW_PAXIS(smw,1) - 1)
+ ct = smw_sctran (smw, "physical", "logical", i)
+ i = max (1, min (int (smw_c1tranr (ct, real (np1))), SMW_LLEN(smw,1)))
+ j = max (1, min (int (smw_c1tranr (ct, real (np2))), SMW_LLEN(smw,1)))
+ call smw_ctfree (ct)
+ np1 = min (i, j)
+ np2 = max (i, j)
+ np = np2 - np1 + 1
+
+ NP1(sh) = np1
+ NP2(sh) = np2
+ SN(sh) = np
+
+
+data_ # Set the coordinate and data arrays if desired otherwise free them.
+ CTLW(sh) = CTLW1(sh)
+ CTWL(sh) = CTWL1(sh)
+
+ # Set linear dispersion terms.
+ W0(sh) = shdr_lw (sh, double(1))
+ W1(sh) = shdr_lw (sh, double(np))
+ WP(sh) = (W1(sh) - W0(sh)) / (np - 1)
+ SN(sh) = np
+
+ if (mode == SHHDR) {
+ do i = 1, SH_NTYPES
+ call mfree (SPEC(sh,i), TY_REAL)
+ call sfree (sp)
+ return
+ }
+
+ # Set WCS array
+ if (SX(sh) == NULL)
+ call malloc (SX(sh), np, TY_REAL)
+ else
+ call realloc (SX(sh), np, TY_REAL)
+ do i = 1, np
+ Memr[SX(sh)+i-1] = shdr_lw (sh, double(i))
+
+ # Set spectrum array in most efficient way.
+ switch (SMW_FORMAT(smw)) {
+ case SMW_ND:
+ if (mode == SHDATA || SY(sh) == NULL) {
+ if (SY(sh) == NULL)
+ call malloc (SY(sh), np, TY_REAL)
+ else
+ call realloc (SY(sh), np, TY_REAL)
+ call aclrr (Memr[SY(sh)], np)
+ if (IS_INDEF(APLOW(sh,1)))
+ aplow[1] = 1
+ else
+ aplow[1] = nint (APLOW(sh,1))
+ if (IS_INDEF(APHIGH(sh,1)))
+ aphigh[1] = 1
+ else
+ aphigh[1] = nint (APHIGH(sh,1))
+ if (IS_INDEF(APLOW(sh,2)))
+ aplow[2] = 1
+ else
+ aplow[2] = nint (APLOW(sh,2))
+ if (IS_INDEF(APHIGH(sh,2)))
+ aphigh[2] = 1
+ else
+ aphigh[2] = nint (APHIGH(sh,2))
+ k = aplow[1]
+ l = aphigh[1]
+ n = aphigh[1] - aplow[1] + 1
+ if (SMW_LAXIS(smw,1) == 1) {
+ do j = aplow[2], aphigh[2] {
+ do i = aplow[1], aphigh[1] {
+ buf = imgs3r (im, np1, np2, i, i, j, j)
+ call aaddr (Memr[buf], Memr[SY(sh)],
+ Memr[SY(sh)], np)
+ }
+ }
+ } else if (SMW_LAXIS(smw,1) == 2) {
+ do j = aplow[2], aphigh[2] {
+ do i = np1, np2 {
+ buf = imgs3r (im, k, l, i, i, j, j)
+ Memr[SY(sh)+i-np1] = Memr[SY(sh)+i-np1] +
+ asumr (Memr[buf], n)
+ }
+ }
+ } else {
+ do i = np1, np2 {
+ do j = aplow[2], aphigh[2] {
+ buf = imgs3r (im, k, l, j, j, i, i)
+ Memr[SY(sh)+i-np1] = Memr[SY(sh)+i-np1] +
+ asumr (Memr[buf], n)
+ }
+ }
+ }
+ }
+ case SMW_ES, SMW_MS:
+ if (mode == SHDATA || SY(sh) == NULL) {
+ if (SY(sh) == NULL)
+ call malloc (SY(sh), np, TY_REAL)
+ else
+ call realloc (SY(sh), np, TY_REAL)
+ i = LINDEX(sh,1)
+ j = LINDEX(sh,2)
+ buf = imgs3r (im, np1, np2, i, i, j, j)
+ call amovr (Memr[buf], Memr[SY(sh)], np)
+ }
+
+ if (mode > SHDATA)
+ call shdr_gtype (sh, mode)
+ }
+
+ # Guess flux scale if necessary.
+ if (FC(sh) == FCYES && FUN_CLASS(FUNIM(sh)) == FUN_UNKNOWN) {
+ if (Memr[SY(sh)+np/2] < 1e-18)
+ call strcpy ("erg/cm2/s/Hz", Memc[str], SZ_LINE)
+ else if (Memr[SY(sh)+np/2] < 1e-5)
+ call strcpy ("erg/cm2/s/A", Memc[str], SZ_LINE)
+ call fun_close (FUNIM(sh))
+ FUNIM(sh) = fun_open (Memc[str])
+ if (FUN_CLASS(FUN(sh)) == FUN_UNKNOWN) {
+ call fun_copy (FUNIM(sh), FUN(sh))
+ call strcpy (FUN_LABEL(FUN(sh)), FLABEL(sh), LEN_SHDRS)
+ call strcpy (FUN_UNITS(FUN(sh)), FUNITS(sh), LEN_SHDRS)
+ }
+ }
+ if (SPEC(sh,mode) != 0)
+ iferr (call fun_ctranr (FUNIM(sh), FUN(sh), UN(sh), Memr[SX(sh)],
+ Memr[SPEC(sh,mode)], Memr[SPEC(sh,mode)], SN(sh)))
+ ;
+
+ call sfree (sp)
+end
+
+
+# SHDR_GTYPE -- Get the selected spectrum type.
+# Currently this only works for multispec data.
+
+procedure shdr_gtype (sh, type)
+
+pointer sh #I SHDR pointer
+int type #I Spectrum type
+
+int i, j, ctowrd(), strdic()
+pointer sp, key, str, im, smw, ct, buf, smw_sctran(), imgs3r()
+real smw_c1tranr()
+
+begin
+ im = IM(sh)
+ smw = MW(sh)
+
+ if (SMW_FORMAT(smw) == SMW_ND)
+ return
+ if (SMW_PDIM(smw) < 3) {
+ if (type != SHDATA && type != SHRAW) {
+ if (SID(sh,type) != NULL)
+ call mfree (SID(sh,type), TY_CHAR)
+ if (SPEC(sh,type) != NULL)
+ call mfree (SPEC(sh,type), TY_REAL)
+ }
+ return
+ }
+
+ # Find the band.
+ call smark (sp)
+ call salloc (key, SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ do i = 1, 5 {
+ call sprintf (Memc[key], SZ_LINE, "BANDID%d")
+ call pargi (i)
+ ifnoerr (call imgstr (im, Memc[key], Memc[str], SZ_LINE)) {
+ j = 1
+ if (ctowrd (Memc[str], j, Memc[key], SZ_LINE) == 0)
+ next
+ if (strdic (Memc[key], Memc[key], SZ_LINE, STYPES) != type)
+ next
+ if (SID(sh,type) == NULL)
+ call malloc (SID(sh,type), LEN_SHDRS, TY_CHAR)
+ call strcpy (Memc[str], Memc[SID(sh,type)], LEN_SHDRS)
+ STYPE(sh,type) = type
+ break
+ }
+ }
+ call sfree (sp)
+ if (i == 6) {
+ if (SID(sh,type) != NULL)
+ call mfree (SID(sh,type), TY_CHAR)
+ if (SPEC(sh,type) != NULL)
+ call mfree (SPEC(sh,type), TY_REAL)
+ return
+ }
+
+ # Map the physical band to logical vector.
+ ct = smw_sctran (smw, "physical", "logical", 4)
+ i = nint (smw_c1tranr (ct, real(i)))
+ call smw_ctfree (ct)
+ if (SMW_PAXIS(smw,2) != 3) {
+ if (i > SMW_LLEN(smw,3))
+ return
+ j = i
+ i = LINDEX(sh,1)
+ } else {
+ if (i > SMW_LLEN(smw,2))
+ return
+ j = 1
+ }
+
+ # Get the spectrum.
+ if (SPEC(sh,type) == NULL)
+ call malloc (SPEC(sh,type), SN(sh), TY_REAL)
+ else
+ call realloc (SPEC(sh,type), SN(sh), TY_REAL)
+ buf = imgs3r (im, NP1(sh), NP2(sh), i, i, j, j)
+ call amovr (Memr[buf], Memr[SPEC(sh,type)], SN(sh))
+end
+
+
+# SHDR_TYPE -- Determine the spectrum type.
+# Currently this only works for multispec data.
+
+procedure shdr_type (sh, index, band)
+
+pointer sh #I SHDR pointer
+int index #I Index
+int band #I Physical band
+
+int i, ctowrd(), strdic()
+pointer sp, key
+
+begin
+ if (SMW_FORMAT(MW(sh)) == SMW_ND)
+ return
+
+ call smark (sp)
+ call salloc (key, SZ_LINE, TY_CHAR)
+
+ if (SID(sh,index) == NULL)
+ call malloc (SID(sh,index), LEN_SHDRS, TY_CHAR)
+
+ call sprintf (Memc[key], SZ_FNAME, "BANDID%d")
+ call pargi (band)
+ iferr (call imgstr (IM(sh), Memc[key], Memc[SID(sh,index)], LEN_SHDRS))
+ Memc[SID(sh,index)] = EOS
+
+ i = 1
+ if (ctowrd (Memc[SID(sh,index)], i, Memc[key], SZ_LINE) > 0)
+ STYPE(sh,index) = strdic (Memc[key], Memc[key], SZ_LINE, STYPES)
+ else
+ STYPE(sh,index) = 0
+
+ call sfree (sp)
+end
+
+
+# SHDR_CLOSE -- Close and free the SHDR structure.
+
+procedure shdr_close (sh)
+
+pointer sh # SHDR structure
+int i
+
+begin
+ if (sh == NULL)
+ return
+ do i = 1, SH_NTYPES {
+ call mfree (SPEC(sh,i), TY_REAL)
+ call mfree (SID(sh,i), TY_CHAR)
+ }
+ call un_close (UN(sh))
+ call un_close (MWUN(sh))
+ call fun_close (FUN(sh))
+ call fun_close (FUNIM(sh))
+ if (MW(sh) != NULL) {
+ call smw_ctfree (CTLW1(sh))
+ call smw_ctfree (CTWL1(sh))
+ }
+ call mfree (sh, TY_STRUCT)
+end
+
+
+# SHDR_COPY -- Make a copy of an SHDR structure.
+# The image pointer is not copied and the MWCS pointer and transform pointers
+# may or may not be copied . The uncopied pointers mean that they will be
+# shared by multiple spectrum structures but it also means that when they are
+# closed the structures will have invalid pointers. The advantage of not
+# copying is that many spectra may come from the same image and the overhead
+# of having copies of the IMIO and MWCS pointers can be avoided.
+
+procedure shdr_copy (sh1, sh2, wcs)
+
+pointer sh1 # SHDR structure to copy
+pointer sh2 # SHDR structure copy
+int wcs # Make copy of wcs?
+
+int i
+pointer un, mwun, fun, funim, spec[SH_NTYPES], sid[SH_NTYPES], smw_newcopy()
+errchk shdr_system
+
+begin
+ if (sh2 == NULL) {
+ call calloc (sh2, LEN_SHDR, TY_STRUCT)
+ call calloc (SID(sh2,1), LEN_SHDRS, TY_CHAR)
+ }
+
+ un = UN(sh2)
+ mwun = MWUN(sh2)
+ fun = FUN(sh2)
+ funim = FUNIM(sh2)
+ call amovi (SPEC(sh2,1), spec, SH_NTYPES)
+ call amovi (SID(sh2,1), sid, SH_NTYPES)
+ call amovi (Memi[sh1], Memi[sh2], LEN_SHDR)
+ call amovi (spec, SPEC(sh2,1), SH_NTYPES)
+ call amovi (sid, SID(sh2,1), SH_NTYPES)
+ UN(sh2) = un
+ MWUN(sh2) = mwun
+ FUN(sh2) = fun
+ FUNIM(sh2) = funim
+ call un_copy (UN(sh1), UN(sh2))
+ call un_copy (MWUN(sh1), MWUN(sh2))
+ call fun_copy (FUN(sh1), FUN(sh2))
+ call fun_copy (FUNIM(sh1), FUNIM(sh2))
+ do i = 1, SH_NTYPES {
+ if (SPEC(sh1,i) != NULL) {
+ if (SPEC(sh2,i) == NULL)
+ call malloc (SPEC(sh2,i), SN(sh1), TY_REAL)
+ else
+ call realloc (SPEC(sh2,i), SN(sh1), TY_REAL)
+ call amovr (Memr[SPEC(sh1,i)], Memr[SPEC(sh2,i)], SN(sh1))
+ }
+ }
+
+ if (wcs == YES && MW(sh1) != NULL) {
+ MW(sh2) = smw_newcopy (MW(sh1))
+ CTLW1(sh2) = NULL
+ CTWL1(sh2) = NULL
+ call shdr_system (sh2, "world")
+ }
+end
+
+
+# SHDR_SYSTEM -- Set or change the WCS system.
+
+procedure shdr_system (sh, system)
+
+pointer sh # SHDR pointer
+char system[ARB] # System
+
+int i, sn
+bool streq()
+double shdr_lw()
+pointer smw, mw, smw_sctran(), un_open()
+errchk smw_sctran, un_open
+
+begin
+ smw = MW(sh)
+ if (smw == NULL)
+ call error (1, "shdr_system: MWCS not defined")
+
+ call smw_ctfree (CTLW1(sh))
+ call smw_ctfree (CTWL1(sh))
+
+ switch (SMW_FORMAT(smw)) {
+ case SMW_ND, SMW_ES:
+ i = 2 ** (SMW_PAXIS(smw,1) - 1)
+ CTLW1(sh) = smw_sctran (smw, "logical", system, i)
+ CTWL1(sh) = smw_sctran (smw, system, "logical", i)
+ case SMW_MS:
+ CTLW1(sh) = smw_sctran (smw, "logical", system, 3B)
+ CTWL1(sh) = smw_sctran (smw, system, "logical", 3B)
+ }
+ CTLW(sh) = CTLW1(sh)
+ CTWL(sh) = CTWL1(sh)
+
+ # Set labels and units
+ call un_close (MWUN(sh))
+ if (streq (system, "physical")) {
+ call strcpy ("Pixel", LABEL(sh), LEN_SHDRS)
+ call strcpy ("", UNITS(sh), LEN_SHDRS)
+ MWUN(sh) = un_open (UNITS(sh))
+ } else {
+ call smw_mw (smw, 1, 1, mw, i, i)
+ iferr (call mw_gwattrs (mw, SMW_PAXIS(smw,1), "label", LABEL(sh),
+ LEN_SHDRS))
+ call strcpy ("", LABEL(sh), LEN_SHDRS)
+ if (streq (LABEL(sh), "equispec") || streq (LABEL(sh), "multispe"))
+ call strcpy ("", LABEL(sh), LEN_SHDRS)
+ iferr (call mw_gwattrs (mw, SMW_PAXIS(smw,1), "units", UNITS(sh),
+ LEN_SHDRS))
+ call strcpy ("", UNITS(sh), LEN_SHDRS)
+ MWUN(sh) = un_open (UNITS(sh))
+ call strcpy (UN_LABEL(UN(sh)), LABEL(sh), LEN_SHDRS)
+ call strcpy (UN_UNITS(UN(sh)), UNITS(sh), LEN_SHDRS)
+ }
+
+ sn = SN(sh)
+ W0(sh) = shdr_lw (sh, double(1))
+ W1(sh) = shdr_lw (sh, double(sn))
+ WP(sh) = (W1(sh) - W0(sh)) / (sn - 1)
+ if (SX(sh) != NULL)
+ do i = 1, sn
+ Memr[SX(sh)+i-1] = shdr_lw (sh, double(i))
+end
+
+
+# SHDR_UNITS -- Set or change the WCS system.
+# This changes W0, W1, WP, and SX.
+
+procedure shdr_units (sh, units)
+
+pointer sh # SHDR pointer
+char units[ARB] # Units
+
+int i, sn
+bool streq()
+double shdr_lw()
+pointer str, un, un_open()
+errchk un_open
+
+begin
+ # Check for unknown units.
+ if (streq (units, "display")) {
+ call malloc (str, SZ_LINE, TY_CHAR)
+ iferr (call mw_gwattrs (SMW_MW(MW(sh),0), SMW_PAXIS(MW(sh),1),
+ "units_display", Memc[str], SZ_LINE)) {
+ un = NULL
+ call un_copy (MWUN(sh), un)
+ } else
+ un = un_open (Memc[str])
+ call mfree (str, TY_CHAR)
+ } else if (streq (units, "default")) {
+ un = NULL
+ call un_copy (MWUN(sh), un)
+ } else
+ un = un_open (units)
+ if (UN_CLASS(un) == UN_UNKNOWN || UN_CLASS(MWUN(sh)) == UN_UNKNOWN) {
+ call un_close (un)
+ call error (1, "Cannot convert to specified units")
+ }
+
+ # Update the coordinates.
+ call un_close (UN(sh))
+ UN(sh) = un
+
+ call strcpy (UN_LABEL(UN(sh)), LABEL(sh), LEN_SHDRS)
+ call strcpy (UN_UNITS(UN(sh)), UNITS(sh), LEN_SHDRS)
+
+ sn = SN(sh)
+ W0(sh) = shdr_lw (sh, double(1))
+ W1(sh) = shdr_lw (sh, double(sn))
+ WP(sh) = (W1(sh) - W0(sh)) / (sn - 1)
+ if (SX(sh) != NULL)
+ do i = 1, sn
+ Memr[SX(sh)+i-1] = shdr_lw (sh, double(i))
+end
+
+
+# SHDR_LW -- Logical to world coordinate transformation.
+# The transformation pointer is generally NULL only after SHDR_LINEAR
+
+double procedure shdr_lw (sh, l)
+
+pointer sh # SHDR pointer
+double l # Logical coordinate
+double w # World coordinate
+
+double l0, l1, l2, w1, smw_c1trand()
+
+begin
+ l0 = l + NP1(sh) - 1
+ if (CTLW(sh) != NULL) {
+ switch (SMW_FORMAT(MW(sh))) {
+ case SMW_ND, SMW_ES:
+ w = smw_c1trand (CTLW(sh), l0)
+ case SMW_MS:
+ call smw_c2trand (CTLW(sh), l0, double (APINDEX(sh)), w, w1)
+ }
+ } else {
+ switch (DC(sh)) {
+ case DCLOG:
+ w = W0(sh) * 10. ** (log10(W1(sh)/W0(sh)) * (l0-1) / (SN(sh)-1))
+ case DCFUNC:
+ w = W0(sh)
+ call smw_c2trand (CTWL1(sh), w, double (AP(sh)), l1, w1)
+ w = W1(sh)
+ call smw_c2trand (CTWL1(sh), w, double (AP(sh)), l2, w1)
+ if (SN(sh) > 1)
+ l1 = (l2 - l1) / (SN(sh) - 1) * (l0 - 1) + l1
+ else
+ l1 = l0 - 1 + l1
+ call smw_c2trand (CTLW1(sh), l1, double (APINDEX(sh)), w, w1)
+ default:
+ w = W0(sh) + (l0 - 1) * WP(sh)
+ }
+ }
+
+ iferr (call un_ctrand (MWUN(sh), UN(sh), w, w, 1))
+ ;
+ return (w)
+end
+
+
+# SHDR_WL -- World to logical coordinate transformation.
+# The transformation pointer is generally NULL only after SHDR_LINEAR
+
+double procedure shdr_wl (sh, w)
+
+pointer sh # SHDR pointer
+double w # World coordinate
+double l # Logical coordinate
+
+double w1, l1, l2, smw_c1trand()
+
+begin
+ iferr (call un_ctrand (UN(sh), MWUN(sh), w, w1, 1))
+ w1 = w
+
+ if (CTWL(sh) != NULL) {
+ switch (SMW_FORMAT(MW(sh))) {
+ case SMW_ND, SMW_ES:
+ l = smw_c1trand (CTWL(sh), w1)
+ case SMW_MS:
+ call smw_c2trand (CTWL(sh), w1, double (AP(sh)),l,l1)
+ }
+ } else {
+ switch (DC(sh)) {
+ case DCLOG:
+ l = log10(w1/W0(sh)) / log10(W1(sh)/W0(sh)) * (SN(sh)-1) + 1
+ case DCFUNC:
+ call smw_c2trand (CTWL1(sh), w1, double (AP(sh)), l, l1)
+
+ w1 = W0(sh)
+ call smw_c2trand (CTWL1(sh), w1, double (AP(sh)), l1, w1)
+ w1 = W1(sh)
+ call smw_c2trand (CTWL1(sh), w1, double (AP(sh)), l2, w1)
+ if (l1 != l2)
+ l = (SN(sh) - 1) / (l2 - l1) * (l - l1) + 1
+ else
+ l = l - l1 + 1
+ default:
+ l = (w1 - W0(sh)) / WP(sh) + 1
+ }
+ }
+
+ return (l-NP1(sh)+1)
+end
+
+
+# SHDR_REBIN -- Rebin spectrum to dispersion of reference spectrum.
+# The interpolation function is set by ONEDINTERP.
+
+procedure shdr_rebin (sh, shref)
+
+pointer sh # Spectrum to be rebinned
+pointer shref # Reference spectrum
+
+char interp[10]
+int i, j, type, ia, ib, n, clgwrd()
+real a, b, sum, asieval(), asigrl()
+double x, w, xmin, xmax, shdr_lw(), shdr_wl()
+pointer unref, unsave, asi, spec
+bool fp_equalr()
+
+begin
+ # Check if rebinning is needed
+ if (DC(sh) == DC(shref) && DC(sh) != DCFUNC &&
+ fp_equalr (W0(sh), W0(shref)) && fp_equalr(WP(sh), WP(shref)) &&
+ SN(sh) == SN(shref))
+ return
+
+ # Do everything in units of reference WCS.
+ unref = UN(shref)
+ unsave = UN(sh)
+ UN(SH) = unref
+
+ call asiinit (asi, clgwrd ("interp", interp, 10, II_FUNCTIONS))
+ do type = 1, SH_NTYPES {
+ if (SPEC(sh,type) == NULL)
+ next
+
+ # Fit the interpolation function to the spectrum.
+ # Extend the interpolation by one pixel at each end.
+
+ n = SN(sh)
+ call malloc (spec, n+2, TY_REAL)
+ call amovr (Memr[SPEC(sh,type)], Memr[spec+1], n)
+ Memr[spec] = Memr[SPEC(sh,type)]
+ Memr[spec+n+1] = Memr[SPEC(sh,type)+n-1]
+ call asifit (asi, Memr[spec], n+2)
+ call mfree (spec, TY_REAL)
+
+ xmin = 0.5
+ xmax = n + 0.5
+
+ # Reallocate data array
+ if (n != SN(shref)) {
+ n = SN(shref)
+ call realloc (SPEC(sh,type), n, TY_REAL)
+ call aclrr (Memr[SPEC(sh,type)], n)
+ }
+ spec = SPEC(sh,type)
+
+ # Compute the average flux in each output pixel.
+
+ x = 0.5
+ w = shdr_lw (shref, x)
+ x = shdr_wl (sh, w)
+ b = max (xmin, min (xmax, x)) + 1
+ do i = 1, n {
+ x = i + 0.5
+ w = shdr_lw (shref, x)
+ x = shdr_wl (sh, w)
+ a = b
+ b = max (xmin, min (xmax, x)) + 1
+ if (a <= b) {
+ ia = nint (a + 0.5)
+ ib = nint (b - 0.5)
+ if (abs (a+0.5-ia) < .00001 && abs (b-0.5-ib) < .00001) {
+ sum = 0.
+ do j = ia, ib
+ sum = sum + asieval (asi, real(j))
+ if (ib - ia > 0)
+ sum = sum / (ib - ia)
+ } else {
+ sum = asigrl (asi, a, b)
+ if (b - a > 0.)
+ sum = sum / (b - a)
+ }
+ } else {
+ ib = nint (b + 0.5)
+ ia = nint (a - 0.5)
+ if (abs (a-0.5-ia) < .00001 && abs (b+0.5-ib) < .00001) {
+ sum = 0.
+ do j = ib, ia
+ sum = sum + asieval (asi, real(j))
+ if (ia - ib > 0)
+ sum = sum / (ia - ib)
+ } else {
+ sum = asigrl (asi, b, a)
+ if (a - b > 0.)
+ sum = sum / (a - b)
+ }
+ }
+
+ Memr[spec] = sum
+ spec = spec + 1
+ }
+ }
+ call asifree (asi)
+
+ # Set the rest of the header. The coordinate transformations are
+ # canceled to indicate they are not valid for the data. They
+ # are not freed because the same pointer may be used in other
+ # spectra from the same image.
+
+ if (SN(sh) != n)
+ call realloc (SX(sh), n, TY_REAL)
+ call amovr (Memr[SX(shref)], Memr[SX(sh)], n)
+ DC(sh) = DC(shref)
+ W0(sh) = W0(shref)
+ W1(sh) = W1(shref)
+ WP(sh) = WP(shref)
+ SN(sh) = SN(shref)
+
+ CTLW(sh) = NULL
+ CTWL(sh) = NULL
+
+ # Restore original units
+ UN(sh) = unsave
+ iferr (call un_ctranr (unref, UN(sh), Memr[SX(sh)], Memr[SX(sh)],
+ SN(sh)))
+ ;
+end
+
+
+# SHDR_LINEAR -- Rebin spectrum to linear dispersion.
+# The interpolation function is set by ONEDINTERP
+
+procedure shdr_linear (sh, w0, w1, sn, dc)
+
+pointer sh # Spectrum to be rebinned
+real w0 # Wavelength of first logical pixel
+real w1 # Wavelength of last logical pixel
+int sn # Number of pixels
+int dc # Dispersion type (DCLINEAR | DCLOG)
+
+char interp[10]
+int i, j, type, ia, ib, n, clgwrd()
+real w0mw, w1mw, a, b, sum, asieval(), asigrl()
+double x, w, w0l, wp, xmin, xmax, shdr_wl()
+pointer unsave, asi, spec
+bool fp_equalr()
+
+begin
+ # Check if rebinning is needed
+ if (DC(sh) == dc && fp_equalr (W0(sh), w0) &&
+ fp_equalr (W1(sh), w1) && SN(sh) == sn)
+ return
+
+ # Do everything in units of MWCS.
+ call un_ctranr (UN(sh), MWUN(sh), w0, w0mw, 1)
+ call un_ctranr (UN(sh), MWUN(sh), w1, w1mw, 1)
+ unsave = UN(sh)
+ UN(SH) = MWUN(sh)
+
+ call asiinit (asi, clgwrd ("interp", interp, 10, II_FUNCTIONS))
+ do type = 1, SH_NTYPES {
+ if (SPEC(sh,type) == NULL)
+ next
+
+ # Fit the interpolation function to the spectrum.
+ # Extend the interpolation by one pixel at each end.
+
+ n = SN(sh)
+ call malloc (spec, n+2, TY_REAL)
+ call amovr (Memr[SPEC(sh,type)], Memr[spec+1], n)
+ Memr[spec] = Memr[SPEC(sh,type)]
+ Memr[spec+n+1] = Memr[SPEC(sh,type)+n-1]
+ call asifit (asi, Memr[spec], n+2)
+ call mfree (spec, TY_REAL)
+
+ xmin = 0.5
+ xmax = n + 0.5
+
+ # Reallocate spectrum data array
+ if (n != sn) {
+ n = sn
+ call realloc (SPEC(sh,type), n, TY_REAL)
+ }
+ spec = SPEC(sh,type)
+
+ # Integrate across pixels using ASIGRL.
+
+ x = 0.5
+ if (dc == DCLOG) {
+ w0l = log10 (w0mw)
+ wp = (log10 (w1mw) - log10(w0mw)) / (n - 1)
+ w = 10. ** (w0l+(x-1)*wp)
+ } else {
+ wp = (w1mw - w0mw) / (n - 1)
+ w = w0mw + (x - 1) * wp
+ }
+ x = shdr_wl (sh, w)
+ b = max (xmin, min (xmax, x)) + 1
+ do i = 1, n {
+ x = i + 0.5
+ if (dc == DCLOG)
+ w = 10. ** (w0l + (x - 1) * wp)
+ else
+ w = w0mw + (x - 1) * wp
+ x = shdr_wl (sh, w)
+ a = b
+ b = max (xmin, min (xmax, x)) + 1
+ if (a <= b) {
+ ia = nint (a + 0.5)
+ ib = nint (b - 0.5)
+ if (abs (a+0.5-ia) < .00001 && abs (b-0.5-ib) < .00001) {
+ sum = 0.
+ do j = ia, ib
+ sum = sum + asieval (asi, real(j))
+ if (ib - ia > 0)
+ sum = sum / (ib - ia)
+ } else {
+ sum = asigrl (asi, a, b)
+ if (b - a > 0.)
+ sum = sum / (b - a)
+ }
+ } else {
+ ib = nint (b + 0.5)
+ ia = nint (a - 0.5)
+ if (abs (a-0.5-ia) < .00001 && abs (b+0.5-ib) < .00001) {
+ sum = 0.
+ do j = ib, ia
+ sum = sum + asieval (asi, real(j))
+ if (ia - ib > 0)
+ sum = sum / (ia - ib)
+ } else {
+ sum = asigrl (asi, b, a)
+ if (a - b > 0.)
+ sum = sum / (a - b)
+ }
+ }
+ Memr[spec] = sum
+ spec = spec + 1
+ }
+ }
+ call asifree (asi)
+
+ # Set the rest of the header. The coordinate transformations are
+ # canceled to indicate they are not valid for the data. They
+ # are not freed because the same pointer may be used in other
+ # spectra from the same image.
+
+ if (SN(sh) != n)
+ call realloc (SX(sh), n, TY_REAL)
+ do i = 0, n-1 {
+ if (dc == DCLOG)
+ w = 10. ** (w0l + i * wp)
+ else
+ w = w0mw + i * wp
+ Memr[SX(sh)+i] = w
+ }
+ W0(sh) = w0
+ W1(sh) = w1
+ WP(sh) = (w1 - w0) / (sn - 1)
+ SN(sh) = sn
+ NP1(sh) = 1
+ NP2(sh) = sn
+ DC(sh) = dc
+
+ CTLW(sh) = NULL
+ CTWL(sh) = NULL
+
+ # Restore original units
+ UN(sh) = unsave
+ iferr (call un_ctranr (MWUN(sh), UN(sh), Memr[SX(sh)], Memr[SX(sh)],
+ sn))
+ ;
+end
+
+
+# SHDR_EXTRACT -- Extract a specific wavelength region.
+
+procedure shdr_extract (sh, w1, w2, rebin)
+
+pointer sh # SHDR structure
+real w1 # Starting wavelength
+real w2 # Ending wavelength
+bool rebin # Rebin wavelength region?
+
+int i, j, i1, i2, n
+double l1, l2
+pointer buf
+bool fp_equald()
+double shdr_wl(), shdr_lw()
+errchk shdr_linear, shdr_lw, shdr_wl
+
+begin
+ l1 = shdr_wl (sh, double (w1))
+ l2 = shdr_wl (sh, double (w2))
+ if (fp_equald(l1,l2) || max(l1,l2) < 1 || min (l1,l2) > SN(sh))
+ call error (1, "No pixels to extract")
+ l1 = max (1D0, min (double (SN(sh)), l1))
+ l2 = max (1D0, min (double (SN(sh)), l2))
+ i1 = nint (l1)
+ i2 = nint (l2)
+ n = abs (i2 - i1) + 1
+
+ if (rebin) {
+ l1 = shdr_lw (sh, l1)
+ l2 = shdr_lw (sh, l2)
+ if (DC(sh) == DCFUNC)
+ call shdr_linear (sh, real (l1), real (l2), n, DCLINEAR)
+ else
+ call shdr_linear (sh, real (l1), real (l2), n, DC(sh))
+ } else {
+ if (i1 == 1 && i2 == SN(sh))
+ return
+
+ if (i1 <= i2) {
+ do j = 1, SH_NTYPES
+ if (SPEC(sh,j) != NULL)
+ call amovr (Memr[SPEC(sh,j)+i1-1], Memr[SPEC(sh,j)], n)
+ } else {
+ call malloc (buf, n, TY_REAL)
+ do j = 1, SH_NTYPES {
+ if (SPEC(sh,j) != NULL) {
+ do i = i1, i2, -1
+ Memr[buf+i1-i] = Memr[SPEC(sh,j)+i-1]
+ call amovr (Memr[buf], Memr[SPEC(sh,j)], n)
+ }
+ }
+ call mfree (buf, TY_REAL)
+ }
+ W0(sh) = Memr[SX(sh)]
+ W1(sh) = Memr[SX(sh)+n-1]
+ SN(sh) = n
+ NP1(sh) = 1
+ NP2(sh) = n
+ if (n > 1)
+ WP(sh) = (W1(sh) - W0(sh)) / (SN(sh) - 1)
+ CTLW(sh) = NULL
+ CTWL(sh) = NULL
+ }
+end
+
+
+# SHDR_GI -- Load an integer value from the header.
+
+procedure shdr_gi (im, field, default, ival)
+
+pointer im
+char field[ARB]
+int default
+int ival
+
+int dummy, imaccf(), imgeti()
+
+begin
+ ival = default
+ if (imaccf (im, field) == YES) {
+ iferr (dummy = imgeti (im, field))
+ call erract (EA_WARN)
+ else
+ ival = dummy
+ }
+end
+
+
+# SHDR_GR -- Load a real value from the header.
+
+procedure shdr_gr (im, field, default, rval)
+
+pointer im
+char field[ARB]
+real default
+real rval
+
+int imaccf()
+real dummy, imgetr()
+
+begin
+ rval = default
+ if (imaccf (im, field) == YES) {
+ iferr (dummy = imgetr (im, field))
+ call erract (EA_WARN)
+ else
+ rval = dummy
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 18:52:55 -0400