diff options
Diffstat (limited to 'noao/onedspec/smw')
31 files changed, 4809 insertions, 0 deletions
diff --git a/noao/onedspec/smw/README b/noao/onedspec/smw/README new file mode 100644 index 00000000..2f2c27b0 --- /dev/null +++ b/noao/onedspec/smw/README @@ -0,0 +1,6 @@ +This directory contains interface routines for the spectral world +coordinate systems. The interface has two functions. The first is to +convert various input formats, including old formats, to one of the WCS +formats used by the ONEDSPEC package. These are MULTISPEC, EQUISPEC, and +NDSPEC. The second is to split large numbers of spectra which exceed the +limits of the MWCS for a single WCS into a number of smaller WCS. diff --git a/noao/onedspec/smw/funits.x b/noao/onedspec/smw/funits.x new file mode 100644 index 00000000..e3c8ddf5 --- /dev/null +++ b/noao/onedspec/smw/funits.x @@ -0,0 +1,445 @@ +include <ctype.h> +include <error.h> +include <funits.h> + + +# FUN_OPEN -- Open funits package +# It is allowed to open an unknown funit type + +pointer procedure fun_open (funits) + +char funits[ARB] # Units string +pointer fun # Units pointer returned + +begin + call calloc (fun, FUN_LEN, TY_STRUCT) + iferr (call fun_decode (fun, funits)) { + call fun_close (fun) + call erract (EA_ERROR) + } + return (fun) +end + + +# FUN_CLOSE -- Close funits package + +procedure fun_close (fun) + +pointer fun # Units pointer + +begin + call mfree (fun, TY_STRUCT) +end + + +# FUN_COPY -- Copy funits pointer + +procedure fun_copy (fun1, fun2) + +pointer fun1, fun2 # Units pointers + +begin + if (fun2 == NULL) + call malloc (fun2, FUN_LEN, TY_STRUCT) + call amovi (Memi[fun1], Memi[fun2], FUN_LEN) +end + + +# FUN_DECODE -- Decode funits string and set up funits structure. +# The main work is done in FUN_DECODE1 so that the funits string may +# be recursive; i.e. the funits string may contain other funits strings. + +procedure fun_decode (fun, funits) + +pointer fun # Units pointer +char funits[ARB] # Units string + +bool streq() +pointer sp, funits1, temp +errchk fun_decode1, fun_ctranr + +begin + if (streq (funits, FUN_USER(fun))) + return + + call smark (sp) + call salloc (funits1, SZ_LINE, TY_CHAR) + call salloc (temp, FUN_LEN, TY_STRUCT) + + # Save a copy to restore in case of an error. + call fun_copy (fun, temp) + + iferr (call fun_decode1 (fun, funits, Memc[funits1], SZ_LINE)) { + call fun_copy (temp, fun) + call sfree (sp) + call erract (EA_ERROR) + } + + call sfree (sp) +end + + +# FUN_DECODE1 -- Decode funits string and set up funits structure. +# Unknown funit strings are allowed. + +procedure fun_decode1 (fun, funits, funits1, sz_funits1) + +pointer fun # Units pointer +char funits[ARB] # Units string +char funits1[sz_funits1] # Secondary funits string to return +int sz_funits1 # Size of secondary funits string + +int funmod, funtype +int i, j, k, nscan(), strdic(), strlen() +real funscale +pointer sp, str + +int class[FUN_NUNITS] +real scale[FUN_NUNITS] +data class /FUN_FREQ,FUN_FREQ,FUN_FREQ,FUN_WAVE/ +data scale /FUN_J,FUN_FU,FUN_CGSH,FUN_CGSA/ + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + call strcpy (funits, Memc[str], SZ_FNAME) + call strlwr (Memc[str]) + call sscan (Memc[str]) + funtype = 0 + funmod = 0 + do i = 1, 2 { + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() != i) + break + + j = strdic (Memc[str], Memc[str], SZ_FNAME, FUN_DIC) + for (k=strlen(Memc[str]); k>0 && + (IS_WHITE(Memc[str+k-1]) || Memc[str+k-1]=='\n'); k=k-1) + Memc[str+k-1] = EOS + + if (j > FUN_NUNITS) { + if (funmod != 0) + break + funmod = j - FUN_NUNITS + } else { + funtype = j + break + } + } + i = nscan() + call gargr (funscale) + if (nscan() != i+1) + funscale = 1 + + if (funtype == 0) { + FUN_TYPE(fun) = 0 + FUN_CLASS(fun) = FUN_UNKNOWN + FUN_LABEL(fun) = EOS + call strcpy (funits, FUN_UNITS(fun), SZ_UNITS) + } else { + FUN_TYPE(fun) = funtype + FUN_CLASS(fun) = class[funtype] + FUN_MOD(fun) = funmod + FUN_SCALE(fun) = scale[funtype] * funscale + FUN_LABEL(fun) = EOS + FUN_UNITS(fun) = EOS + call strcpy (funits, FUN_USER(fun), SZ_UNITS) + switch (funmod) { + case FUN_LOG: + call strcat ("Log ", FUN_LABEL(fun), SZ_UNITS) + case FUN_MAG: + call strcat ("Mag ", FUN_LABEL(fun), SZ_UNITS) + } + call strcat ("Flux", FUN_LABEL(fun), SZ_UNITS) + if (funscale != 1) { + call sprintf (FUN_UNITS(fun), SZ_UNITS, "%sx%.1g") + call pargstr (Memc[str]) + call pargr (funscale) + } else { + call sprintf (FUN_UNITS(fun), SZ_UNITS, "%s") + call pargstr (Memc[str]) + } + } + + call sfree (sp) +end + + +# FUN_COMPARE -- Compare two funits + +bool procedure fun_compare (fun1, fun2) + +pointer fun1, fun2 # Units pointers to compare +bool strne() + +begin + if (strne (FUN_UNITS(fun1), FUN_UNITS(fun2))) + return (false) + if (strne (FUN_LABEL(fun1), FUN_LABEL(fun2))) + return (false) + return (true) +end + + +# FUN_CTRANR -- Transform funits +# Error is returned if the transform cannot be made + +procedure fun_ctranr (fun1, fun2, dun, dval, fval1, fval2, nvals) + +pointer fun1 # Input funits pointer +pointer fun2 # Output funits pointer +pointer dun # Input units pointer +real dval[nvals] # Input dispersion values +real fval1[nvals] # Input flux values +real fval2[nvals] # Output flux values +int nvals # Number of values + +int i +real s, lambda +pointer ang, un_open() +bool fun_compare() +errchk un_open, un_ctranr + +begin + if (fun_compare (fun1, fun2)) { + call amovr (fval1, fval2, nvals) + return + } + + if (FUN_CLASS(fun1) == FUN_UNKNOWN || FUN_CLASS(fun2) == FUN_UNKNOWN) + call error (1, "Cannot convert between selected funits") + + call amovr (fval1, fval2, nvals) + + s = FUN_SCALE(fun1) + switch (FUN_MOD(fun1)) { + case FUN_LOG: + do i = 1, nvals + fval2[i] = 10. ** fval2[i] + case FUN_MAG: + do i = 1, nvals + fval2[i] = 10. ** (-0.4 * fval2[i]) + } + switch (FUN_CLASS(fun1)) { + case FUN_FREQ: + do i = 1, nvals + fval2[i] = fval2[i] / s + case FUN_WAVE: + if (FUN_CLASS(fun2) != FUN_WAVE) { + s = s * FUN_VLIGHT + ang = un_open ("angstroms") + do i = 1, nvals { + call un_ctranr (dun, ang, dval[i], lambda, 1) + fval2[i] = fval2[i] / s * lambda**2 + } + call un_close (ang) + } else { + do i = 1, nvals + fval2[i] = fval2[i] / s + } + } + + s = FUN_SCALE(fun2) + switch (FUN_CLASS(fun2)) { + case FUN_FREQ: + do i = 1, nvals + fval2[i] = fval2[i] * s + case FUN_WAVE: + if (FUN_CLASS(fun1) != FUN_WAVE) { + s = s * FUN_VLIGHT + ang = un_open ("angstroms") + do i = 1, nvals { + call un_ctranr (dun, ang, dval[i], lambda, 1) + fval2[i] = fval2[i] * s / lambda**2 + } + call un_close (ang) + } else { + do i = 1, nvals + fval2[i] = fval2[i] * s + } + } + switch (FUN_MOD(fun2)) { + case FUN_LOG: + do i = 1, nvals + fval2[i] = log10 (fval2[i]) + case FUN_MAG: + do i = 1, nvals + fval2[i] = -2.5 * log10 (fval2[i]) + } +end + + +# FUN_CHANGER -- Change funits +# Error is returned if the conversion cannot be made + +procedure fun_changer (fun, funits, dun, dvals, fvals, nvals, update) + +pointer fun # Units pointer (may be changed) +char funits[ARB] # Desired funits +pointer dun # Dispersion units pointer +real dvals[nvals] # Dispersion values +real fvals[nvals] # Flux Values +int nvals # Number of values +int update # Update funits pointer? + +bool streq(), fun_compare() +pointer fun1, fun_open() +errchk fun_open, fun_ctranr + +begin + + # Check for same funit string + if (streq (funits, FUN_USER(fun))) + return + + # Check for error in funits string, or the same funits. + fun1 = fun_open (funits) + if (fun_compare (fun1, fun)) { + call strcpy (funits, FUN_USER(fun), SZ_UNITS) + call fun_close (fun1) + return + } + + iferr { + call fun_ctranr (fun, fun1, dun, dvals, fvals, fvals, nvals) + if (update == YES) + call fun_copy (fun1, fun) + call fun_close(fun1) + } then { + call fun_close(fun1) + call erract (EA_ERROR) + } +end + + +# FUN_CTRAND -- Transform funits +# Error is returned if the transform cannot be made + +procedure fun_ctrand (fun1, fun2, dun, dval, fval1, fval2, nvals) + +pointer fun1 # Input funits pointer +pointer fun2 # Output funits pointer +pointer dun # Input dispersion units pointer +double dval[nvals] # Input dispersion values +double fval1[nvals] # Input flux values +double fval2[nvals] # Output flux values +int nvals # Number of values + +int i +double s, lambda +pointer ang, un_open() +bool fun_compare() +errchk un_open, un_ctrand + +begin + if (fun_compare (fun1, fun2)) { + call amovd (fval1, fval2, nvals) + return + } + + if (FUN_CLASS(fun1) == FUN_UNKNOWN || FUN_CLASS(fun2) == FUN_UNKNOWN) + call error (1, "Cannot convert between selected funits") + + call amovd (fval1, fval2, nvals) + + s = FUN_SCALE(fun1) + switch (FUN_MOD(fun1)) { + case FUN_LOG: + do i = 1, nvals + fval2[i] = 10. ** fval2[i] + case FUN_MAG: + do i = 1, nvals + fval2[i] = 10. ** (-0.4 * fval2[i]) + } + switch (FUN_CLASS(fun1)) { + case FUN_FREQ: + do i = 1, nvals + fval2[i] = fval2[i] / s + case FUN_WAVE: + if (FUN_CLASS(fun2) != FUN_WAVE) { + s = s * FUN_VLIGHT + ang = un_open ("angstroms") + do i = 1, nvals { + call un_ctrand (dun, ang, dval[i], lambda, 1) + fval2[i] = fval2[i] / s * lambda**2 + } + call un_close (ang) + } else { + do i = 1, nvals + fval2[i] = fval2[i] / s + } + } + + s = FUN_SCALE(fun2) + switch (FUN_CLASS(fun2)) { + case FUN_FREQ: + do i = 1, nvals + fval2[i] = fval2[i] * s + case FUN_WAVE: + if (FUN_CLASS(fun1) != FUN_WAVE) { + s = s * FUN_VLIGHT + ang = un_open ("angstroms") + do i = 1, nvals { + call un_ctrand (dun, ang, dval[i], lambda, 1) + fval2[i] = fval2[i] * s / lambda**2 + } + call un_close (ang) + } else { + do i = 1, nvals + fval2[i] = fval2[i] * s + } + } + switch (FUN_MOD(fun2)) { + case FUN_LOG: + do i = 1, nvals + fval2[i] = log10 (fval2[i]) + case FUN_MAG: + do i = 1, nvals + fval2[i] = -2.5 * log10 (fval2[i]) + } + +end + + +# FUN_CHANGED -- Change funits +# Error is returned if the conversion cannot be made + +procedure fun_changed (fun, funits, dun, dvals, fvals, nvals, update) + +pointer fun # Units pointer (may be changed) +char funits[ARB] # Desired funits +pointer dun # Input dispersion pointer +double dvals[nvals] # Input dispersion values +double fvals[nvals] # Flux values +int nvals # Number of values +int update # Update funits pointer? + +bool streq(), fun_compare() +pointer fun1, fun_open() +errchk fun_open, fun_ctrand + +begin + + # Check for same funit string + if (streq (funits, FUN_USER(fun))) + return + + # Check for error in funits string, or the same funits. + fun1 = fun_open (funits) + if (fun_compare (fun1, fun)) { + call strcpy (funits, FUN_USER(fun), SZ_UNITS) + call fun_close (fun1) + return + } + + iferr { + call fun_ctrand (fun, fun1, dun, dvals, fvals, fvals, nvals) + if (update == YES) + call fun_copy (fun1, fun) + call fun_close(fun1) + } then { + call fun_close(fun1) + call erract (EA_ERROR) + } +end diff --git a/noao/onedspec/smw/mkpkg b/noao/onedspec/smw/mkpkg new file mode 100644 index 00000000..64326969 --- /dev/null +++ b/noao/onedspec/smw/mkpkg @@ -0,0 +1,48 @@ +# SMW/SHDR Interface + +update: + $checkout libsmw.a noaolib$ + $update libsmw.a + $checkin libsmw.a noaolib$ + ; + +generic: + $set GEN = "$$generic -k" + + $ifolder (smwctran.x, smwctran.gx) + $(GEN) smwctran.gx -o smwctran.x $endif + ; + +libsmw.a: + $ifeq (USE_GENERIC, yes) $call generic $endif + + funits.x <ctype.h> <error.h> <funits.h> + shdr.x <error.h> <funits.h> <imset.h> <math/iminterp.h>\ + <smw.h> <units.h> <imhdr.h> + smwclose.x <smw.h> + smwct.x <smw.h> + smwctfree.x <smw.h> + smwctran.x <smw.h> + smwdaxis.x <smw.h> + smwequispec.x <mwset.h> <smw.h> <imhdr.h> + smwesms.x <mwset.h> <smw.h> + smwgapid.x <smw.h> + smwgwattrs.x <error.h> <smw.h> + smwmerge.x <mwset.h> <smw.h> + smwmultispec.x <smw.h> + smwmw.x <smw.h> + smwnd.x <imhdr.h> <smw.h> + smwndes.x <imhdr.h> <smw.h> + smwnewcopy.x <smw.h> + smwoldms.x <mwset.h> <smw.h> + smwonedspec.x <smw.h> <imhdr.h> + smwopen.x <smw.h> + smwopenim.x <imio.h> <mwset.h> <imhdr.h> + smwsapid.x <smw.h> + smwsaveim.x <imio.h> <smw.h> <imhdr.h> + smwsaxes.x <imhdr.h> <mwset.h> <smw.h> + smwsctran.x <smw.h> + smwsmw.x <smw.h> + smwswattrs.x <error.h> <smw.h> + units.x <ctype.h> <error.h> <units.h> + ; 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 diff --git a/noao/onedspec/smw/smwclose.x b/noao/onedspec/smw/smwclose.x new file mode 100644 index 00000000..339ebd98 --- /dev/null +++ b/noao/onedspec/smw/smwclose.x @@ -0,0 +1,46 @@ +include <smw.h> + + +# SMW_CLOSE -- Close the SMW data structure. +# This includes closing the MWCS pointers. + +procedure smw_close (smw) + +pointer smw # SMW pointer + +int i +pointer apids + +begin + if (smw == NULL) + return + + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call mfree (SMW_APID(smw), TY_CHAR) + call mw_close (SMW_MW(smw,0)) + case SMW_ES: + call mfree (SMW_APS(smw), TY_INT) + call mfree (SMW_BEAMS(smw), TY_INT) + call mfree (SMW_APLOW(smw), TY_REAL) + call mfree (SMW_APHIGH(smw), TY_REAL) + call mfree (SMW_APID(smw), TY_CHAR) + apids = SMW_APIDS(smw) - 1 + do i = 1, SMW_NSPEC(smw) + call mfree (Memi[apids+i], TY_CHAR) + call mfree (SMW_APIDS(smw), TY_POINTER) + call mw_close (SMW_MW(smw,0)) + case SMW_MS: + call mfree (SMW_APS(smw), TY_INT) + call mfree (SMW_BEAMS(smw), TY_INT) + call mfree (SMW_APLOW(smw), TY_REAL) + call mfree (SMW_APHIGH(smw), TY_REAL) + call mfree (SMW_APID(smw), TY_CHAR) + apids = SMW_APIDS(smw) - 1 + do i = 1, SMW_NSPEC(smw) + call mfree (Memi[apids+i], TY_CHAR) + do i = 0, SMW_NMW(smw)-1 + call mw_close (SMW_MW(smw,i)) + } + call mfree (smw, TY_STRUCT) +end diff --git a/noao/onedspec/smw/smwct.x b/noao/onedspec/smw/smwct.x new file mode 100644 index 00000000..b568f759 --- /dev/null +++ b/noao/onedspec/smw/smwct.x @@ -0,0 +1,19 @@ +include <smw.h> + + +# SMW_CT -- Get MCWS CT pointer for the specified physical line. + +pointer procedure smw_ct (sct, line) + +pointer sct #I SMW pointer +int line #I Physical line + +begin + if (SMW_NCT(sct) == 1) + return (SMW_CT(sct,0)) + + if (line < 1 || line > SMW_NSPEC(SMW_SMW(sct))) + call error (1, "smw_ct: aperture not found") + + return (SMW_CT(sct,(line-1)/SMW_NSPLIT)) +end diff --git a/noao/onedspec/smw/smwctfree.x b/noao/onedspec/smw/smwctfree.x new file mode 100644 index 00000000..90a506d7 --- /dev/null +++ b/noao/onedspec/smw/smwctfree.x @@ -0,0 +1,19 @@ +include <smw.h> + + +# SMW_CTFREE -- Free a spectral SMW coordinate transform pointer. + +procedure smw_ctfree (ct) + +pointer ct # SMW CT pointer +int i + +begin + if (ct == NULL) + return + + do i = 0, SMW_NCT(ct)-1 + call mw_ctfree (SMW_CT(ct,i)) + call mw_ctfree (SMW_CTL(ct)) + call mfree (ct, TY_STRUCT) +end diff --git a/noao/onedspec/smw/smwctran.gx b/noao/onedspec/smw/smwctran.gx new file mode 100644 index 00000000..4029aaab --- /dev/null +++ b/noao/onedspec/smw/smwctran.gx @@ -0,0 +1,166 @@ +include <smw.h> + + +# Evaluate SMW coordinate transform. These procedures call the +# MWCS procedures unless the WCS is a split MULTISPEC WCS. In that +# case the appropriate piece needs to be determined and the physical +# line numbers manipulated. For log sampled spectra conversions +# must be made for EQUISPEC/NDSPEC. The convention is that coordinates +# are always input and output and linear. Note that the MULTISPEC +# function driver already takes care of this. +# +# SMW_CTRANR -- N dimensional real coordinate transformation. +# SMW_CTRAND -- N dimensional double coordinate transformation. +# SMW_C1TRANR -- One dimensional real coordinate transformation. +# SMW_C1TRAND -- One dimensional double coordinate transformation. +# SMW_C2TRANR -- Two dimensional real coordinate transformation. +# SMW_C2TRAND -- Two dimensional double coordinate transformation. + + +$for (rd) +# SMW_CTRAN -- N dimensional coordinate transformation. + +procedure smw_ctran$t (ct, p1, p2, ndim) + +pointer ct #I SMW CT pointer +PIXEL p1[ndim] #I Input coordinate +PIXEL p2[ndim] #O Output coordinate +int ndim #I Dimensionality + +int i, j, format, daxis, aaxis, dtype, naps +pointer smw, aps +errchk mw_ctran$t + +begin + if (SMW_NCT(ct) != 1) + call error (1, "SMW_CTRAN: Wrong WCS type") + + call amov$t (p1, p2, ndim) + + smw = SMW_SMW(ct) + format = SMW_FORMAT(smw) + daxis = SMW_DAXIS(ct) + aaxis = SMW_AAXIS(ct) + dtype = SMW_DTYPE(smw) + naps = SMW_NSPEC(smw) + aps = SMW_APS(smw) + switch (format) { + case SMW_ND, SMW_ES: + switch (SMW_CTTYPE(ct)) { + case SMW_LW, SMW_PW: + call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim) + if (daxis > 0 && dtype == DCLOG) + p2[daxis] = 10. ** max (-20$F, min (20$F, p2[daxis])) + if (aaxis > 0 && format == SMW_ES) { + i = max (1, min (naps, nint (p2[aaxis]))) + p2[aaxis] = Memi[aps+i-1] + } + case SMW_WL, SMW_WP: + if (daxis > 0 && dtype == DCLOG) + p2[daxis] = log10 (p2[daxis]) + if (aaxis > 0 && format == SMW_ES) { + j = nint (p2[aaxis]) + p2[aaxis] = 1 + do i = 1, naps { + if (j == Memi[aps+i-1]) { + p2[aaxis] = i + break + } + } + } + call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim) + default: + call mw_ctran$t (SMW_CT(ct,0), p2, p2, ndim) + } + case SMW_MS: + call mw_ctran$t (SMW_CT(ct,0), p1, p2, ndim) + } +end + + +# SMW_C1TRAN -- One dimensional coordinate transformation. + +PIXEL procedure smw_c1tran$t (ct, x1) + +pointer ct #I SMW CT pointer +PIXEL x1 #I Input coordinate +PIXEL x2 #O Output coordinate + +errchk mw_ctran$t + +begin + call smw_ctran$t (ct, x1, x2, 1) + return (x2) +end + + +# SMW_C2TRAN -- Two dimensional coordinate transformation. + +procedure smw_c2tran$t (ct, x1, y1, x2, y2) + +pointer ct #I SMW CT pointer +PIXEL x1, y1 #I Input coordinates +PIXEL x2, y2 #O Output coordinates + +PIXEL p1[2], p2[2] +int i, j, naps +PIXEL xp, yp +pointer aps, smw_ct() +errchk smw_ct, mw_c2tran$t + +begin + # Unsplit WCS. + if (SMW_NCT(ct) == 1) { + p1[1] = x1 + p1[2] = y1 + call smw_ctran$t (ct, p1, p2, 2) + x2 = p2[1] + y2 = p2[2] + return + } + + # If we get here then we are dealing with a split MULTISPEC WCS. + # Depending on the systems being transformed there may need to + # transformation made on the physical coordinate system. + + switch (SMW_CTTYPE(ct)) { + case SMW_LW: + call mw_c2tran$t (SMW_CTL(ct), x1, y1, xp, yp) + i = nint (yp) + yp = mod (i-1, SMW_NSPLIT) + 1 + call mw_c2tran$t (smw_ct(ct,i), xp, yp, x2, y2) + case SMW_PW: + i = nint (y1) + yp = mod (i-1, SMW_NSPLIT) + 1 + call mw_c2tran$t (smw_ct(ct,i), x1, yp, x2, y2) + case SMW_WL: + aps = SMW_APS(SMW_SMW(ct)) + naps = SMW_NSPEC(SMW_SMW(ct)) + j = nint (y1) + do i = 1, naps { + if (Memi[aps+i-1] == j) { + call mw_c2tran$t (smw_ct(ct,i), x1, y1, xp, yp) + yp = i + call mw_c2tran$t (SMW_CTL(ct), xp, yp, x2, y2) + return + } + } + call error (1, "Aperture not found") + case SMW_WP: + aps = SMW_APS(SMW_SMW(ct)) + naps = SMW_NSPEC(SMW_SMW(ct)) + j = nint (y1) + do i = 1, naps { + if (Memi[aps+i-1] == j) { + call mw_c2tran$t (smw_ct(ct,i), x1, y1, x2, y2) + y2 = i + return + } + } + call error (1, "Aperture not found") + default: + x2 = x1 + y2 = y1 + } +end +$endfor diff --git a/noao/onedspec/smw/smwctran.x b/noao/onedspec/smw/smwctran.x new file mode 100644 index 00000000..38967be2 --- /dev/null +++ b/noao/onedspec/smw/smwctran.x @@ -0,0 +1,312 @@ +include <smw.h> + + +# Evaluate SMW coordinate transform. These procedures call the +# MWCS procedures unless the WCS is a split MULTISPEC WCS. In that +# case the appropriate piece needs to be determined and the physical +# line numbers manipulated. For log sampled spectra conversions +# must be made for EQUISPEC/NDSPEC. The convention is that coordinates +# are always input and output and linear. Note that the MULTISPEC +# function driver already takes care of this. +# +# SMW_CTRANR -- N dimensional real coordinate transformation. +# SMW_CTRAND -- N dimensional double coordinate transformation. +# SMW_C1TRANR -- One dimensional real coordinate transformation. +# SMW_C1TRAND -- One dimensional double coordinate transformation. +# SMW_C2TRANR -- Two dimensional real coordinate transformation. +# SMW_C2TRAND -- Two dimensional double coordinate transformation. + + + +# SMW_CTRAN -- N dimensional coordinate transformation. + +procedure smw_ctranr (ct, p1, p2, ndim) + +pointer ct #I SMW CT pointer +real p1[ndim] #I Input coordinate +real p2[ndim] #O Output coordinate +int ndim #I Dimensionality + +int i, j, format, daxis, aaxis, dtype, naps +pointer smw, aps +errchk mw_ctranr + +begin + if (SMW_NCT(ct) != 1) + call error (1, "SMW_CTRAN: Wrong WCS type") + + call amovr (p1, p2, ndim) + + smw = SMW_SMW(ct) + format = SMW_FORMAT(smw) + daxis = SMW_DAXIS(ct) + aaxis = SMW_AAXIS(ct) + dtype = SMW_DTYPE(smw) + naps = SMW_NSPEC(smw) + aps = SMW_APS(smw) + switch (format) { + case SMW_ND, SMW_ES: + switch (SMW_CTTYPE(ct)) { + case SMW_LW, SMW_PW: + call mw_ctranr (SMW_CT(ct,0), p2, p2, ndim) + if (daxis > 0 && dtype == DCLOG) + p2[daxis] = 10. ** max (-20.0, min (20.0, p2[daxis])) + if (aaxis > 0 && format == SMW_ES) { + i = max (1, min (naps, nint (p2[aaxis]))) + p2[aaxis] = Memi[aps+i-1] + } + case SMW_WL, SMW_WP: + if (daxis > 0 && dtype == DCLOG) + p2[daxis] = log10 (p2[daxis]) + if (aaxis > 0 && format == SMW_ES) { + j = nint (p2[aaxis]) + p2[aaxis] = 1 + do i = 1, naps { + if (j == Memi[aps+i-1]) { + p2[aaxis] = i + break + } + } + } + call mw_ctranr (SMW_CT(ct,0), p2, p2, ndim) + default: + call mw_ctranr (SMW_CT(ct,0), p2, p2, ndim) + } + case SMW_MS: + call mw_ctranr (SMW_CT(ct,0), p1, p2, ndim) + } +end + + +# SMW_C1TRAN -- One dimensional coordinate transformation. + +real procedure smw_c1tranr (ct, x1) + +pointer ct #I SMW CT pointer +real x1 #I Input coordinate +real x2 #O Output coordinate + +errchk mw_ctranr + +begin + call smw_ctranr (ct, x1, x2, 1) + return (x2) +end + + +# SMW_C2TRAN -- Two dimensional coordinate transformation. + +procedure smw_c2tranr (ct, x1, y1, x2, y2) + +pointer ct #I SMW CT pointer +real x1, y1 #I Input coordinates +real x2, y2 #O Output coordinates + +real p1[2], p2[2] +int i, j, naps +real xp, yp +pointer aps, smw_ct() +errchk smw_ct, mw_c2tranr + +begin + # Unsplit WCS. + if (SMW_NCT(ct) == 1) { + p1[1] = x1 + p1[2] = y1 + call smw_ctranr (ct, p1, p2, 2) + x2 = p2[1] + y2 = p2[2] + return + } + + # If we get here then we are dealing with a split MULTISPEC WCS. + # Depending on the systems being transformed there may need to + # transformation made on the physical coordinate system. + + switch (SMW_CTTYPE(ct)) { + case SMW_LW: + call mw_c2tranr (SMW_CTL(ct), x1, y1, xp, yp) + i = nint (yp) + yp = mod (i-1, SMW_NSPLIT) + 1 + call mw_c2tranr (smw_ct(ct,i), xp, yp, x2, y2) + case SMW_PW: + i = nint (y1) + yp = mod (i-1, SMW_NSPLIT) + 1 + call mw_c2tranr (smw_ct(ct,i), x1, yp, x2, y2) + case SMW_WL: + aps = SMW_APS(SMW_SMW(ct)) + naps = SMW_NSPEC(SMW_SMW(ct)) + j = nint (y1) + do i = 1, naps { + if (Memi[aps+i-1] == j) { + call mw_c2tranr (smw_ct(ct,i), x1, y1, xp, yp) + yp = i + call mw_c2tranr (SMW_CTL(ct), xp, yp, x2, y2) + return + } + } + call error (1, "Aperture not found") + case SMW_WP: + aps = SMW_APS(SMW_SMW(ct)) + naps = SMW_NSPEC(SMW_SMW(ct)) + j = nint (y1) + do i = 1, naps { + if (Memi[aps+i-1] == j) { + call mw_c2tranr (smw_ct(ct,i), x1, y1, x2, y2) + y2 = i + return + } + } + call error (1, "Aperture not found") + default: + x2 = x1 + y2 = y1 + } +end + +# SMW_CTRAN -- N dimensional coordinate transformation. + +procedure smw_ctrand (ct, p1, p2, ndim) + +pointer ct #I SMW CT pointer +double p1[ndim] #I Input coordinate +double p2[ndim] #O Output coordinate +int ndim #I Dimensionality + +int i, j, format, daxis, aaxis, dtype, naps +pointer smw, aps +errchk mw_ctrand + +begin + if (SMW_NCT(ct) != 1) + call error (1, "SMW_CTRAN: Wrong WCS type") + + call amovd (p1, p2, ndim) + + smw = SMW_SMW(ct) + format = SMW_FORMAT(smw) + daxis = SMW_DAXIS(ct) + aaxis = SMW_AAXIS(ct) + dtype = SMW_DTYPE(smw) + naps = SMW_NSPEC(smw) + aps = SMW_APS(smw) + switch (format) { + case SMW_ND, SMW_ES: + switch (SMW_CTTYPE(ct)) { + case SMW_LW, SMW_PW: + call mw_ctrand (SMW_CT(ct,0), p2, p2, ndim) + if (daxis > 0 && dtype == DCLOG) + p2[daxis] = 10. ** max (-20.0D0, min (20.0D0, p2[daxis])) + if (aaxis > 0 && format == SMW_ES) { + i = max (1, min (naps, nint (p2[aaxis]))) + p2[aaxis] = Memi[aps+i-1] + } + case SMW_WL, SMW_WP: + if (daxis > 0 && dtype == DCLOG) + p2[daxis] = log10 (p2[daxis]) + if (aaxis > 0 && format == SMW_ES) { + j = nint (p2[aaxis]) + p2[aaxis] = 1 + do i = 1, naps { + if (j == Memi[aps+i-1]) { + p2[aaxis] = i + break + } + } + } + call mw_ctrand (SMW_CT(ct,0), p2, p2, ndim) + default: + call mw_ctrand (SMW_CT(ct,0), p2, p2, ndim) + } + case SMW_MS: + call mw_ctrand (SMW_CT(ct,0), p1, p2, ndim) + } +end + + +# SMW_C1TRAN -- One dimensional coordinate transformation. + +double procedure smw_c1trand (ct, x1) + +pointer ct #I SMW CT pointer +double x1 #I Input coordinate +double x2 #O Output coordinate + +errchk mw_ctrand + +begin + call smw_ctrand (ct, x1, x2, 1) + return (x2) +end + + +# SMW_C2TRAN -- Two dimensional coordinate transformation. + +procedure smw_c2trand (ct, x1, y1, x2, y2) + +pointer ct #I SMW CT pointer +double x1, y1 #I Input coordinates +double x2, y2 #O Output coordinates + +double p1[2], p2[2] +int i, j, naps +double xp, yp +pointer aps, smw_ct() +errchk smw_ct, mw_c2trand + +begin + # Unsplit WCS. + if (SMW_NCT(ct) == 1) { + p1[1] = x1 + p1[2] = y1 + call smw_ctrand (ct, p1, p2, 2) + x2 = p2[1] + y2 = p2[2] + return + } + + # If we get here then we are dealing with a split MULTISPEC WCS. + # Depending on the systems being transformed there may need to + # transformation made on the physical coordinate system. + + switch (SMW_CTTYPE(ct)) { + case SMW_LW: + call mw_c2trand (SMW_CTL(ct), x1, y1, xp, yp) + i = nint (yp) + yp = mod (i-1, SMW_NSPLIT) + 1 + call mw_c2trand (smw_ct(ct,i), xp, yp, x2, y2) + case SMW_PW: + i = nint (y1) + yp = mod (i-1, SMW_NSPLIT) + 1 + call mw_c2trand (smw_ct(ct,i), x1, yp, x2, y2) + case SMW_WL: + aps = SMW_APS(SMW_SMW(ct)) + naps = SMW_NSPEC(SMW_SMW(ct)) + j = nint (y1) + do i = 1, naps { + if (Memi[aps+i-1] == j) { + call mw_c2trand (smw_ct(ct,i), x1, y1, xp, yp) + yp = i + call mw_c2trand (SMW_CTL(ct), xp, yp, x2, y2) + return + } + } + call error (1, "Aperture not found") + case SMW_WP: + aps = SMW_APS(SMW_SMW(ct)) + naps = SMW_NSPEC(SMW_SMW(ct)) + j = nint (y1) + do i = 1, naps { + if (Memi[aps+i-1] == j) { + call mw_c2trand (smw_ct(ct,i), x1, y1, x2, y2) + y2 = i + return + } + } + call error (1, "Aperture not found") + default: + x2 = x1 + y2 = y1 + } +end + diff --git a/noao/onedspec/smw/smwdaxis.x b/noao/onedspec/smw/smwdaxis.x new file mode 100644 index 00000000..0bea9375 --- /dev/null +++ b/noao/onedspec/smw/smwdaxis.x @@ -0,0 +1,109 @@ +include <smw.h> + +define CTYPES "|LAMBDA|FREQ|WAVELENGTH|VELO|VELO-LSR|VELO-HEL|VELO-OBS|" + +# SMW_DAXIS -- Set physical dispersion axis and summing factors. +# A default value of zero for the dispersion axis will cause the dispersion +# axis to be sought in the image header and, if not found, from the CL +# "dispaxis" parameter. A default value of zero for the summing factors will +# cause them to be queried from the CL "nsum" parameter. A default value of +# INDEFI in either parameter will leave the current default unchanged. +# +# When this procedure is called with an SMW and IMIO pointer the SMW +# pointer is updated to desired default dispersion axis and summing +# parameters. + +procedure smw_daxis (smw, im, daxisp, nsum1, nsum2) + +pointer smw #I SMW pointer +pointer im #I IMIO pointer +int daxisp #I Default dispersion axis +int nsum1, nsum2 #I Default summing factors + +int i, da, ns[2] +int imgeti(), clgeti(), clscan(), nscan(), nowhite(), strdic() +pointer sp, key, val +data da/0/, ns/0,0/ +errchk clgeti + +begin + # Set defaults. + # A value of 0 will use the image DISPAXIS or query the CL and + # a value of INDEFI will leave the current default unchanged. + + if (!IS_INDEFI (daxisp)) + da = daxisp + if (!IS_INDEFI (nsum1)) + ns[1] = nsum1 + if (!IS_INDEFI (nsum2)) + ns[2] = nsum2 + + if (smw == NULL) + return + + # This procedure is specific to the NDSPEC format. + if (SMW_FORMAT(smw) != SMW_ND) + return + + # Set dispersion axis. + if (da == 0) { + if (im == NULL) + SMW_PAXIS(smw,1) = clgeti ("dispaxis") + else { + iferr (SMW_PAXIS(smw,1) = imgeti (im, "DISPAXIS")) { + SMW_PAXIS(smw,1) = clgeti ("dispaxis") + call smark (sp) + call salloc (key, 8, TY_CHAR) + call salloc (val, SZ_FNAME, TY_CHAR) + do i = 1, 7 { + call sprintf (Memc[key], 8, "CTYPE%d") + call pargi (i) + iferr (call imgstr (im, Memc[key], Memc[val], SZ_FNAME)) + break + if (nowhite (Memc[val], Memc[val], SZ_FNAME) > 0) { + call strupr (Memc[val]) + if (strdic(Memc[val],Memc[val],SZ_FNAME,CTYPES)>0) { + SMW_PAXIS(smw,1) = i + break + } + } + } + call sfree (sp) + } + if (SMW_PAXIS(smw,1) < 1 || SMW_PAXIS(smw,1) > 7) { + i = SMW_PAXIS(smw,1) + SMW_PAXIS(smw,1) = clgeti ("dispaxis") + call eprintf ( + "WARNING: Image header dispersion axis %d invalid. Using axis %d.\n") + call pargi (i) + call pargi (SMW_PAXIS(smw,1)) + } + } + } else + SMW_PAXIS(smw,1) = da + + # Set summing parameters. + if (ns[1] == 0 || ns[2] == 0) { + if (clscan("nsum") == EOF) + call error (1, "smw_daxis: Error in 'nsum' parameter") + call gargi (i) + if (ns[1] == 0) { + if (nscan() == 1) + SMW_NSUM(smw,1) = max (1, i) + else + call error (1, "smw_daxis: Error in 'nsum' parameter") + } else + SMW_NSUM(smw,1) = ns[1] + call gargi (i) + if (ns[2] == 0) { + if (nscan() == 2) + SMW_NSUM(smw,2) = max (1, i) + else + SMW_NSUM(smw,2) = SMW_NSUM(smw,1) + } else + SMW_NSUM(smw,2) = ns[2] + } else { + SMW_NSUM(smw,1) = ns[1] + SMW_NSUM(smw,2) = ns[2] + } +end diff --git a/noao/onedspec/smw/smwequispec.x b/noao/onedspec/smw/smwequispec.x new file mode 100644 index 00000000..5cda57c0 --- /dev/null +++ b/noao/onedspec/smw/smwequispec.x @@ -0,0 +1,86 @@ +include <imhdr.h> +include <mwset.h> +include <smw.h> + + +# SMW_EQUISPEC -- Setup the EQUISPEC SMW parameters. +# The aperture information is in the APNUM and APID keywords and the +# WCS information is in the linear MWCS. + +procedure smw_equispec (im, smw) + +pointer im #I IMIO pointer +pointer smw #U MWCS pointer input SMW pointer output + +int i, j, k, nchar, ap, beam, dtype, nw +double w1, dw, z +real aplow[2], aphigh[2], mw_c1tranr() +pointer sp, key, val, wterm, mw, ct, mw_sctran() +int imgeti(), mw_stati(), ctoi(), ctor() +errchk imgstr, mw_gwtermd, mw_sctran +errchk smw_open, smw_saxes, smw_swattrs, smw_sapid + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (val, SZ_LINE, TY_CHAR) + + # Determine the dispersion parameters + mw = smw + i = mw_stati (mw, MW_NDIM) + call salloc (wterm, 2*i+i*i, TY_DOUBLE) + call mw_gwtermd (mw, Memd[wterm], Memd[wterm+i], Memd[wterm+2*i], i) + w1 = Memd[wterm+i] + dw = Memd[wterm+2*i] + + # Determine the number of physical pixels. + ct = mw_sctran (mw, "logical", "physical", 1) + nw = max (mw_c1tranr (ct, 1.), mw_c1tranr (ct, real(IM_LEN(im,1)))) + call mw_ctfree (ct) + + # Determine the dispersion type. + iferr (dtype = imgeti (im, "DC-FLAG")) + dtype = DCNO + if (dtype==DCLOG) { + if (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.) + dtype = DCLINEAR + else { + w1 = 10D0 ** w1 + dw = w1 * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1) + } + } + + # Set the SMW data structure. + call smw_open (smw, NULL, im) + do i = 1, SMW_NSPEC(smw) { + call smw_mw (smw, i, 1, mw, j, k) + if (j < 1000) + call sprintf (Memc[key], SZ_FNAME, "APNUM%d") + else + call sprintf (Memc[key], SZ_FNAME, "AP%d") + call pargi (j) + call imgstr (im, Memc[key], Memc[val], SZ_LINE) + k = 1 + nchar = ctoi (Memc[val], k, ap) + nchar = ctoi (Memc[val], k, beam) + if (ctor (Memc[val], k, aplow[1]) == 0) + aplow[1] = INDEF + if (ctor (Memc[val], k, aphigh[1]) == 0) + aphigh[1] = INDEF + if (ctor (Memc[val], k, aplow[2]) == 0) + aplow[2] = INDEF + if (ctor (Memc[val], k, aphigh[2]) == 0) + aphigh[2] = INDEF + z = 0. + + call smw_swattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, "") + + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (j) + ifnoerr (call imgstr (im, Memc[key], Memc[val], SZ_LINE)) + call smw_sapid (smw, i, 1, Memc[val]) + } + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwesms.x b/noao/onedspec/smw/smwesms.x new file mode 100644 index 00000000..d1b8a368 --- /dev/null +++ b/noao/onedspec/smw/smwesms.x @@ -0,0 +1,96 @@ +include <mwset.h> +include <smw.h> + + +# SMW_ESMS -- Convert EQUISPEC WCS into MULTISPEC WCS. +# This is called by SMW_SWATTRS when the equal linear coordinate system +# requirement of the EQUISPEC WCS is violated. + +procedure smw_esms (smw) + +pointer smw #U SMW pointer + +int i, j, k, pdim1, pdim2, ap, beam, dtype, nw, axes[2] +double w1, dw, z +real aplow, aphigh +pointer sp, key, str, lterm, mw, mw1, mw2, apid, mw_open() +data axes/1,2/ + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (lterm, 12, TY_DOUBLE) + + # Set the basic MWCS + mw1 = SMW_MW(smw,0) + pdim1 = SMW_PDIM(smw) + pdim2 = max (2, pdim1) + mw2 = mw_open (NULL, pdim2) + call mw_newsystem (mw2, "multispec", pdim2) + call mw_swtype (mw2, axes, 2, "multispec", "") + if (pdim2 > 2) + call mw_swtype (mw2, 3, 1, "linear", "") + call mw_gltermd (mw1, Memd[lterm+pdim2], Memd[lterm], pdim1) + if (pdim1 == 1) { + Memd[lterm+1] = 0. + Memd[lterm+3] = 0. + Memd[lterm+4] = 0. + Memd[lterm+5] = 1. + } + call mw_sltermd (mw2, Memd[lterm+pdim2], Memd[lterm], pdim2) + ifnoerr (call mw_gwattrs (mw1, 1, "label", Memc[str], SZ_LINE)) + call mw_swattrs (mw2, 1, "label", Memc[str]) + ifnoerr (call mw_gwattrs (mw1, 1, "units", Memc[str], SZ_LINE)) + call mw_swattrs (mw2, 1, "units", Memc[str]) + ifnoerr (call mw_gwattrs (mw1, 1, "units_display", Memc[str], SZ_LINE)) + call mw_swattrs (mw2, 1, "units_display", Memc[str]) + + # Write the MULTISPEC WCS + dtype = SMW_DTYPE(smw) + w1 = SMW_W1(smw) + dw = SMW_DW(smw) + nw = SMW_NW(smw) + z = SMW_Z(smw) + if (dtype == DCLOG) { + dw = log10 ((w1+(nw-1)*dw)/w1)/(nw-1) + w1 = log10 (w1) + } + + call smw_open (mw2, smw, 0) + do i = 1, SMW_NSPEC(smw) { + ap = Memi[SMW_APS(smw)+i-1] + beam = Memi[SMW_BEAMS(smw)+i-1] + aplow = Memr[SMW_APLOW(smw)+2*i-2] + aphigh = Memr[SMW_APHIGH(smw)+2*i-2] + apid = Memi[SMW_APIDS(smw)+i-1] + + call smw_mw (mw2, i, 1, mw, j, k) + call sprintf (Memc[key], SZ_FNAME, "spec%d") + call pargi (j) + call sprintf (Memc[str], SZ_LINE, + "%d %d %d %.14g %.14g %d %.14g %.2f %.2f") + call pargi (ap) + call pargi (beam) + call pargi (dtype) + call pargd (w1) + call pargd (dw) + call pargi (nw) + call pargd (z) + call pargr (aplow) + call pargr (aphigh) + call mw_swattrs (mw, 2, Memc[key], Memc[str]) + + if (apid != NULL) + call smw_sapid (mw2, i, 1, Memc[apid]) + + # This is used if there is a split MULTISPEC WCS. + if (SMW_APS(mw2) != NULL) + Memi[SMW_APS(mw2)+i-1] = ap + } + + call smw_close (smw) + smw = mw2 + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwgapid.x b/noao/onedspec/smw/smwgapid.x new file mode 100644 index 00000000..214be533 --- /dev/null +++ b/noao/onedspec/smw/smwgapid.x @@ -0,0 +1,30 @@ +include <smw.h> + + +# SMW_GAPID -- Get aperture id. + +procedure smw_gapid (smw, index1, index2, apid, maxchar) + +pointer smw #I SMW pointer +int index1 #I Spectrum index +int index2 #I Spectrum index +char apid[maxchar] #O Aperture id +int maxchar #I Maximum number of characters + +pointer ptr + +begin + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call strcpy (Memc[SMW_APID(smw)], apid, maxchar) + case SMW_ES, SMW_MS: + if (index1 < 0 || index1 > SMW_NSPEC(smw)) + call error (1, "smw_gapid: index out of bounds") + + ptr = Memi[SMW_APIDS(smw)+index1-1] + if (index1 == 0 || ptr == NULL) + call strcpy (Memc[SMW_APID(smw)], apid, maxchar) + else + call strcpy (Memc[ptr], apid, maxchar) + } +end diff --git a/noao/onedspec/smw/smwgwattrs.x b/noao/onedspec/smw/smwgwattrs.x new file mode 100644 index 00000000..4084fd4c --- /dev/null +++ b/noao/onedspec/smw/smwgwattrs.x @@ -0,0 +1,134 @@ +include <error.h> +include <smw.h> + + +# SMW_GWATTRS -- Get spectrum attribute parameters. +# BE CAREFUL OF OUTPUT VARIABLES BEING THE SAME MEMORY ADDRESS! + +procedure smw_gwattrs (smw, index1, index2, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + +pointer smw # SMW pointer +int index1 # Spectrum index +int index2 # Spectrum index +int ap # Aperture number +int beam # Beam number +int dtype # Dispersion type +double w1 # Starting coordinate +double dw # Coordinate interval +int nw # Number of valid pixels +double z # Redshift factor +real aplow[2], aphigh[2] # Aperture limits +pointer coeff # Nonlinear coeff string (input/output) + +int i, j, n, ip, sz_coeff, strlen(), ctoi(), ctor(), ctod() +double a, b +pointer sp, key, mw +errchk smw_mw, mw_gwattrs + +data sz_coeff /SZ_LINE/ + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + + if (coeff == NULL) + call malloc (coeff, sz_coeff, TY_CHAR) + else + call realloc (coeff, sz_coeff, TY_CHAR) + + # Determine parameters based on the SMW format. + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call smw_mw (smw, index1, index2, mw, i, j) + + dtype = SMW_DTYPE(smw) + nw = SMW_NW(smw) + w1 = SMW_W1(smw) + dw = SMW_DW(smw) + z = SMW_Z(smw) + + ap = index1 + beam = 0 + aplow[1] = 1 + aphigh[1] = 1 + aplow[2] = 1 + aphigh[2] = 1 + if (SMW_LDIM(smw) > 1) { + aplow[1] = i - (SMW_NSUM(smw,1)-1) / 2 + aphigh[1] = nint (aplow[1]) + SMW_NSUM(smw,1) - 1 + aplow[1] = max (1, nint (aplow[1])) + aphigh[1] = min (SMW_LLEN(smw,2), nint (aphigh[1])) + } + if (SMW_LDIM(smw) > 2) { + aplow[2] = j - (SMW_NSUM(smw,2)-1) / 2 + aphigh[2] = nint (aplow[2]) + SMW_NSUM(smw,2) - 1 + aplow[2] = max (1, nint (aplow[2])) + aphigh[2] = min (SMW_LLEN(smw,3), nint (aphigh[2])) + } + + Memc[coeff] = EOS + case SMW_ES: + call smw_mw (smw, index1, index2, mw, i, j) + + dtype = SMW_DTYPE(smw) + nw = SMW_NW(smw) + w1 = SMW_W1(smw) + dw = SMW_DW(smw) + z = SMW_Z(smw) + + ap = Memi[SMW_APS(smw)+index1-1] + beam = Memi[SMW_BEAMS(smw)+index1-1] + aplow[1] = Memr[SMW_APLOW(smw)+2*index1-2] + aphigh[1] = Memr[SMW_APHIGH(smw)+2*index1-2] + aplow[2] = Memr[SMW_APLOW(smw)+2*index1-1] + aphigh[2] = Memr[SMW_APHIGH(smw)+2*index1-1] + + Memc[coeff] = EOS + case SMW_MS: + call smw_mw (smw, index1, index2, mw, i, j) + + call sprintf (Memc[key], SZ_FNAME, "spec%d") + call pargi (i) + + call mw_gwattrs (mw, 2, Memc[key], Memc[coeff], sz_coeff) + while (strlen (Memc[coeff]) == sz_coeff) { + sz_coeff = 2 * sz_coeff + call realloc (coeff, sz_coeff, TY_CHAR) + call mw_gwattrs (mw, 2, Memc[key], Memc[coeff], sz_coeff) + } + + ip = 1 + i = ctoi (Memc[coeff], ip, ap) + i = ctoi (Memc[coeff], ip, beam) + i = ctoi (Memc[coeff], ip, j) + i = ctod (Memc[coeff], ip, a) + i = ctod (Memc[coeff], ip, b) + i = ctoi (Memc[coeff], ip, n) + i = ctod (Memc[coeff], ip, z) + i = ctor (Memc[coeff], ip, aplow[1]) + i = ctor (Memc[coeff], ip, aphigh[1]) + aplow[2] = INDEF + aphigh[2] = INDEF + if (Memc[coeff+ip-1] != EOS) + call strcpy (Memc[coeff+ip], Memc[coeff], sz_coeff) + else + Memc[coeff] = EOS + + if (j==DCLOG) { + if (abs(a)>20. || abs(a+(n-1)*b)>20.) + j = DCLINEAR + else { + a = 10D0 ** a + b = a * (10D0 ** ((n-1)*b) - 1) / (n - 1) + } + } + + dtype = j + w1 = a + dw = b + nw = n + } + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwmerge.x b/noao/onedspec/smw/smwmerge.x new file mode 100644 index 00000000..d3e09bd1 --- /dev/null +++ b/noao/onedspec/smw/smwmerge.x @@ -0,0 +1,102 @@ +include <mwset.h> +include <smw.h> + + +# SMW_MERGE -- Merge split MWCS array to a single MWCS. + +procedure smw_merge (smw) + +pointer smw #U Input split WCS, output single WCS + +int i, pdim, naps, format, beam, dtype, dtype1, nw, nw1 +int ap, axes[3] +double w1, dw, z, w11, dw1, z1 +real aplow[2], aphigh[2] +pointer sp, key, val, term, coeff, mw, mw1, mw_open() +data axes/1,2,3/ + +begin + if (SMW_NMW(smw) == 1) + return + + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (val, SZ_LINE, TY_CHAR) + call salloc (term, 15, TY_DOUBLE) + coeff = NULL + + pdim = SMW_PDIM(smw) + naps = SMW_NSPEC(smw) + mw1 = SMW_MW(smw,0) + + # Determine output WCS format. + format = SMW_ES + do i = 1, naps { + call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + if (i == 1) { + dtype1 = dtype + w11 = w1 + dw1 = dw + z1 = z + nw1 = nw + } + if (dtype>1||dtype!=dtype1||w1!=w11||dw!=dw1||nw!=nw1||z!=z1) { + format = SMW_MS + break + } + } + + # Setup WCS. + switch (format) { + case SMW_ES: + mw = mw_open (NULL, pdim) + call mw_newsystem (mw, "equispec", pdim) + call mw_swtype (mw, axes, pdim, "linear", "") + + case SMW_MS: + mw = mw_open (NULL, pdim) + call mw_newsystem (mw, "multispec", pdim) + call mw_swtype (mw, axes, pdim, "multispec", "") + if (pdim > 2) + call mw_swtype (mw, 3, 1, "linear", "") + } + + ifnoerr (call mw_gwattrs (mw1, 1, "label", Memc[val], SZ_LINE)) + call mw_swattrs (mw, 1, "label", Memc[val]) + ifnoerr (call mw_gwattrs (mw1, 1, "units", Memc[val], SZ_LINE)) + call mw_swattrs (mw, 1, "units", Memc[val]) + ifnoerr (call mw_gwattrs (mw1, 1, "units_display", Memc[val], SZ_LINE)) + call mw_swattrs (mw, 1, "units_display", Memc[val]) + call mw_gltermd (mw1, Memd[term+pdim], Memd[term], pdim) + call mw_sltermd (mw, Memd[term+pdim], Memd[term], pdim) + call mw_gwtermd (mw1, Memd[term], Memd[term+pdim], + Memd[term+2*pdim], pdim) + Memd[term] = 1. + Memd[term+pdim] = w1 / (1 + z) + Memd[term+2*pdim] = dw / (1 + z) + call mw_swtermd (mw, Memd[term], Memd[term+pdim], + Memd[term+2*pdim], pdim) + + # Set the SMW structure. + call smw_open (mw, smw, NULL) + if (format == SMW_MS) { + do i = 1, SMW_NMW(mw) - 1 + call mw_close (SMW_MW(mw,i)) + SMW_NMW(mw) = 1 + } + do i = 1, naps { + call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + call smw_swattrs (mw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, Memc[coeff]) + call smw_gapid (smw, i, 1, Memc[val], SZ_LINE) + call smw_sapid (mw, i, 1, Memc[val]) + } + + call smw_close (smw) + smw = mw + + call mfree (coeff, TY_CHAR) + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwmultispec.x b/noao/onedspec/smw/smwmultispec.x new file mode 100644 index 00000000..18e2dbd0 --- /dev/null +++ b/noao/onedspec/smw/smwmultispec.x @@ -0,0 +1,30 @@ +include <smw.h> + + +# SMW_MULTISPEC -- Setup the MULTISPEC SMW parameters. + +procedure smw_multispec (im, smw) + +pointer im #I IMIO pointer +pointer smw #U MWCS pointer input SMW pointer output + +int i, j, k +pointer sp, key, val, mw +errchk smw_open, smw_saxes, smw_sapid + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (val, SZ_LINE, TY_CHAR) + + call smw_open (smw, NULL, im) + do i = 1, SMW_NSPEC(smw) { + call smw_mw (smw, i, 1, mw, j, k) + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (j) + ifnoerr (call imgstr (im, Memc[key], Memc[val], SZ_LINE)) + call smw_sapid (smw, i, 1, Memc[val]) + } + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwmw.x b/noao/onedspec/smw/smwmw.x new file mode 100644 index 00000000..a79aaf98 --- /dev/null +++ b/noao/onedspec/smw/smwmw.x @@ -0,0 +1,38 @@ +include <smw.h> + + +# SMW_MW -- Get MWCS pointer and coordinates from spectrum line and band + +procedure smw_mw (smw, line, band, mw, x, y) + +pointer smw #I SMW pointer +int line #I Spectrum line +int band #I Spectrum band +pointer mw #O MWCS pointer +int x, y #O MWCS coordinates + +real mw_c1tranr() + +begin + if (line < 1 || line > SMW_NSPEC(smw)) + call error (1, "smw_mw: spectrum not found") + + switch (SMW_FORMAT(smw)) { + case SMW_ND: + mw = SMW_MW(smw,0) + x = mod (line - 1, SMW_LLEN(smw,2)) + 1 + y = (line - 1) / SMW_LLEN(smw,2) + band + default: + if (SMW_NMW(smw) == 1) { + mw = SMW_MW(smw,0) + x = line + y = band + if (SMW_CTLP(smw) != NULL) + x = nint (mw_c1tranr (SMW_CTLP(smw), real(line))) + } else { + mw = SMW_MW(smw,(line-1)/SMW_NSPLIT) + x = mod (line - 1, SMW_NSPLIT) + 1 + y = band + } + } +end diff --git a/noao/onedspec/smw/smwnd.x b/noao/onedspec/smw/smwnd.x new file mode 100644 index 00000000..10b48079 --- /dev/null +++ b/noao/onedspec/smw/smwnd.x @@ -0,0 +1,19 @@ +include <imhdr.h> +include <smw.h> + + +# SMW_ND -- Setup the NDSPEC SMW. +# If there is only one spectrum convert it to EQUISPEC if possible. + +procedure smw_nd (im, smw) + +pointer im #I IMIO pointer +pointer smw #U MWCS pointer input SMW pointer output + +errchk smw_open, smw_daxis, smw_ndes + +begin + call smw_open (smw, NULL, im) + if (SMW_NSPEC(smw) == 1) + call smw_ndes (im, smw) +end diff --git a/noao/onedspec/smw/smwndes.x b/noao/onedspec/smw/smwndes.x new file mode 100644 index 00000000..9d35ca6d --- /dev/null +++ b/noao/onedspec/smw/smwndes.x @@ -0,0 +1,82 @@ +include <imhdr.h> +include <smw.h> + + +# SMW_NDES -- Convert NDSPEC WCS into EQUISPEC WCS. +# This requires that the logical dispersion axis be 1. + +procedure smw_ndes (im, smw) + +pointer im #I IMIO pointer +pointer smw #U Input NDSPEC SMW, output EQUISPEC SMW + +int i, pdim1, pdim2, daxis, ap, beam, dtype, nw, axes[2] +real aplow[2], aphigh[2] +double w1, dw, z +pointer sp, key, str, lterm1, lterm2, coeff, mw1, mw2, mw_open() +errchk mw_open, mw_gltermd, mw_gwtermd, smw_open, smw_saxes, smw_gwattrs +data axes/1,2/, coeff/NULL/ + +begin + # Require the dispersion to be along the first logical axis. + if (SMW_LAXIS(smw,1) != 1) + return + + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (lterm1, 15, TY_DOUBLE) + call salloc (lterm2, 15, TY_DOUBLE) + + # Set the MWCS. Only the logical and world transformations along + # the dispersion axis are transfered. + + pdim1 = SMW_PDIM(smw) + pdim2 = IM_NDIM(im) + daxis = SMW_PAXIS(smw,1) + mw1 = SMW_MW(smw,0) + + mw2 = mw_open (NULL, pdim2) + call mw_newsystem (mw2, "equispec", pdim2) + call mw_swtype (mw2, axes, pdim2, "linear", "") + ifnoerr (call mw_gwattrs (mw1, daxis, "label", Memc[str], SZ_LINE)) + call mw_swattrs (mw2, 1, "label", Memc[str]) + ifnoerr (call mw_gwattrs (mw1, daxis, "units", Memc[str], SZ_LINE)) + call mw_swattrs (mw2, 1, "units", Memc[str]) + ifnoerr (call mw_gwattrs (mw1, daxis, "units_display", Memc[str], + SZ_LINE)) + call mw_swattrs (mw2, 1, "units_display", Memc[str]) + + call mw_gltermd (mw1, Memd[lterm1+pdim1], Memd[lterm1], pdim1) + call mw_gltermd (mw2, Memd[lterm2+pdim2], Memd[lterm2], pdim2) + Memd[lterm2] = Memd[lterm1+daxis-1] + Memd[lterm2+pdim2] = Memd[lterm1+pdim1+(pdim1+1)*(daxis-1)] + call mw_sltermd (mw2, Memd[lterm2+pdim2], Memd[lterm2], pdim2) + + call mw_gwtermd (mw1, Memd[lterm1], Memd[lterm1+pdim1], + Memd[lterm1+2*pdim1], pdim1) + call mw_gwtermd (mw2, Memd[lterm2], Memd[lterm2+pdim2], + Memd[lterm2+2*pdim2], pdim2) + Memd[lterm2] = Memd[lterm1+daxis-1] + Memd[lterm2+pdim2] = Memd[lterm1+pdim1+daxis-1] + Memd[lterm2+2*pdim2] = Memd[lterm1+2*pdim1+(pdim1+1)*(daxis-1)] + call mw_swtermd (mw2, Memd[lterm2], Memd[lterm2+pdim2], + Memd[lterm2+2*pdim2], pdim2) + + # Set the EQUISPEC SMW. + IM_LEN(im,2) = SMW_NSPEC(smw) + IM_LEN(im,3) = SMW_NBANDS(smw) + call smw_open (mw2, NULL, im) + do i = 1, SMW_NSPEC(smw) { + call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + call smw_swattrs (mw2, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, Memc[coeff]) + } + call mfree (coeff, TY_CHAR) + + call smw_close (smw) + smw = mw2 + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwnewcopy.x b/noao/onedspec/smw/smwnewcopy.x new file mode 100644 index 00000000..230ed205 --- /dev/null +++ b/noao/onedspec/smw/smwnewcopy.x @@ -0,0 +1,58 @@ +include <smw.h> + + +# SMW_NEWCOPY -- Make a new copy of an SMW structure. + +pointer procedure smw_newcopy (smw) + +pointer smw #I SMW pointer to copy +pointer new #O SMW copy + +int i, nspec +pointer mw_newcopy(), mw_sctran() + +begin + call calloc (new, SMW_LEN(SMW_NMW(smw)), TY_STRUCT) + call amovi (Memi[smw], Memi[new], SMW_LEN(SMW_NMW(smw))) + + if (SMW_APID(smw) != NULL) { + call malloc (SMW_APID(new), SZ_LINE, TY_CHAR) + call strcpy (Memc[SMW_APID(smw)], Memc[SMW_APID(new)], SZ_LINE) + } + + nspec = SMW_NSPEC(smw) + if (SMW_APS(smw) != NULL) { + call malloc (SMW_APS(new), nspec, TY_INT) + call amovi (Memi[SMW_APS(smw)], Memi[SMW_APS(new)], nspec) + } + if (SMW_BEAMS(smw) != NULL) { + call malloc (SMW_BEAMS(new), nspec, TY_INT) + call amovi (Memi[SMW_BEAMS(smw)], Memi[SMW_BEAMS(new)], nspec) + } + if (SMW_APLOW(smw) != NULL) { + call malloc (SMW_APLOW(new), 2*nspec, TY_REAL) + call amovr (Memr[SMW_APLOW(smw)], Memr[SMW_APLOW(new)], 2*nspec) + } + if (SMW_APHIGH(smw) != NULL) { + call malloc (SMW_APHIGH(new), 2*nspec, TY_REAL) + call amovr (Memr[SMW_APHIGH(smw)], Memr[SMW_APHIGH(new)], 2*nspec) + } + if (SMW_APIDS(smw) != NULL) { + call calloc (SMW_APIDS(new), nspec, TY_POINTER) + do i = 0, nspec-1 { + if (Memi[SMW_APIDS(smw)+i] != NULL) { + call malloc (Memi[SMW_APIDS(new)+i], SZ_LINE, TY_CHAR) + call strcpy (Memc[Memi[SMW_APIDS(smw)+i]], + Memc[Memi[SMW_APIDS(new)+i]], SZ_LINE) + } + } + } + + do i = 0, SMW_NMW(smw)-1 + SMW_MW(new,i) = mw_newcopy (SMW_MW(smw,i)) + + if (SMW_PDIM(smw) > 1) + SMW_CTLP(new) = mw_sctran (SMW_MW(new,0), "logical", "physical", 2) + + return (new) +end diff --git a/noao/onedspec/smw/smwoldms.x b/noao/onedspec/smw/smwoldms.x new file mode 100644 index 00000000..6640641f --- /dev/null +++ b/noao/onedspec/smw/smwoldms.x @@ -0,0 +1,101 @@ +include <mwset.h> +include <smw.h> + + +# SMW_OLDMS -- Convert old multispec format into MULTISPEC SMW. + +procedure smw_oldms (im, smw) + +pointer im #I IMIO pointer +pointer smw #U Input MWCS pointer, output SMW pointer + +int i, j, k, nchar, ap, beam, dtype, nw, axes[2] +double w1, dw, z +real aplow[2], aphigh[2] +pointer sp, key, val, lterm, mw, mw_open() +int imgeti(), mw_stati(), ctoi(), ctor(), ctod(), imofnlu(), imgnfn() +errchk imgstr, mw_gltermd, mw_sltermd +data axes/1,2/ + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (val, SZ_LINE, TY_CHAR) + call salloc (lterm, 12, TY_DOUBLE) + + # Set the basic multispec MWCS + i = mw_stati (smw, MW_NDIM) + j = max (2, i) + mw = mw_open (NULL, j) + call mw_newsystem (mw, "multispec", j) + call mw_swtype (mw, axes, 2, "multispec", "") + if (j > 2) + call mw_swtype (mw, 3, 1, "linear", "") + call mw_gltermd (smw, Memd[lterm+j], Memd[lterm], i) + if (i == 1) { + Memd[lterm+1] = 0. + Memd[lterm+3] = 0. + Memd[lterm+4] = 0. + Memd[lterm+5] = 1. + } + call mw_sltermd (mw, Memd[lterm+j], Memd[lterm], j) + + iferr (dtype = imgeti (im, "DC-FLAG")) + dtype = -1 + else { + call mw_swattrs (mw, 1, "label", "Wavelength") + call mw_swattrs (mw, 1, "units", "Angstroms") + } + + call mw_close (smw) + smw = mw + + # Set the SMW data structure. + call smw_open (smw, NULL, im) + do i = 1, SMW_NSPEC(smw) { + call smw_mw (smw, i, 1, mw, j, k) + call sprintf (Memc[key], SZ_FNAME, "APNUM%d") + call pargi (j) + call imgstr (im, Memc[key], Memc[val], SZ_LINE) + call imdelf (im, Memc[key]) + + k = 1 + nchar = ctoi (Memc[val], k, ap) + nchar = ctoi (Memc[val], k, beam) + nchar = ctod (Memc[val], k, w1) + nchar = ctod (Memc[val], k, dw) + nchar = ctoi (Memc[val], k, nw) + if (ctor (Memc[val], k, aplow[1]) == 0) + aplow[1] = INDEF + if (ctor (Memc[val], k, aphigh[1]) == 0) + aphigh[1] = INDEF + z = 0. + + k = dtype + if (k==1 && (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.)) + k = 0 + call smw_swattrs (smw, i, 1, ap, beam, k, w1, dw, nw, z, + aplow, aphigh, "") + + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (j) + ifnoerr (call imgstr (im, Memc[key], Memc[val], SZ_LINE)) { + call smw_sapid (smw, i, 1, Memc[val]) + call imdelf (im, Memc[key]) + } + } + + # Delete old parameters + i = imofnlu (im, + "DISPAXIS,APFORMAT,BEAM-NUM,DC-FLAG,W0,WPC,NP1,NP2") + while (imgnfn (i, Memc[key], SZ_FNAME) != EOF) { + iferr (call imdelf (im, Memc[key])) + ; + } + call imcfnl (i) + + # Update MWCS + call smw_saveim (smw, im) + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwonedspec.x b/noao/onedspec/smw/smwonedspec.x new file mode 100644 index 00000000..b7d2fa6a --- /dev/null +++ b/noao/onedspec/smw/smwonedspec.x @@ -0,0 +1,109 @@ +include <imhdr.h> +include <smw.h> + + +# SMW_ONEDSPEC -- Convert old "onedspec" format to EQUISPEC. + +procedure smw_onedspec (im, smw) + +pointer im #I IMIO pointer +pointer smw #U MWCS pointer input SMW pointer output + +int i, dtype, ap, beam, nw, imgeti(), imofnlu(), imgnfn() +real aplow[2], aphigh[2], imgetr(), mw_c1tranr() +double ltm, ltv, r, w, dw, z, imgetd() +pointer sp, key, mw, ct, mw_openim(), mw_sctran() +bool fp_equald() +errchk smw_open, smw_saxes, mw_gwtermd, mw_sctran + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + + # Convert old W0/WPC keywords if needed. + mw = smw + iferr (w = imgetd (im, "CRVAL1")) { + ifnoerr (w = imgetd (im, "W0")) { + dw = imgetd (im, "WPC") + iferr (ltm = imgetd (im, "LTM1_1")) + ltm = 1 + iferr (ltv = imgetd (im, "LTV1")) + ltv = 0 + r = ltm + ltv + dw = dw / ltm + call imaddd (im, "CRPIX1", r) + call imaddd (im, "CRVAL1", w) + call imaddd (im, "CD1_1", dw) + call imaddd (im, "CDELT1", dw) + call mw_close(mw) + mw = mw_openim (im) + } + } + + # Get dispersion and determine number of valid pixels. + call mw_gwtermd (mw, r, w, dw, 1) + w = w - (r - 1) * dw + r = 1 + call mw_swtermd (mw, r, w, dw, 1) + ct = mw_sctran (mw, "logical", "physical", 1) + nw = max (mw_c1tranr (ct, 1.), mw_c1tranr (ct, real (IM_LEN(im,1)))) + call mw_ctfree (ct) + + iferr (dtype = imgeti (im, "DC-FLAG")) { + if (fp_equald (1D0, w) || fp_equald (1D0, dw)) + dtype = DCNO + else + dtype = DCLINEAR + } + if (dtype==DCLOG) { + if (abs(w)>20. || abs(w+(nw-1)*dw)>20.) + dtype = DCLINEAR + else { + w = 10D0 ** w + dw = w * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1) + } + } + + # Convert to EQUISPEC system. + call mw_swattrs (mw, 0, "system", "equispec") + if (dtype != DCNO) { + iferr (call mw_gwattrs (mw, 1, "label", Memc[key], SZ_FNAME)) { + iferr (call mw_gwattrs (mw, 1, "units", Memc[key], SZ_FNAME)) { + call mw_swattrs (mw, 1, "units", "angstroms") + call mw_swattrs (mw, 1, "label", "Wavelength") + } + } + } + + # Set the SMW data structure. + call smw_open (smw, NULL, im) + + # Determine the aperture parameters. + iferr (beam = imgeti (im, "BEAM-NUM")) + beam = 1 + iferr (ap = imgeti (im, "APNUM")) + ap = beam + iferr (aplow[1] = imgetr (im, "APLOW")) + aplow[1] = INDEF + iferr (aphigh[1] = imgetr (im, "APHIGH")) + aphigh[1] = INDEF + iferr (z = imgetd (im, "DOPCOR")) + z = 0. + + call smw_swattrs (smw, 1, 1, ap, beam, dtype, w, dw, nw, z, + aplow, aphigh, "") + + # Delete old parameters + i = imofnlu (im, + "BEAM-NUM,APNUM,APLOW,APHIGH,DOPCOR,DC-FLAG,W0,WPC,NP1,NP2") + while (imgnfn (i, Memc[key], SZ_FNAME) != EOF) { + iferr (call imdelf (im, Memc[key])) + ; + } + call imcfnl (i) + + # Update MWCS + call smw_saveim (smw, im) + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwopen.x b/noao/onedspec/smw/smwopen.x new file mode 100644 index 00000000..782c8749 --- /dev/null +++ b/noao/onedspec/smw/smwopen.x @@ -0,0 +1,70 @@ +include <smw.h> + + +# SMW_OPEN -- Open SMW structure. +# The basic MWCS pointer and a template SMW pointer or image is input +# and the SMW pointer is returned in its place. + +procedure smw_open (mw, smw1, im) + +pointer mw #U MWCS pointer input and SMW pointer output +pointer smw1 #I Template SMW pointer +pointer im #I Template IMIO pointer + +int i, nspec, nmw, format, strdic() +pointer sp, sys, smw, mw_sctran(), mw_newcopy() +errchk smw_daxis, smw_saxes, mw_sctran + +begin + call smark (sp) + call salloc (sys, SZ_FNAME, TY_CHAR) + + call mw_gwattrs (mw, 0, "system", Memc[sys], SZ_FNAME) + format = strdic (Memc[sys], Memc[sys], SZ_FNAME, SMW_FORMATS) + + call calloc (smw, SMW_LEN(1), TY_STRUCT) + call malloc (SMW_APID(smw), SZ_LINE, TY_CHAR) + SMW_FORMAT(smw) = format + SMW_DTYPE(smw) = INDEFI + SMW_NMW(smw) = 1 + SMW_MW(smw,0) = mw + + switch (format) { + case SMW_ND: + call smw_daxis (smw, im, INDEFI, INDEFI, INDEFI) + call smw_saxes (smw, smw1, im) + + case SMW_ES: + call smw_saxes (smw, smw1, im) + + nspec = SMW_NSPEC(smw) + call calloc (SMW_APS(smw), nspec, TY_INT) + call calloc (SMW_BEAMS(smw), nspec, TY_INT) + call calloc (SMW_APLOW(smw), 2*nspec, TY_REAL) + call calloc (SMW_APHIGH(smw), 2*nspec, TY_REAL) + call calloc (SMW_APIDS(smw), nspec, TY_POINTER) + if (SMW_PDIM(smw) > 1) + SMW_CTLP(smw) = mw_sctran (mw, "logical", "physical", 2) + + case SMW_MS: + call smw_saxes (smw, smw1, im) + + nspec = SMW_NSPEC(smw) + call calloc (SMW_APIDS(smw), nspec, TY_POINTER) + if (SMW_PDIM(smw) > 1) + SMW_CTLP(smw) = mw_sctran (mw, "logical", "physical", 2) + + nmw = 1 + (nspec - 1) / SMW_NSPLIT + if (nmw > 1) { + call realloc (smw, SMW_LEN(nmw), TY_STRUCT) + call calloc (SMW_APS(smw), nspec, TY_INT) + } + do i = 1, nmw-1 + SMW_MW(smw,i) = mw_newcopy (mw) + SMW_NMW(smw) = nmw + } + + mw = smw + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwopenim.x b/noao/onedspec/smw/smwopenim.x new file mode 100644 index 00000000..468f09a7 --- /dev/null +++ b/noao/onedspec/smw/smwopenim.x @@ -0,0 +1,69 @@ +include <imhdr.h> +include <imio.h> +include <mwset.h> + +define SYSTEMS "|equispec|multispec|physical|image|world|linear|" + + +# SMW_OPENIM -- Open the spectral MWCS for various input formats. + +pointer procedure smw_openim (im) + +pointer im #I Image pointer +pointer mw #O MWCS pointer + +pointer sp, system, mw_openim() +bool streq() +int i, wcsdim, sys, strdic(), mw_stati() +errchk mw_openim, smw_oldms, smw_linear + +begin + call smark (sp) + call salloc (system, SZ_FNAME, TY_CHAR) + + # Workaround for truncation of header during image header copy. + IM_HDRLEN(im) = IM_LENHDRMEM(im) + + # Force higher dimensions to length 1. + do i = IM_NDIM(im) + 1, 3 + IM_LEN(im,i) = 1 + + mw = mw_openim (im) + call mw_seti (mw, MW_USEAXMAP, NO) + wcsdim = mw_stati (mw, MW_NDIM) + call mw_gwattrs (mw, 0, "system", Memc[system], SZ_FNAME) + sys = strdic (Memc[system], Memc[system], SZ_FNAME, SYSTEMS) + + # Set various input systems. + switch (sys) { + case 1: + call smw_equispec (im, mw) + case 2: + call smw_multispec (im, mw) + default: + if (sys == 0) { + call eprintf ( + "WARNING: Unknown coordinate system `%s' - assuming `linear'.\n") + call pargstr (Memc[system]) + } else if (sys == 3) + call mw_newsystem (mw, "image", wcsdim) + + # Old "multispec" format. + ifnoerr (call imgstr (im, "APFORMAT", Memc[system], SZ_FNAME)) { + if (streq (Memc[system], "onedspec")) + call smw_onedspec (im, mw) + else + call smw_oldms (im, mw) + + # Old "onedspec" format or other 1D image. + } else if (wcsdim == 1) { + call smw_onedspec (im, mw) + + # N-dimensional image. + } else + call smw_nd (im, mw) + } + + call sfree (sp) + return (mw) +end diff --git a/noao/onedspec/smw/smwsapid.x b/noao/onedspec/smw/smwsapid.x new file mode 100644 index 00000000..1bdf30a8 --- /dev/null +++ b/noao/onedspec/smw/smwsapid.x @@ -0,0 +1,40 @@ +include <smw.h> + + +# SMW_SAPID -- Set aperture id. + +procedure smw_sapid (smw, index1, index2, apid) + +pointer smw #I SMW pointer +int index1 #I Spectrum index +int index2 #I Spectrum index +char apid[ARB] #I Aperture id + +pointer ptr +bool streq() +errchk malloc + +begin + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call strcpy (apid, Memc[SMW_APID(smw)], SZ_LINE) + case SMW_ES, SMW_MS: + if (index1 < 0 || index1 > SMW_NSPEC(smw)) + call error (1, "smw_sapid: index out of bounds") + + if (index1 == 0) + call strcpy (apid, Memc[SMW_APID(smw)], SZ_LINE) + + else { + ptr = Memi[SMW_APIDS(smw)+index1-1] + if (streq (apid, Memc[SMW_APID(smw)])) + call mfree (ptr, TY_CHAR) + else { + if (ptr == NULL) + call malloc (ptr, SZ_LINE, TY_CHAR) + call strcpy (apid, Memc[ptr], SZ_LINE) + } + Memi[SMW_APIDS(smw)+index1-1] = ptr + } + } +end diff --git a/noao/onedspec/smw/smwsaveim.x b/noao/onedspec/smw/smwsaveim.x new file mode 100644 index 00000000..892a3319 --- /dev/null +++ b/noao/onedspec/smw/smwsaveim.x @@ -0,0 +1,251 @@ +include <imhdr.h> +include <imio.h> +include <smw.h> + + +# SMW_SAVEIM -- Save spectral WCS in image header. +# The input and output formats are EQUISPEC and MULTISPEC. A split input +# MULTISPEC WCS is first merged to a single EQUISPEC or MULTISPEC WCS. +# An input MULTISPEC WCS is converted to EQUISPEC output if possible. + +procedure smw_saveim (smw, im) + +pointer smw # SMW pointer +pointer im # Image pointer + +int i, j, format, nl, pdim, pdim1, beam, dtype, dtype1, nw, nw1 +int ap, axes[3] +real aplow[2], aphigh[2] +double v, m, w1, dw, z, w11, dw1, z1 +pointer sp, key, str1, str2, axmap, lterm, coeff, mw, mw1 + +bool strne(), fp_equald() +int imaccf(), imgeti() +pointer mw_open() +errchk smw_merge, imdelf +data axes/1,2,3/ + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (str2, SZ_LINE, TY_CHAR) + call salloc (axmap, 6, TY_INT) + call salloc (lterm, 15, TY_DOUBLE) + coeff = NULL + + # Merge split WCS into a single WCS. + call smw_merge (smw) + + mw = SMW_MW(smw,0) + pdim = SMW_PDIM(smw) + format = SMW_FORMAT(smw) + if (IM_NDIM(im) == 1) + nl = 1 + else + nl = IM_LEN(im,2) + + # If writing to an existing image we must follow IM_NPHYSDIM + # but in a NEW_COPY header we may really want a lower dimension. + # Since IM_NPHYSDIM is outside the interface we only violate + # it here and use a temporary keyword to communicate from the + # routine setting up the WCS. + + pdim1 = max (IM_NDIM(im), IM_NPHYSDIM(im)) + ifnoerr (i = imgeti (im, "SMW_NDIM")) { + pdim1 = i + call imdelf (im, "SMW_NDIM") + } + + # Check if MULTISPEC WCS can be converted to EQUISPEC. + if (format == SMW_MS) { + format = SMW_ES + do i = 1, nl { + call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + if (i == 1) { + dtype1 = dtype + w11 = w1 + dw1 = dw + z1 = z + nw1 = nw + } + if (dtype>1||dtype!=dtype1||!fp_equald(w1,w11)|| + !fp_equald(dw,dw1)||nw!=nw1||!fp_equald(z,z1)) { + format = SMW_MS + break + } + } + } + + # Save WCS in desired format. + switch (format) { + case SMW_ND: + if (SMW_DTYPE(smw) != -1) + call imaddi (im, "DC-FLAG", SMW_DTYPE(smw)) + else if (imaccf (im, "DC-FLAG") == YES) + call imdelf (im, "DC-FLAG") + if (imaccf (im, "DISPAXIS") == YES) + call imaddi (im, "DISPAXIS", SMW_PAXIS(smw,1)) + + call smw_gapid (smw, 1, 1, IM_TITLE(im), SZ_IMTITLE) + call mw_saveim (mw, im) + + case SMW_ES: + # Save aperture information. + do i = 1, nl { + call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + if (i < 1000) + call sprintf (Memc[key], SZ_FNAME, "APNUM%d") + else + call sprintf (Memc[key], SZ_FNAME, "AP%d") + call pargi (i) + call sprintf (Memc[str1], SZ_LINE, "%d %d") + call pargi (ap) + call pargi (beam) + if (!IS_INDEF(aplow[1]) || !IS_INDEF(aphigh[1])) { + call sprintf (Memc[str2], SZ_LINE, " %.2f %.2f") + call pargr (aplow[1]) + call pargr (aphigh[1]) + call strcat (Memc[str2], Memc[str1], SZ_LINE) + if (!IS_INDEF(aplow[2]) || !IS_INDEF(aphigh[2])) { + call sprintf (Memc[str2], SZ_LINE, " %.2f %.2f") + call pargr (aplow[2]) + call pargr (aphigh[2]) + call strcat (Memc[str2], Memc[str1], SZ_LINE) + } + } + call imastr (im, Memc[key], Memc[str1]) + if (i == 1) { + iferr (call imdelf (im, "APID1")) + ; + } + call smw_gapid (smw, i, 1, Memc[str1], SZ_LINE) + if (Memc[str1] != EOS) { + if (strne (Memc[str1], IM_TITLE(im))) { + if (nl == 1) { + call imastr (im, "MSTITLE", IM_TITLE(im)) + call strcpy (Memc[str1], IM_TITLE(im), SZ_IMTITLE) + } else { + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (i) + call imastr (im, Memc[key], Memc[str1]) + } + } + } + } + + # Delete unnecessary aperture information. + do i = nl+1, ARB { + if (i < 1000) + call sprintf (Memc[key], SZ_FNAME, "APNUM%d") + else + call sprintf (Memc[key], SZ_FNAME, "AP%d") + call pargi (i) + iferr (call imdelf (im, Memc[key])) + break + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (i) + iferr (call imdelf (im, Memc[key])) + ; + } + + # Add dispersion parameters to image. + if (dtype != -1) + call imaddi (im, "DC-FLAG", dtype) + else if (imaccf (im, "DC-FLAG") == YES) + call imdelf (im, "DC-FLAG") + if (nw < IM_LEN(im,1)) + call imaddi (im, "NP2", nw) + else if (imaccf (im, "NP2") == YES) + call imdelf (im, "NP2") + + # Setup EQUISPEC WCS. + + mw1 = mw_open (NULL, pdim1) + call mw_newsystem (mw1, "equispec", pdim1) + call mw_swtype (mw1, axes, pdim1, "linear", "") + ifnoerr (call mw_gwattrs (mw, 1, "label", Memc[str1], SZ_LINE)) + call mw_swattrs (mw1, 1, "label", Memc[str1]) + ifnoerr (call mw_gwattrs (mw, 1, "units", Memc[str1], SZ_LINE)) + call mw_swattrs (mw1, 1, "units", Memc[str1]) + ifnoerr (call mw_gwattrs (mw, 1, "units_display", Memc[str1], + SZ_LINE)) + call mw_swattrs (mw1, 1, "units_display", Memc[str1]) + call mw_gltermd (mw, Memd[lterm+pdim], Memd[lterm], pdim) + v = Memd[lterm] + m = Memd[lterm+pdim] + call mw_gltermd (mw1, Memd[lterm+pdim1], Memd[lterm], pdim1) + Memd[lterm] = v + Memd[lterm+pdim1] = m + call mw_sltermd (mw1, Memd[lterm+pdim1], Memd[lterm], pdim1) + call mw_gwtermd (mw1, Memd[lterm], Memd[lterm+pdim1], + Memd[lterm+2*pdim1], pdim1) + Memd[lterm] = 1. + w1 = w1 / (1 + z) + dw = dw / (1 + z) + if (dtype == DCLOG) { + dw = log10 ((w1 + (nw - 1) * dw) / w1) / (nw - 1) + w1 = log10 (w1) + } + Memd[lterm+pdim1] = w1 + Memd[lterm+2*pdim1] = dw + call mw_swtermd (mw1, Memd[lterm], Memd[lterm+pdim1], + Memd[lterm+2*pdim1], pdim1) + call mw_saveim (mw1, im) + call mw_close (mw1) + + case SMW_MS: + # Delete any APNUM keywords. If there is only one spectrum + # define the axis mapping. + + do j = 1, ARB { + if (j < 1000) + call sprintf (Memc[key], SZ_FNAME, "APNUM%d") + else + call sprintf (Memc[key], SZ_FNAME, "AP%d") + call pargi (j) + iferr (call imdelf (im, Memc[key])) + break + } + if (IM_NDIM(im) == 1) { + call aclri (Memi[axmap], 2*pdim) + Memi[axmap] = 1 + call mw_saxmap (mw, Memi[axmap], Memi[axmap+pdim], pdim) + } + + # Set aperture ids. + do i = 1, nl { + if (i == 1) { + iferr (call imdelf (im, "APID1")) + ; + } + call smw_gapid (smw, i, 1, Memc[str1], SZ_LINE) + if (Memc[str1] != EOS) { + if (strne (Memc[str1], IM_TITLE(im))) { + if (nl == 1) { + call imastr (im, "MSTITLE", IM_TITLE(im)) + call strcpy (Memc[str1], IM_TITLE(im), SZ_IMTITLE) + } else { + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (i) + call imastr (im, Memc[key], Memc[str1]) + } + } + } + } + + do i = nl+1, ARB { + call sprintf (Memc[key], SZ_FNAME, "APID%d") + call pargi (i) + iferr (call imdelf (im, Memc[key])) + break + } + + call mw_saveim (mw, im) + } + + call mfree (coeff, TY_CHAR) + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwsaxes.x b/noao/onedspec/smw/smwsaxes.x new file mode 100644 index 00000000..f5e31b63 --- /dev/null +++ b/noao/onedspec/smw/smwsaxes.x @@ -0,0 +1,247 @@ +include <imhdr.h> +include <mwset.h> +include <smw.h> + + +# SMW_SAXES -- Set axes parameters based on previously set dispersion axis. +# If the dispersion axis has been excluded for NDSPEC allow another axis to +# be chosen with a warning. For EQUISPEC and MULTISPEC require the dispersion +# to be 1 and also to be present. + +procedure smw_saxes (smw, smw1, im) + +pointer smw #I SMW pointer +pointer smw1 #I Template SMW pointer +pointer im #I Template IMIO pointer + +int i, pdim, ldim, paxis, laxis, nw, dtype, nspec +real smw_c1tranr() +double w1, dw +pointer sp, str, axno, axval, r, w, cd, mw, ct, smw_sctran() +int mw_stati(), imgeti() +bool streq(), fp_equald() +errchk smw_sctran + +begin + # If a template SMW pointer is specified just copy the axes parameters. + if (smw1 != NULL) { + call strcpy (Memc[SMW_APID(smw1)], Memc[SMW_APID(smw)], SZ_LINE) + SMW_NSPEC(smw) = SMW_NSPEC(smw1) + SMW_NBANDS(smw) = SMW_NBANDS(smw1) + SMW_TRANS(smw) = SMW_TRANS(smw1) + call amovi (SMW_PAXIS(smw1,1), SMW_PAXIS(smw,1), 3) + SMW_LDIM(smw) = SMW_LDIM(smw1) + call amovi (SMW_LAXIS(smw1,1), SMW_LAXIS(smw,1), 3) + call amovi (SMW_LLEN(smw1,1), SMW_LLEN(smw,1), 3) + call amovi (SMW_NSUM(smw1,1), SMW_NSUM(smw,1), 2) + + mw = SMW_MW(smw,0) + SMW_PDIM(smw) = mw_stati (mw, MW_NDIM) + if (SMW_PDIM(smw) > SMW_PDIM(smw1)) + do i = SMW_PDIM(smw1)+1, SMW_PDIM(smw) + SMW_PAXIS(smw,i) = i + + return + } + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (axno, 3, TY_INT) + call salloc (axval, 3, TY_INT) + call aclri (Memi[axno], 3) + + # Determine the dimensions. + mw = SMW_MW(smw,0) + pdim = mw_stati (mw, MW_NDIM) + ldim = IM_NDIM(im) + call mw_gaxmap (mw, Memi[axno], Memi[axval], pdim) + + # Set the physical dispersion axis. + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call salloc (r, pdim, TY_DOUBLE) + call salloc (w, pdim, TY_DOUBLE) + call salloc (cd, pdim*pdim, TY_DOUBLE) + + # Check for a transposed or rotated 2D image. + SMW_TRANS(smw) = NO + if (pdim == 2) { + call mw_gltermd (mw, Memd[cd], Memd[w], pdim) + if (Memd[cd] == 0D0 && Memd[cd+3] == 0D0) { + Memd[cd] = Memd[cd+1] + Memd[cd+1] = 0. + Memd[cd+3] = Memd[cd+2] + Memd[cd+2] = 0. + call mw_sltermd (mw, Memd[cd], Memd[w], pdim) + paxis = SMW_PAXIS(smw,1) + if (paxis == 1) + SMW_PAXIS(smw,1) = 2 + else + SMW_PAXIS(smw,1) = 1 + SMW_TRANS(smw) = YES + } else if (Memd[cd+1] != 0D0 || Memd[cd+2] != 0D0) { + Memd[w] = 0 + Memd[w+1] = 0 + Memd[cd] = 1 + Memd[cd+1] = 0 + Memd[cd+2] = 0 + Memd[cd+3] = 1 + call mw_sltermd (mw, Memd[cd], Memd[w], pdim) + } + } + + # If the dispersion axis is of length 1 or has been excluded find + # the first longer axis and print a warning. + + paxis = SMW_PAXIS(smw,1) + i = max (1, min (pdim, paxis)) + laxis = max (1, Memi[axno+i-1]) + if (IM_LEN(im,laxis) == 1) + do laxis = 1, ldim + if (IM_LEN(im,laxis) != 1) + break + + # Determine the number of spectra. + nspec = 1 + do i = 1, ldim + if (i != laxis) + nspec = nspec * IM_LEN(im,i) + SMW_NSPEC(smw) = nspec + SMW_NBANDS(smw) = 1 + + i = paxis + do paxis = 1, pdim + if (Memi[axno+paxis-1] == laxis) + break + + if (i != paxis && nspec > 1) { + call eprintf ( + "WARNING: Dispersion axis %d not found. Using axis %d.\n") + call pargi (i) + call pargi (paxis) + } + + # Set the dispersion system. + call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) + if (SMW_TRANS(smw) == YES) { + Memd[cd] = Memd[cd+1] + Memd[cd+1] = 0. + Memd[cd+3] = Memd[cd+2] + Memd[cd+2] = 0. + } + if (pdim == 2 && (Memd[cd+1] != 0D0 || Memd[cd+2] != 0D0)) { + iferr (dtype = imgeti (im, "DC-FLAG")) + dtype = DCNO + if (dtype != DCNO) { + call sfree (sp) + call error (1, + "Rotated, dispersion calibrated spectra are not allowed") + } + Memd[r] = 0 + Memd[r+1] = 0 + Memd[w] = 0 + Memd[w+1] = 0 + Memd[cd] = 1 + Memd[cd+1] = 0 + Memd[cd+2] = 0 + Memd[cd+3] = 1 + } + do i = 0, pdim-1 { + dw = Memd[cd+i*(pdim+1)] + if (dw == 0D0) + Memd[cd+i*(pdim+1)] = 1D0 + } + call mw_swtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) + + dw = Memd[cd+(paxis-1)*(pdim+1)] + w1 = Memd[w+paxis-1] - (Memd[r+paxis-1] - 1) * dw + nw = IM_LEN(im,laxis) + + i = 2 ** (paxis - 1) + ct = smw_sctran (smw, "logical", "physical", i) + nw = max (smw_c1tranr (ct, 0.5), smw_c1tranr (ct, nw+0.5)) + call smw_ctfree (ct) + + iferr (dtype = imgeti (im, "DC-FLAG")) { + iferr (call mw_gwattrs (mw,paxis,"axtype",Memc[str],SZ_LINE)) + Memc[str] = EOS + if (streq (Memc[str], "ra") || streq (Memc[str], "dec")) + dtype = DCNO + else if (fp_equald (1D0, w1) || fp_equald (1D0, dw)) + dtype = DCNO + else + dtype = DCLINEAR + } + if (dtype==DCLOG) { + if (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.) + dtype = DCLINEAR + else { + w1 = 10D0 ** w1 + dw = w1 * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1) + } + } + + if (dtype != DCNO) { + + + iferr (call mw_gwattrs (mw,paxis,"label",Memc[str],SZ_LINE)) { + iferr (call mw_gwattrs(mw,paxis,"units",Memc[str],SZ_LINE)) { + call mw_swattrs (mw, paxis, "units", "angstroms") + call mw_swattrs (mw, paxis, "label", "Wavelength") + } + } + } + + SMW_DTYPE(smw) = INDEFI + call smw_swattrs (smw, 1, 1, INDEFI, INDEFI, + dtype, w1, dw, nw, 0D0, INDEFR, INDEFR, "") + case SMW_ES, SMW_MS: + paxis = 1 + i = Memi[axno+1] + if (i == 0) + SMW_NSPEC(smw) = 1 + else + SMW_NSPEC(smw) = IM_LEN(im,i) + i = Memi[axno+2] + if (i == 0) + SMW_NBANDS(smw) = 1 + else + SMW_NBANDS(smw) = IM_LEN(im,i) + } + + # Check and set the physical and logical dispersion axes. + laxis = Memi[axno+paxis-1] + if (laxis == 0) { + if (Memi[axval+paxis-1] == 0) + laxis = paxis + else + call error (1, "No dispersion axis") + } + + SMW_PDIM(smw) = pdim + SMW_LDIM(smw) = ldim + SMW_PAXIS(smw,1) = paxis + SMW_LAXIS(smw,1) = laxis + SMW_LLEN(smw,1) = IM_LEN(im,laxis) + SMW_LLEN(smw,2) = 1 + SMW_LLEN(smw,3) = 1 + + # Set the spatial axes. + i = 2 + do laxis = 1, ldim { + if (laxis != SMW_LAXIS(smw,1)) { + do paxis = 1, pdim + if (Memi[axno+paxis-1] == laxis) + break + SMW_PAXIS(smw,i) = paxis + SMW_LAXIS(smw,i) = laxis + SMW_LLEN(smw,i) = IM_LEN(im,laxis) + i = i + 1 + } + } + + # Set the default title. + call smw_sapid (smw, 0, 1, IM_TITLE(im)) + + call sfree (sp) +end diff --git a/noao/onedspec/smw/smwsctran.x b/noao/onedspec/smw/smwsctran.x new file mode 100644 index 00000000..06f240db --- /dev/null +++ b/noao/onedspec/smw/smwsctran.x @@ -0,0 +1,96 @@ +include <error.h> +include <smw.h> + + +# SMW_SCTRAN -- Set the SMW coordinate system transformation. +# This routine sets up the SMW_CT structure which may consist of multiple +# pointers if the MWCS is split. If the MWCS is not split then the MWCS +# transformation routine is used directly. However if the MWCS is split then +# there may need to be an intermediate mapping from the input coordinate to +# the physical coordinate in which the split MWCS is defined as well as a +# transformation for each WCS piece. + +pointer procedure smw_sctran (smw, system1, system2, axbits) + +pointer smw #I SMW pointer +char system1[ARB] #I Input coordinate system +char system2[ARB] #I Output coordinate system +int axbits #I Bitmap defining axes to be transformed +pointer ct #O SMW CT pointer + +int i, cttype, nct, axes[3], naxes, strdic() +pointer mw_sctran() +errchk mw_sctran + +begin + # Determine the coordinate transformation type and setup the structure. + cttype = 10 * (strdic (system1, system1, ARB, SMW_CTTYPES) + 1) + + strdic (system2, system2, ARB, SMW_CTTYPES) + 1 + + nct = SMW_NMW(smw) + if (cttype == SMW_LP || cttype == SMW_PL) + nct = 1 + + call calloc (ct, SMW_CTLEN(nct), TY_STRUCT) + SMW_SMW(ct) = smw + SMW_CTTYPE(ct) = cttype + SMW_NCT(ct) = nct + + # Determine dispersion and aperture axes. + call mw_gaxlist (SMW_MW(smw,0), axbits, axes, naxes) + do i = 1, naxes { + if (axes[i] == SMW_PAXIS(smw,1)) + SMW_DAXIS(ct) = i + if (axes[i] == SMW_PAXIS(smw,2)) + SMW_AAXIS(ct) = i + } + + # If the MWCS is not split use the MWCS transformation directly. + if (nct == 1) { + iferr (SMW_CT(ct,0) = mw_sctran (SMW_MW(smw,0), system1, system2, + axbits)) { + switch (cttype) { + case SMW_WL, SMW_WP: + SMW_CT(ct,0) = mw_sctran (SMW_MW(smw,0), "physical", + system2, axbits) + case SMW_LW, SMW_PW: + SMW_CT(ct,0) = mw_sctran (SMW_MW(smw,0), system1, + "physical", axbits) + default: + call erract (EA_ERROR) + } + } + return(ct) + } + + # If there is a split MWCS then setup the intermediary transformation. + switch (cttype) { + case SMW_LW: + SMW_CTL(ct) = mw_sctran (SMW_MW(smw,0), system1, "physical", + axbits) + do i = 0, nct-1 { + iferr (SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical", + system2, axbits)) + SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical", + "physical", axbits) + } + case SMW_WL: + do i = 0, nct-1 { + iferr (SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), system1, + "physical", axbits)) + SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical", + "physical", axbits) + } + SMW_CTL(ct) = mw_sctran (SMW_MW(smw,0), "physical", system2, + axbits) + case SMW_PW, SMW_WP: + do i = 0, nct-1 { + iferr (SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), system1, + system2, axbits)) + SMW_CT(ct,i) = mw_sctran (SMW_MW(smw,i), "physical", + system2, axbits) + } + } + + return (ct) +end diff --git a/noao/onedspec/smw/smwsmw.x b/noao/onedspec/smw/smwsmw.x new file mode 100644 index 00000000..c3870c4a --- /dev/null +++ b/noao/onedspec/smw/smwsmw.x @@ -0,0 +1,21 @@ +include <smw.h> + + +# SMW_SMW -- Set MCWS pointer + +procedure smw_smw (smw, line, mw) + +pointer smw #I SMW pointer +int line #I Physical line +pointer mw #I MWCS pointer + +begin + if (SMW_NMW(smw) == 1) + SMW_MW(smw,0) = mw + + else { + if (line < 1 || line > SMW_NSPEC(smw)) + call error (1, "smw_smw: aperture not found") + SMW_MW(smw,(line-1)/SMW_NSPLIT) = mw + } +end diff --git a/noao/onedspec/smw/smwswattrs.x b/noao/onedspec/smw/smwswattrs.x new file mode 100644 index 00000000..ff859cfc --- /dev/null +++ b/noao/onedspec/smw/smwswattrs.x @@ -0,0 +1,162 @@ +include <error.h> +include <smw.h> + + +# SMW_SWATTRS -- Set spectrum attribute parameters +# This routine has the feature that if the coordinate system of a single +# spectrum in an EQUISPEC WCS is changed then the image WCS is changed +# to a MULTISPEC WCS. + +procedure smw_swattrs (smw, index1, index2, ap, beam, dtype, w1, dw, nw, z, + aplow, aphigh, coeff) + +pointer smw # SMW pointer +int index1 # Spectrum index +int index2 # Spectrum index +int ap # Aperture number +int beam # Beam number +int dtype # Dispersion type +double w1 # Starting coordinate +double dw # Coordinate interval +int nw # Number of valid pixels +double z # Redshift factor +real aplow[2], aphigh[2] # Aperture limits +char coeff[ARB] # Nonlinear coeff string + +bool fp_equald() +int i, j, sz_val, strlen() +double a, b +pointer sp, str, val, mw +errchk smw_mw + +define start_ 10 + +begin + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + +start_ + switch (SMW_FORMAT(smw)) { + case SMW_ND: + if (!IS_INDEFI(SMW_DTYPE(smw)) && (!fp_equald(w1,SMW_W1(smw)) || + !fp_equald(dw,SMW_DW(smw)) || !fp_equald(z,SMW_Z(smw)))) { + call malloc (val, 15, TY_DOUBLE) + mw = SMW_MW(smw,0) + i = SMW_PDIM(smw) + j = SMW_PAXIS(smw,1) + call mw_gwtermd (mw, Memd[val], Memd[val+i], Memd[val+2*i], i) + Memd[val+j-1] = 1. + switch (dtype) { + case DCNO, DCLINEAR: + a = w1 / (1 + z) + b = dw / (1 + z) + case DCLOG: + a = log10 (w1 / (1 + z)) + b = log10 ((w1 + (nw - 1) * dw) / w1) / (nw - 1) + case DCFUNC: + call error (1, + "Nonlinear functions not allowed for NSPEC format") + } + Memd[val+i+j-1] = a + Memd[val+2*i+(i+1)*(j-1)] = b + call mw_swtermd (mw, Memd[val], Memd[val+i], Memd[val+2*i], i) + call mfree (val, TY_DOUBLE) + } + SMW_DTYPE(smw) = dtype + SMW_NW(smw) = nw + SMW_W1(smw) = w1 + SMW_DW(smw) = dw + SMW_Z(smw) = z + + case SMW_ES: + # Check for any changes to the dispersion system. + if (dtype == DCFUNC) { + call smw_esms(smw) + goto start_ + } + if (!IS_INDEFI(SMW_DTYPE(smw)) && (dtype != SMW_DTYPE(smw) || + nw != SMW_NW(smw) || !fp_equald(w1,SMW_W1(smw)) || + !fp_equald(dw,SMW_DW(smw)) || !fp_equald(z,SMW_Z(smw)))) { + if (SMW_NSPEC(smw) > 1 && index1 > 0) { + call smw_esms(smw) + goto start_ + } + if (!fp_equald(w1,SMW_W1(smw)) || !fp_equald(dw,SMW_DW(smw)) || + !fp_equald(z,SMW_Z(smw))) { + call malloc (val, 15, TY_DOUBLE) + mw = SMW_MW(smw,0) + i = SMW_PDIM(smw) + j = SMW_PAXIS(smw,1) + call mw_gwtermd (mw, Memd[val], Memd[val+i], + Memd[val+2*i], i) + Memd[val+j-1] = 1. + switch (dtype) { + case DCNO, DCLINEAR: + a = w1 / (1 + z) + b = dw / (1 + z) + case DCLOG: + a = log10 (w1 / (1 + z)) + b = log10 ((w1 + (nw - 1) * dw) / w1) / (nw - 1) + } + Memd[val+i+j-1] = a + Memd[val+2*i+(i+1)*(j-1)] = b + call mw_swtermd (mw, Memd[val], Memd[val+i], + Memd[val+2*i], i) + call mfree (val, TY_DOUBLE) + } + } + + SMW_DTYPE(smw) = dtype + SMW_NW(smw) = nw + SMW_W1(smw) = w1 + SMW_DW(smw) = dw + SMW_Z(smw) = z + + if (index1 > 0) { + Memi[SMW_APS(smw)+index1-1] = ap + Memi[SMW_BEAMS(smw)+index1-1] = beam + Memr[SMW_APLOW(smw)+2*index1-2] = aplow[1] + Memr[SMW_APHIGH(smw)+2*index1-2] = aphigh[1] + Memr[SMW_APLOW(smw)+2*index1-1] = aplow[2] + Memr[SMW_APHIGH(smw)+2*index1-1] = aphigh[2] + } + + case SMW_MS: + # We can't use SPRINTF for the whole string because it can only + # handle a limited length and trucates long coefficient strings. + # Use STRCAT instead. + + call smw_mw (smw, index1, index2, mw, i, j) + sz_val = strlen (coeff) + SZ_LINE + call salloc (val, sz_val, TY_CHAR) + call sprintf (Memc[str], SZ_LINE, "spec%d") + call pargi (i) + call sprintf (Memc[val], sz_val, + "%d %d %d %.14g %.14g %d %.14g %.2f %.2f") + call pargi (ap) + call pargi (beam) + call pargi (dtype) + if (dtype == DCLOG) { + call pargd (log10 (w1)) + call pargd (log10 ((w1+(nw-1)*dw)/w1)/(nw-1)) + } else { + call pargd (w1) + call pargd (dw) + } + call pargi (nw) + call pargd (z) + call pargr (aplow[1]) + call pargr (aphigh[1]) + if (coeff[1] != EOS) { + call strcat (" ", Memc[val], sz_val) + call strcat (coeff, Memc[val], sz_val) + } + call mw_swattrs (mw, 2, Memc[str], Memc[val]) + + if (SMW_APS(smw) != NULL) + Memi[SMW_APS(smw)+index1-1] = ap + } + + call sfree (sp) +end diff --git a/noao/onedspec/smw/units.x b/noao/onedspec/smw/units.x new file mode 100644 index 00000000..f44abb57 --- /dev/null +++ b/noao/onedspec/smw/units.x @@ -0,0 +1,529 @@ +include <ctype.h> +include <error.h> +include <units.h> + + +# UN_OPEN -- Open units package +# It is allowed to open an unknown unit type + +pointer procedure un_open (units) + +char units[ARB] # Units string +pointer un # Units pointer returned + +begin + call calloc (un, UN_LEN, TY_STRUCT) + iferr (call un_decode (un, units)) { + call un_close (un) + call erract (EA_ERROR) + } + return (un) +end + + +# UN_CLOSE -- Close units package + +procedure un_close (un) + +pointer un # Units pointer + +begin + call mfree (un, TY_STRUCT) +end + + +# UN_COPY -- Copy units pointer + +procedure un_copy (un1, un2) + +pointer un1, un2 # Units pointers + +begin + if (un2 == NULL) + call malloc (un2, UN_LEN, TY_STRUCT) + call amovi (Memi[un1], Memi[un2], UN_LEN) +end + + +# UN_DECODE -- Decode units string and set up units structure. +# The main work is done in UN_DECODE1 so that the units string may +# be recursive; i.e. the units string may contain other units strings. +# In particular, this is required for the velocity units to specify +# a reference wavelength. + +procedure un_decode (un, units) + +pointer un # Units pointer +char units[ARB] # Units string + +bool streq() +pointer sp, units1, temp, un1, un2 +errchk un_decode1, un_ctranr + +begin + if (streq (units, UN_USER(un))) + return + + call smark (sp) + call salloc (units1, SZ_LINE, TY_CHAR) + call salloc (temp, UN_LEN, TY_STRUCT) + + # Save a copy to restore in case of an error. + call un_copy (un, temp) + + iferr { + # Decode the primary units + call un_decode1 (un, units, Memc[units1], SZ_LINE) + + # Decode velocity reference wavelength if necessary. + if (UN_CLASS(un) == UN_VEL || UN_CLASS(un) == UN_DOP) { + call salloc (un1, UN_LEN, TY_STRUCT) + call un_decode1 (un1, Memc[units1], Memc[units1], SZ_LINE) + if (UN_CLASS(un1) == UN_VEL || UN_CLASS(un1) == UN_DOP) + call error (1, + "Velocity reference units may not be velocity") + call salloc (un2, UN_LEN, TY_STRUCT) + call un_decode1 (un2, "angstroms", Memc[units1], SZ_LINE) + call un_ctranr (un1, un2, UN_VREF(un), UN_VREF(un), 1) + } + } then { + call un_copy (temp, un) + call sfree (sp) + call erract (EA_ERROR) + } + + call sfree (sp) +end + + +# UN_DECODE1 -- Decode units string and set up units structure. +# Return any secondary units string. Unknown unit strings are allowed. + +procedure un_decode1 (un, units, units1, sz_units1) + +pointer un # Units pointer +char units[ARB] # Units string +char units1[sz_units1] # Secondary units string to return +int sz_units1 # Size of secondary units string + +int unlog, uninv, untype +int i, j, k, nscan(), strdic(), strlen() +pointer sp, str +pointer stp, sym, stfind(), strefsbuf() + +int class[UN_NUNITS] +real scale[UN_NUNITS] +data stp/NULL/ +data class /UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE,UN_WAVE, + UN_FREQ,UN_FREQ,UN_FREQ,UN_FREQ,UN_VEL,UN_VEL, + UN_ENERGY,UN_ENERGY,UN_ENERGY,UN_DOP/ +data scale /UN_ANG,UN_NM,UN_MMIC,UN_MIC,UN_MM,UN_CM,UN_M,UN_HZ,UN_KHZ, + UN_MHZ,UN_GHZ,UN_MPS,UN_KPS,UN_EV,UN_KEV,UN_MEV,UN_Z/ + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + iferr (call un_abbr (stp)) + ; + + call strcpy (units, Memc[str], SZ_FNAME) + if (stp != NULL) { + sym = stfind (stp, Memc[str]) + if (sym != NULL) + call strcpy (Memc[strefsbuf(stp,Memi[sym])], + Memc[str], SZ_FNAME) + } + call strlwr (Memc[str]) + call sscan (Memc[str]) + untype = 0 + unlog = NO + uninv = NO + do i = 1, 3 { + call gargwrd (Memc[str], SZ_FNAME) + if (nscan() != i) + break + + j = strdic (Memc[str], Memc[str], SZ_FNAME, UN_DIC) + for (k=strlen(Memc[str]); k>0 && + (IS_WHITE(Memc[str+k-1]) || Memc[str+k-1]=='\n'); k=k-1) + Memc[str+k-1] = EOS + + if (j > UN_NUNITS) { + j = j - UN_NUNITS + if (j == 1) { + if (unlog == YES) + break + unlog = YES + } else if (j == 2) { + if (uninv == YES) + break + uninv = YES + } + } else { + if (class[j] == UN_VEL || class[j] == UN_DOP) { + call gargr (UN_VREF(un)) + call gargstr (units1, sz_units1) + if (nscan() != i+2) + call error (1, "Error in velocity reference wavelength") + } else + UN_VREF(un) = 0. + untype = j + break + } + } + + if (untype == 0) { + UN_TYPE(un) = 0 + UN_CLASS(un) = UN_UNKNOWN + UN_LABEL(un) = EOS + call strcpy (units, UN_UNITS(un), SZ_UNITS) + } else { + UN_TYPE(un) = untype + UN_CLASS(un) = class[untype] + UN_LOG(un) = unlog + UN_INV(un) = uninv + UN_SCALE(un) = scale[untype] + UN_LABEL(un) = EOS + UN_UNITS(un) = EOS + call strcpy (units, UN_USER(un), SZ_UNITS) + + if (unlog == YES) + call strcat ("Log ", UN_LABEL(un), SZ_UNITS) + if (uninv == YES) + call strcat ("inverse ", UN_UNITS(un), SZ_UNITS) + call strcat (Memc[str], UN_UNITS(un), SZ_UNITS) + switch (class[j]) { + case UN_WAVE: + if (uninv == NO) + call strcat ("Wavelength", UN_LABEL(un), SZ_UNITS) + else + call strcat ("Wavenumber", UN_LABEL(un), SZ_UNITS) + case UN_FREQ: + call strcat ("Frequency", UN_LABEL(un), SZ_UNITS) + case UN_VEL: + call strcat ("Velocity", UN_LABEL(un), SZ_UNITS) + case UN_ENERGY: + call strcat ("Energy", UN_LABEL(un), SZ_UNITS) + case UN_DOP: + call strcat ("Redshift", UN_LABEL(un), SZ_UNITS) + } + } + + call sfree (sp) +end + + +# UN_COMPARE -- Compare two units + +bool procedure un_compare (un1, un2) + +pointer un1, un2 # Units pointers to compare +bool strne() + +begin + if (strne (UN_UNITS(un1), UN_UNITS(un2))) + return (false) + if (strne (UN_LABEL(un1), UN_LABEL(un2))) + return (false) + if (UN_VREF(un1) != UN_VREF(un2)) + return (false) + return (true) +end + + +# UN_CTRANR -- Transform units +# Error is returned if the transform cannot be made + +procedure un_ctranr (un1, un2, val1, val2, nvals) + +pointer un1 # Input units pointer +pointer un2 # Output units pointer +real val1[nvals] # Input values +real val2[nvals] # Output values +int nvals # Number of values + +int i +real s, v, z +bool un_compare() + +begin + if (un_compare (un1, un2)) { + call amovr (val1, val2, nvals) + return + } + + if (UN_CLASS(un1) == UN_UNKNOWN || UN_CLASS(un2) == UN_UNKNOWN) + call error (1, "Cannot convert between selected units") + + call amovr (val1, val2, nvals) + + s = UN_SCALE(un1) + if (UN_LOG(un1) == YES) + do i = 1, nvals + val2[i] = 10. ** val2[i] + if (UN_INV(un1) == YES) + do i = 1, nvals + val2[i] = 1. / val2[i] + switch (UN_CLASS(un1)) { + case UN_WAVE: + do i = 1, nvals + val2[i] = val2[i] / s + case UN_FREQ: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_VEL: + v = UN_VREF(un1) + do i = 1, nvals { + z = val2[i] / s + val2[i] = sqrt ((1 + z) / (1 - z)) * v + } + case UN_ENERGY: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_DOP: + v = UN_VREF(un1) + do i = 1, nvals + val2[i] = (val2[i] / s + 1) * v + } + + s = UN_SCALE(un2) + switch (UN_CLASS(un2)) { + case UN_WAVE: + do i = 1, nvals + val2[i] = val2[i] * s + case UN_FREQ: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_VEL: + v = UN_VREF(un2) + do i = 1, nvals { + z = (val2[i] / v) ** 2 + val2[i] = (z - 1) / (z + 1) * s + } + case UN_ENERGY: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_DOP: + v = UN_VREF(un2) + do i = 1, nvals + val2[i] = (val2[i] / v - 1) * s + } + if (UN_INV(un2) == YES) + do i = 1, nvals + val2[i] = 1. / val2[i] + if (UN_LOG(un2) == YES) + do i = 1, nvals + val2[i] = log10 (val2[i]) +end + + +# UN_CHANGER -- Change units +# Error is returned if the conversion cannot be made + +procedure un_changer (un, units, vals, nvals, update) + +pointer un # Units pointer (may be changed) +char units[ARB] # Desired units +real vals[nvals] # Values +int nvals # Number of values +int update # Update units pointer? + +bool streq(), un_compare() +pointer un1, un_open() +errchk un_open, un_ctranr + +begin + + # Check for same unit string + if (streq (units, UN_USER(un))) + return + + # Check for error in units string, or the same units. + un1 = un_open (units) + if (un_compare (un1, un)) { + call strcpy (units, UN_USER(un), SZ_UNITS) + call un_close (un1) + return + } + + iferr { + call un_ctranr (un, un1, vals, vals, nvals) + if (update == YES) + call un_copy (un1, un) + call un_close(un1) + } then { + call un_close(un1) + call erract (EA_ERROR) + } +end + + +# UN_CTRAND -- Transform units +# Error is returned if the transform cannot be made + +procedure un_ctrand (un1, un2, val1, val2, nvals) + +pointer un1 # Input units pointer +pointer un2 # Output units pointer +double val1[nvals] # Input values +double val2[nvals] # Output values +int nvals # Number of values + +int i +double s, v, z +bool un_compare() + +begin + if (un_compare (un1, un2)) { + call amovd (val1, val2, nvals) + return + } + + if (UN_CLASS(un1) == UN_UNKNOWN || UN_CLASS(un2) == UN_UNKNOWN) + call error (1, "Cannot convert between selected units") + + call amovd (val1, val2, nvals) + + s = UN_SCALE(un1) + if (UN_LOG(un1) == YES) + do i = 1, nvals + val2[i] = 10. ** val2[i] + if (UN_INV(un1) == YES) + do i = 1, nvals + val2[i] = 1. / val2[i] + switch (UN_CLASS(un1)) { + case UN_WAVE: + do i = 1, nvals + val2[i] = val2[i] / s + case UN_FREQ: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_VEL: + v = UN_VREF(un1) + do i = 1, nvals { + z = val2[i] / s + val2[i] = sqrt ((1 + z) / (1 - z)) * v + } + case UN_ENERGY: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_DOP: + v = UN_VREF(un1) + do i = 1, nvals + val2[i] = (val2[i] / s + 1) * v + } + + s = UN_SCALE(un2) + switch (UN_CLASS(un2)) { + case UN_WAVE: + do i = 1, nvals + val2[i] = val2[i] * s + case UN_FREQ: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_VEL: + v = UN_VREF(un2) + do i = 1, nvals { + z = (val2[i] / v) ** 2 + val2[i] = (z - 1) / (z + 1) * s + } + case UN_ENERGY: + do i = 1, nvals + val2[i] = s / val2[i] + case UN_DOP: + v = UN_VREF(un2) + do i = 1, nvals + val2[i] = (val2[i] / v - 1) * s + } + if (UN_INV(un2) == YES) + do i = 1, nvals + val2[i] = 1. / val2[i] + if (UN_LOG(un2) == YES) + do i = 1, nvals + val2[i] = log10 (val2[i]) +end + + +# UN_CHANGED -- Change units +# Error is returned if the conversion cannot be made + +procedure un_changed (un, units, vals, nvals, update) + +pointer un # Units pointer (may be changed) +char units[ARB] # Desired units +double vals[nvals] # Values +int nvals # Number of values +int update # Update units pointer? + +bool streq(), un_compare() +pointer un1, un_open() +errchk un_open, un_ctrand + +begin + + # Check for same unit string + if (streq (units, UN_USER(un))) + return + + # Check for error in units string, or the same units. + un1 = un_open (units) + if (un_compare (un1, un)) { + call strcpy (units, UN_USER(un), SZ_UNITS) + call un_close (un1) + return + } + + iferr { + call un_ctrand (un, un1, vals, vals, nvals) + if (update == YES) + call un_copy (un1, un) + call un_close(un1) + } then { + call un_close(un1) + call erract (EA_ERROR) + } +end + + +# UN_ABBR -- Load abbreviations into a symbol table. + +procedure un_abbr (stp) + +pointer stp #U Symbol table + +int fd, open(), fscan(), nscan(), stpstr() +pointer sp, key, val +pointer sym, stopen(), stfind(), stenter(), strefsbuf() +errchk open + +begin + if (stp != NULL) + return + + fd = open (ABBREVIATIONS, READ_ONLY, TEXT_FILE) + stp = stopen ("unabbr", 20, 20, 40*SZ_LINE) + + call smark (sp) + call salloc (key, SZ_LINE, TY_CHAR) + call salloc (val, SZ_LINE, TY_CHAR) + + while (fscan (fd) != EOF) { + call gargwrd (Memc[key], SZ_LINE) + call gargwrd (Memc[val], SZ_LINE) + if (nscan() != 2) + next + if (Memc[key] == '#') + next + + sym = stfind (stp, Memc[key]) + if (sym == NULL) { + sym = stenter (stp, Memc[key], 1) + Memi[sym] = stpstr (stp, Memc[val], SZ_LINE) + } else + call strcpy (Memc[val], Memc[strefsbuf(stp,Memi[sym])], SZ_LINE) + } + + call close (fd) + call sfree (sp) +end |