Initial commit

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