aboutsummaryrefslogtreecommitdiff
path: root/noao/twodspec/apextract/apmw.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/twodspec/apextract/apmw.x')
-rw-r--r--noao/twodspec/apextract/apmw.x280
1 files changed, 280 insertions, 0 deletions
diff --git a/noao/twodspec/apextract/apmw.x b/noao/twodspec/apextract/apmw.x
new file mode 100644
index 00000000..9fcd35d7
--- /dev/null
+++ b/noao/twodspec/apextract/apmw.x
@@ -0,0 +1,280 @@
+include <error.h>
+include <imhdr.h>
+include <imio.h>
+include <mwset.h>
+
+
+# APMW_OPEN -- Open APMW structure.
+# APMW_CLOSE -- Close APMW structure.
+# APMW_SETAP -- Set aperture values in APMW structure.
+# APMW_SAVEIM -- Set WCS in image header.
+# APMW_WCSFIX -- Fix up WCS
+
+# Output formats
+define ONEDSPEC 1 # Individual 1D spectra
+define MULTISPEC 2 # Multiple spectra
+define ECHELLE 3 # Echelle spectra
+define STRIP 4 # Strip spectra
+define NORM 5 # Normalized spectra
+define FLAT 6 # Flat spectra
+define RATIO 7 # Ratio of data to model
+define DIFF 8 # Difference of data and model
+define FIT 9 # Model
+define NOISE 10 # Noise calculation
+
+
+# Data structure for the apertures. This version assumes the coordinates
+# are the same for all the apertures.
+
+define APMW_LEN (8 + $1 * 6) # Structure length
+
+define APMW_LABEL Memi[$1] # WCS label
+define APMW_UNITS Memi[$1+1] # WCS units
+define APMW_DTYPE Memi[$1+2] # Dispersion type
+define APMW_NW Memi[$1+3] # Number of pixels
+define APMW_W1 Memd[P2D($1+4)] # Starting coordinate
+define APMW_DW Memd[P2D($1+6)] # Coordinate per pixel
+define APMW_AP Memi[$1+6*($2-1)+8] # Aperture
+define APMW_BEAM Memi[$1+6*($2-1)+9] # Beam
+define APMW_APLOW Memd[P2D($1+6*($2-1)+10)] # Aperture low
+define APMW_APHIGH Memd[P2D($1+6*($2-1)+12)] # Aperture high
+
+
+# APMW_OPEN -- Open APMW structure.
+
+pointer procedure apmw_open (in, out, dispaxis, naps, nw)
+
+pointer in #I Input IMIO pointer
+pointer out #I Output IMIO pointer
+int dispaxis #I Input dispersion axis
+int naps #I Number of apertures
+int nw #I Number of dispersion pixels
+pointer apmw #O Returned APMW pointer
+
+int imgeti()
+double mw_c1trand()
+pointer mw, ct, mw_openim(), mw_sctran()
+errchk mw_openim, mw_sctran, mw_c1trand, apmw_wcsfix
+
+begin
+ # Allocate data structure.
+ call malloc (apmw, APMW_LEN(naps), TY_STRUCT)
+ call malloc (APMW_LABEL(apmw), SZ_LINE, TY_CHAR)
+ call malloc (APMW_UNITS(apmw), SZ_LINE, TY_CHAR)
+
+ # Set defaults.
+ call strcpy ("Pixel", Memc[APMW_LABEL(apmw)], SZ_LINE)
+ Memc[APMW_UNITS(apmw)] = EOS
+ APMW_DTYPE(apmw,i) = -1
+ APMW_NW(apmw,i) = nw
+ APMW_W1(apmw,i) = 1.
+ APMW_DW(apmw,i) = 1.
+
+ # Get WCS info from input image.
+ iferr {
+ mw = mw_openim (in)
+ iferr (APMW_DTYPE(apmw) = imgeti (in, "DC-FLAG"))
+ APMW_DTYPE(apmw) = -1
+ iferr (call mw_gwattrs (mw, dispaxis, "label",
+ Memc[APMW_LABEL(apmw)], SZ_LINE)) {
+ if (APMW_DTYPE(apmw) == -1)
+ call strcpy ("Pixel", Memc[APMW_LABEL(apmw)], SZ_LINE)
+ else
+ call strcpy ("Wavelength", Memc[APMW_LABEL(apmw)], SZ_LINE)
+ }
+ iferr (call mw_gwattrs (mw, dispaxis, "units",
+ Memc[APMW_UNITS(apmw)], SZ_LINE)) {
+ if (APMW_DTYPE(apmw) == -1)
+ Memc[APMW_UNITS(apmw)] = EOS
+ else
+ call strcpy ("Angstroms", Memc[APMW_UNITS(apmw)], SZ_LINE)
+ }
+
+ call apmw_wcsfix (in, mw)
+ iferr (ct = mw_sctran (mw, "logical", "world", dispaxis))
+ call error (1,
+ "Coordinate system ignored (rotated?). Using pixel coordinates.")
+ APMW_W1(apmw) = mw_c1trand (ct, 1D0)
+ APMW_DW(apmw) = mw_c1trand (ct, double (nw))
+ APMW_DW(apmw) = (APMW_DW(apmw)-APMW_W1(apmw))/(nw-1)
+ } then
+ call erract (EA_WARN)
+
+ call mw_close (mw)
+
+ return (apmw)
+end
+
+
+# APMW_CLOSE -- Close APMW structure.
+
+procedure apmw_close (apmw)
+
+pointer apmw # APMW pointer
+
+begin
+ call mfree (APMW_LABEL(apmw), TY_CHAR)
+ call mfree (APMW_UNITS(apmw), TY_CHAR)
+ call mfree (apmw, TY_STRUCT)
+end
+
+
+# APMW_SETAP -- Set aperture values in APMW structure.
+
+procedure apmw_setap (apmw, line, ap, beam, aplow, aphigh)
+
+pointer apmw # APMW pointer
+int line # Image line
+int ap # Aperture
+int beam # Beam
+real aplow # Aperture lower limit
+real aphigh # Aperture upper limit
+
+begin
+ APMW_AP(apmw,line) = ap
+ APMW_BEAM(apmw,line) = beam
+ APMW_APLOW(apmw,line) = aplow
+ APMW_APHIGH(apmw,line) = aphigh
+end
+
+
+# APMW_SAVEIM -- Save WCS in image header.
+
+procedure apmw_saveim (apmw, im, fmt)
+
+pointer apmw #I APMW pointer
+pointer im #I IMIO pointer
+int fmt #I Output format
+
+int i, naps, wcsdim, axes[3], imaccf()
+double r[3], w[3], cd[9]
+bool strne()
+pointer sp, key, str, mw, list, mw_open(), imofnlu(), imgnfn()
+errchk imdelf
+data axes/1,2,3/
+
+begin
+ call smark (sp)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ if (fmt == STRIP)
+ naps = 1
+ else
+ naps = IM_LEN(im, 2)
+
+ # Workaround for truncation of header during image header copy.
+ IM_HDRLEN(im) = IM_LENHDRMEM(im)
+
+ # Delete keywords.
+ list = imofnlu (im, "SLFIB[0-9]*")
+ while (imgnfn (list, Memc[key], SZ_FNAME) != EOF)
+ call imdelf (im, Memc[key])
+ call imcfnl (list)
+
+ # Add aperture parameters to image header.
+ do i = 1, naps {
+ call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
+ call pargi (i)
+ call sprintf (Memc[str], SZ_LINE, "%d %d %.2f %.2f")
+ call pargi (APMW_AP(apmw,i))
+ call pargi (APMW_BEAM(apmw,i))
+ call pargd (APMW_APLOW(apmw,i))
+ call pargd (APMW_APHIGH(apmw,i))
+ call imastr (im, Memc[key], Memc[str])
+ if (naps == 1) {
+ call sprintf (Memc[key], SZ_FNAME, "APID%d")
+ call pargi (i)
+ ifnoerr (call imgstr (im, Memc[key], Memc[str], SZ_LINE)) {
+ if (strne (Memc[str], IM_TITLE(im))) {
+ call imastr (im, "MSTITLE", IM_TITLE(im))
+ call strcpy (Memc[str], IM_TITLE(im), SZ_IMTITLE)
+ }
+ call imdelf (im, Memc[key])
+ }
+ }
+ }
+
+ # Add dispersion parameters to image header.
+ if (APMW_DTYPE(apmw) != -1)
+ call imaddi (im, "DC-FLAG", APMW_DTYPE(apmw))
+ else if (imaccf (im, "DC-FLAG") == YES)
+ call imdelf (im, "DC-FLAG")
+ if (APMW_NW(apmw) < IM_LEN(im,1))
+ call imaddi (im, "NP2", APMW_NW(apmw))
+ else if (imaccf (im, "NP2") == YES)
+ call imdelf (im, "NP2")
+ iferr (call imdelf (im, "dispaxis"))
+ ;
+ if (fmt == STRIP)
+ call imaddi (im, "dispaxis", 1)
+
+ # Set WCS in image header.
+ wcsdim = IM_NPHYSDIM(im)
+ mw = mw_open (NULL, wcsdim)
+ if (fmt == STRIP)
+ call mw_newsystem (mw, "linear", wcsdim)
+ else
+ call mw_newsystem (mw, "equispec", wcsdim)
+ call mw_swtype (mw, axes, wcsdim, "linear", "")
+ if (Memc[APMW_LABEL(apmw)] != EOS)
+ call mw_swattrs (mw, 1, "label", Memc[APMW_LABEL(apmw)])
+ if (Memc[APMW_UNITS(apmw)] != EOS)
+ call mw_swattrs (mw, 1, "units", Memc[APMW_UNITS(apmw)])
+
+ call aclrd (r, 3)
+ call aclrd (w, 3)
+ call aclrd (cd, 9)
+ r[1] = 1.
+ w[1] = APMW_W1(apmw)
+ cd[1] = APMW_DW(apmw)
+ if (wcsdim == 2)
+ cd[4] = 1.
+ if (wcsdim == 3) {
+ cd[5] = 1.
+ cd[9] = 1.
+ }
+ call mw_swtermd (mw, r, w, cd, wcsdim)
+
+ call mw_saveim (mw, im)
+ call mw_close (mw)
+
+ call sfree (sp)
+end
+
+
+# APMW_WCSFIX -- Fix up WCS to avoid CDELT=0 which occurs if there are WCS
+# keywords in the header but no CDELT.
+
+procedure apmw_wcsfix (im, mw)
+
+pointer im # IMIO pointer
+pointer mw # MWCS pointer
+
+int i, ndim, mw_stati()
+double val
+pointer sp, r, w, cd
+errchk mw_gwtermd, mw_swtermd
+
+begin
+ call mw_seti (mw, MW_USEAXMAP, NO)
+ ndim = mw_stati (mw, MW_NDIM)
+
+ call smark (sp)
+ call salloc (r, ndim, TY_DOUBLE)
+ call salloc (w, ndim, TY_DOUBLE)
+ call salloc (cd, ndim*ndim, TY_DOUBLE)
+
+ # Check cd terms. Assume no rotation.
+ call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], ndim)
+ do i = 0, ndim-1 {
+ val = Memd[cd+i*(ndim+1)]
+ if (val == 0D0) {
+ Memd[w+i] = 1D0
+ Memd[cd+i*(ndim+1)] = 1D0
+ }
+ }
+ call mw_swtermd (mw, Memd[r], Memd[w], Memd[cd], ndim)
+
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 11:07:44 -0400