1 files changed, 117 insertions, 0 deletions
diff --git a/sys/mwcs/mwtransd.x b/sys/mwcs/mwtransd.x
new file mode 100644
index 00000000..267317ce
--- /dev/null
+++ b/sys/mwcs/mwtransd.x
@@ -0,0 +1,117 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<syserr.h>
+include	"mwcs.h"
+
+# MW_TRANSLATE -- Translate the logical system, i.e., perform a linear
+# transformation of the logical system by modifying the Lterm of the MWCS.
+# The transformation is defined in terms of the CURRENT LOGICAL SYSTEM,
+# subject to axis mapping, dimensional reduction, etc.  This is unlike
+# MW_SLTERM, which defines the Lterm relative to the physical system in
+# physical terms (no axis mapping, full dimensionality, etc.).
+#
+#	p' = ltm * (p - ltv_1) + ltv_2
+#	
+# For convenience the transformation is specified using separate translation
+# vectors for the input and output systems.  If ltv_1 is set to zero a
+# "fully reduced" transformation of the form used internally may be entered.
+
+procedure mw_translated (mw, ltv_1, ltm, ltv_2, ndim)
+
+pointer	mw			#I pointer to MWCS descriptor
+double	ltv_1[ndim]		#I input translation vector
+double	ltm[ndim,ndim]		#I linear transformation matrix
+double	ltv_2[ndim]		#I output translation vector
+int	ndim			#I dimensionality of transform
+
+double	v
+pointer	sp, o_ltm, o_ltv, n_ltm, n_ltv, ltv
+int	pdim, nelem, axis[MAX_DIM], i, j
+errchk	syserrs
+define	err_ 91
+
+begin
+	pdim = MI_NDIM(mw)
+	nelem = pdim * pdim
+
+	call smark (sp)
+	call salloc (ltv, ndim, TY_DOUBLE)
+	call salloc (o_ltm, nelem, TY_DOUBLE)
+	call salloc (o_ltv, pdim, TY_DOUBLE)
+	call salloc (n_ltm, nelem, TY_DOUBLE)
+	call salloc (n_ltv, pdim, TY_DOUBLE)
+
+	# Combine the input and output translation vectors.
+	do j = 1, ndim {
+	    v = ltv_2[j]
+	    do i = 1, ndim
+		v = v + ltm[i,j] * (-ltv_1[i])
+	    Memd[ltv+j-1] = v
+	}
+
+	# Get axis map.
+	if (MI_USEAXMAP(mw) == NO) {
+	    if (ndim > MI_NDIM(mw))
+		goto err_
+	    do i = 1, ndim
+		axis[i] = i
+	} else {
+	    if (ndim > MI_NLOGDIM(mw))
+err_		call syserrs (SYS_MWNDIM, "mw_translate")
+	    do i = 1, ndim
+		axis[i] = MI_PHYSAX(mw,i)
+	}
+
+	# Perform the transformation.  Use a procedure call to dereference
+	# the pointers to simplify the notation.
+
+	call mw_axtran (D(mw,MI_LTM(mw)), D(mw,MI_LTV(mw)),
+	    Memd[n_ltm], Memd[n_ltv], pdim, ltm, Memd[ltv], axis, ndim)
+
+	# Update the Lterm.
+	call amovd (Memd[n_ltm], D(mw,MI_LTM(mw)), nelem)
+	call amovd (Memd[n_ltv], D(mw,MI_LTV(mw)), pdim)
+
+	call sfree (sp)
+end
+
+
+# MW_AXTRAN -- Axis mapped linear transformation.  Matrix or vector elements
+# not included in the axis map are propagated unchanged.
+
+procedure mw_axtran (o_ltm,o_ltv, n_ltm,n_ltv, pdim, ltm,ltv, ax, ndim)
+
+double	o_ltm[pdim,pdim]	#I matrix to be transformed
+double	o_ltv[pdim]		#I vector to be transformed
+double	n_ltm[pdim,pdim]	#O transformed matrix
+double	n_ltv[pdim]		#O transformed vector
+int	pdim			#I dimension of these guys
+double	ltm[ndim,ndim]		#I transform matrix
+double	ltv[ndim]		#I transform vector
+int	ax[ndim]		#I transform axis map: physax=axis[logax]
+int	ndim			#I dimension of these guys
+
+double	v
+int	i, j, k
+
+begin
+	# Transform the matrix.
+	call amovd (o_ltm, n_ltm, pdim * pdim)
+	do j = 1, ndim
+	    do i = 1, ndim {
+		v = 0
+		do k = 1, ndim
+		    # v = v + o_ltm[ax[k],ax[j]] * ltm[i,k]
+		    v = v + ltm[k,j] * o_ltm[ax[i],ax[k]]
+		n_ltm[ax[i],ax[j]] = v
+	    }
+
+	# Transform the vector.	    
+	call amovd (o_ltv, n_ltv, pdim)
+	do j = 1, ndim {
+	    v = ltv[j]
+	    do i = 1, ndim
+		v = v + ltm[i,j] * o_ltv[ax[i]]
+	    n_ltv[ax[j]] = v
+	}
+end