Initial commit

1 files changed, 222 insertions, 0 deletions

diff --git a/noao/onedspec/t_specshift.x b/noao/onedspec/t_specshift.x
new file mode 100644
index 00000000..e5b8ea0e
--- /dev/null
+++ b/noao/onedspec/t_specshift.x
@@ -0,0 +1,222 @@
+include	<error.h>
+include	<smw.h>
+
+# Function types.
+define  CHEBYSHEV       1       # CURFIT Chebyshev polynomial
+define  LEGENDRE        2       # CURFIT Legendre polynomial
+define  SPLINE3         3       # CURFIT cubic spline
+define  SPLINE1         4       # CURFIT linear spline
+define  PIXEL           5       # pixel coordinate array
+define  SAMPLE          6       # sampled coordinates
+
+
+# T_SSHIFT -- Shift the spectral coordinates
+ 
+procedure t_sshift ()
+ 
+int	list			# Input list of spectra
+double	shift			# Shift to apply
+pointer	aps			# Aperture list
+bool	verbose			# Verbose?
+ 
+int	ap, beam, dtype, nw
+double	w1, dw, z
+real	aplow[2], aphigh[2]
+pointer	sp, image, coeff, tmp, im, mw
+
+bool	clgetb()
+double	clgetd()
+int	imtopenp(), imtgetim()
+pointer	rng_open(), immap(), smw_openim()
+errchk	immap, smw_openim, smw_gwattrs, smw_swattrs, sshift
+ 
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	coeff = NULL
+
+	list = imtopenp ("spectra")
+	shift = clgetd ("shift")
+	call clgstr ("apertures", Memc[image], SZ_FNAME)
+	verbose = clgetb ("verbose")
+
+	iferr (aps = rng_open (Memc[image], INDEF, INDEF, INDEF))
+	    call error (0, "Bad aperture list")
+
+	while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+	    im = NULL
+	    mw = NULL
+	    iferr {
+	        tmp = immap (Memc[image], READ_WRITE, 0); im = tmp
+		tmp = smw_openim (im); mw = tmp
+
+		switch (SMW_FORMAT(mw)) {
+		case SMW_ND:
+		    call smw_gwattrs (mw, 1, 1, ap, beam, dtype,
+			w1, dw, nw, z, aplow, aphigh, coeff)
+		    w1 = w1 + shift
+		    call smw_swattrs (mw, 1, 1, ap, beam, dtype,
+			w1, dw, nw, z, aplow, aphigh, Memc[coeff])
+		    if (verbose) {
+			call printf ("%s: shift = %g, %g --> %g\n")
+			    call pargstr (Memc[image])
+			    call pargd (shift)
+			    call pargd (w1 - shift)
+			    call pargd (w1)
+		    }
+		case SMW_ES, SMW_MS:
+		    call sshift (im, mw, Memc[image], aps, shift,
+			verbose)
+		}
+	    } then
+		call erract (EA_WARN)
+
+	    if (mw != NULL) {
+		call smw_saveim (mw, im)
+		call smw_close (mw)
+	    }
+	    if (im != NULL)
+		call imunmap (im)
+	}
+
+	call rng_close (aps)
+	call imtclose (list)
+	call mfree (coeff, TY_CHAR)
+	call sfree (sp)
+end
+
+
+# SSHIFT -- Shift coordinate zero point of selected aperture in
+# MULTISPEC and EQUISPEC images.
+
+procedure sshift (im, mw, image, aps, shift, verbose)
+
+pointer	im			# IMIO pointer
+pointer	mw			# MWCS pointer
+char	image[ARB]		# Image name
+pointer	aps			# Aperture range list
+double	shift			# Shift to add
+bool	verbose			# Verbose?
+
+int	i, ap, beam, dtype, nw, naps
+double	w1, dw, z
+real	aplow[2], aphigh[2]
+pointer	coeff, coeffs
+bool	rng_elementi()
+errchk	sshift1
+
+begin
+	coeff = NULL
+	coeffs = NULL
+ 
+	# Go through each spectrum and change the selected apertures.
+	naps = 0
+	do i = 1, SMW_NSPEC(mw) {
+	    # Get aperture info
+	    iferr (call smw_gwattrs (mw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, coeff))
+		break
+
+	    # Check if aperture is to be changed
+	    if (!rng_elementi (aps, ap))
+		next
+
+	    # Apply shift
+	    w1 = w1 + shift
+	    if (dtype == 2)
+		call sshift1 (shift, coeff)
+
+	    call smw_swattrs (mw, i, 1, ap, beam, dtype, w1, dw, nw, z,
+		aplow, aphigh, Memc[coeff])
+
+	    # Make record
+	    if (verbose) {
+		if (naps == 1) {
+		    call printf ("%s: shift = %g\n")
+			call pargstr (image)
+			call pargd (shift)
+		}
+		call printf ("  Aperture %d:  %g --> %g\n")
+		    call pargi (ap)
+		    call pargd (w1 - shift)
+		    call pargd (w1)
+	    }
+	}
+
+	call mfree (coeff, TY_CHAR)
+	call mfree (coeffs, TY_DOUBLE)
+end
+
+
+# SSHIFT1 -- Shift coordinate zero point of nonlinear functions.
+
+procedure sshift1 (shift, coeff)
+
+double	shift			# Shift to add
+pointer	coeff			# Attribute function coefficients
+
+int	i, j, ip, nalloc, ncoeff, type, order, fd
+double	dval
+pointer	coeffs
+int	ctod(), stropen()
+errchk	stropen
+
+begin
+	if (coeff == NULL)
+	    return
+	if (Memc[coeff] == EOS)
+	    return
+
+	coeffs = NULL
+	ncoeff = 0
+	ip = 1
+	while (ctod (Memc[coeff], ip, dval) > 0) {
+	    if (coeffs == NULL) {
+		nalloc = 10
+		call malloc (coeffs, nalloc, TY_DOUBLE)
+	    } else if (ncoeff == nalloc) {
+		nalloc = nalloc + 10
+		call realloc (coeffs, nalloc, TY_DOUBLE)
+	    }
+	    Memd[coeffs+ncoeff] = dval
+	    ncoeff = ncoeff + 1
+	}
+	ip = ip + SZ_LINE
+	call realloc (coeff, ip, TY_CHAR)
+	call aclrc (Memc[coeff], ip)
+	fd = stropen (Memc[coeff], ip, NEW_FILE)
+
+	ip = 0
+	while (ip < ncoeff) {
+	    if (ip > 0)
+		call fprintf (fd, " ")
+	    Memd[coeffs+ip+1] = Memd[coeffs+ip+1] + shift
+	    type = nint (Memd[coeffs+ip+2])
+	    order = nint (Memd[coeffs+ip+3])
+	    call fprintf (fd, "%.3g %g %d %d")
+		call pargd (Memd[coeffs+ip])
+		call pargd (Memd[coeffs+ip+1])
+		call pargi (type)
+		call pargi (order)
+	    switch (type) {
+	    case CHEBYSHEV, LEGENDRE:
+		j = 6 + order
+	    case SPLINE3:
+		j = 9 + order
+	    case SPLINE1:
+		j = 7 + order
+	    case PIXEL:
+		j = 4 + order
+	    case SAMPLE:
+		j = 5 + order
+	    }
+	    do i = 4, j-1 {
+		call fprintf (fd, " %g")
+		    call pargd (Memd[coeffs+ip+i])
+	    }
+	    ip = ip + j
+	}
+	call strclose (fd)
+
+	call mfree (coeffs, TY_DOUBLE)
+end