From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- sys/mwcs/gen/mkpkg | 29 +++++++++++++++ sys/mwcs/gen/mwc1trand.x | 24 ++++++++++++ sys/mwcs/gen/mwc1tranr.x | 24 ++++++++++++ sys/mwcs/gen/mwc2trand.x | 38 +++++++++++++++++++ sys/mwcs/gen/mwc2tranr.x | 38 +++++++++++++++++++ sys/mwcs/gen/mwctrand.x | 97 ++++++++++++++++++++++++++++++++++++++++++++++++ sys/mwcs/gen/mwctranr.x | 97 ++++++++++++++++++++++++++++++++++++++++++++++++ sys/mwcs/gen/mwgctrand.x | 44 ++++++++++++++++++++++ sys/mwcs/gen/mwgctranr.x | 44 ++++++++++++++++++++++ sys/mwcs/gen/mwltrand.x | 26 +++++++++++++ sys/mwcs/gen/mwltranr.x | 26 +++++++++++++ sys/mwcs/gen/mwmmuld.x | 21 +++++++++++ sys/mwcs/gen/mwmmulr.x | 21 +++++++++++ sys/mwcs/gen/mwv1trand.x | 32 ++++++++++++++++ sys/mwcs/gen/mwv1tranr.x | 32 ++++++++++++++++ sys/mwcs/gen/mwv2trand.x | 49 ++++++++++++++++++++++++ sys/mwcs/gen/mwv2tranr.x | 49 ++++++++++++++++++++++++ sys/mwcs/gen/mwvmuld.x | 20 ++++++++++ sys/mwcs/gen/mwvmulr.x | 20 ++++++++++ sys/mwcs/gen/mwvtrand.x | 18 +++++++++ sys/mwcs/gen/mwvtranr.x | 18 +++++++++ 21 files changed, 767 insertions(+) create mode 100644 sys/mwcs/gen/mkpkg create mode 100644 sys/mwcs/gen/mwc1trand.x create mode 100644 sys/mwcs/gen/mwc1tranr.x create mode 100644 sys/mwcs/gen/mwc2trand.x create mode 100644 sys/mwcs/gen/mwc2tranr.x create mode 100644 sys/mwcs/gen/mwctrand.x create mode 100644 sys/mwcs/gen/mwctranr.x create mode 100644 sys/mwcs/gen/mwgctrand.x create mode 100644 sys/mwcs/gen/mwgctranr.x create mode 100644 sys/mwcs/gen/mwltrand.x create mode 100644 sys/mwcs/gen/mwltranr.x create mode 100644 sys/mwcs/gen/mwmmuld.x create mode 100644 sys/mwcs/gen/mwmmulr.x create mode 100644 sys/mwcs/gen/mwv1trand.x create mode 100644 sys/mwcs/gen/mwv1tranr.x create mode 100644 sys/mwcs/gen/mwv2trand.x create mode 100644 sys/mwcs/gen/mwv2tranr.x create mode 100644 sys/mwcs/gen/mwvmuld.x create mode 100644 sys/mwcs/gen/mwvmulr.x create mode 100644 sys/mwcs/gen/mwvtrand.x create mode 100644 sys/mwcs/gen/mwvtranr.x (limited to 'sys/mwcs/gen') diff --git a/sys/mwcs/gen/mkpkg b/sys/mwcs/gen/mkpkg new file mode 100644 index 00000000..bc8fe837 --- /dev/null +++ b/sys/mwcs/gen/mkpkg @@ -0,0 +1,29 @@ +# Make the generic portion of MWCS. + +$checkout libex.a lib$ +$udate libex.a +$checkin libex.a lib$ +$exit + +libex.a: + mwc1trand.x ../mwcs.h + mwc1tranr.x ../mwcs.h + mwc2trand.x ../mwcs.h + mwc2tranr.x ../mwcs.h + mwctrand.x ../mwcs.h + mwctranr.x ../mwcs.h + mwgctrand.x ../mwcs.h + mwgctranr.x ../mwcs.h + mwltrand.x + mwltranr.x + mwmmuld.x + mwmmulr.x + mwv1trand.x ../mwcs.h + mwv1tranr.x ../mwcs.h + mwv2trand.x ../mwcs.h + mwv2tranr.x ../mwcs.h + mwvmuld.x + mwvmulr.x + mwvtrand.x + mwvtranr.x + ; diff --git a/sys/mwcs/gen/mwc1trand.x b/sys/mwcs/gen/mwc1trand.x new file mode 100644 index 00000000..af46e02d --- /dev/null +++ b/sys/mwcs/gen/mwc1trand.x @@ -0,0 +1,24 @@ +include "../mwcs.h" + +# MW_C1TRAN -- Optimized 1D coordinate transformation. + +double procedure mw_c1trand (a_ct, x) + +pointer a_ct #I pointer to CTRAN descriptor +double x #I coordinates in input system + +double y +pointer ct + +begin + # Get real or double version of descriptor. + ct = CT_D(a_ct) + + # Perform the transformation; LNR is a simple linear transformation. + if (CT_TYPE(ct) == LNR) { + return (Memd[CT_LTM(ct)] * x + Memd[CT_LTV(ct)]) + } else { + call mw_ctrand (a_ct, x, y, 1) + return (y) + } +end diff --git a/sys/mwcs/gen/mwc1tranr.x b/sys/mwcs/gen/mwc1tranr.x new file mode 100644 index 00000000..06ad0bf7 --- /dev/null +++ b/sys/mwcs/gen/mwc1tranr.x @@ -0,0 +1,24 @@ +include "../mwcs.h" + +# MW_C1TRAN -- Optimized 1D coordinate transformation. + +real procedure mw_c1tranr (a_ct, x) + +pointer a_ct #I pointer to CTRAN descriptor +real x #I coordinates in input system + +real y +pointer ct + +begin + # Get real or double version of descriptor. + ct = CT_R(a_ct) + + # Perform the transformation; LNR is a simple linear transformation. + if (CT_TYPE(ct) == LNR) { + return (Memr[CT_LTM(ct)] * x + Memr[CT_LTV(ct)]) + } else { + call mw_ctranr (a_ct, x, y, 1) + return (y) + } +end diff --git a/sys/mwcs/gen/mwc2trand.x b/sys/mwcs/gen/mwc2trand.x new file mode 100644 index 00000000..0cb156bd --- /dev/null +++ b/sys/mwcs/gen/mwc2trand.x @@ -0,0 +1,38 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_C2TRAN -- Optimized 2D coordinate transformation. + +procedure mw_c2trand (a_ct, x1,y1, x2,y2) + +pointer a_ct #I pointer to CTRAN descriptor +double x1,y1 #I coordinates in input system +double x2,y2 #O coordinates in output system + +pointer ct, ltm, ltv +double p1[2], p2[2] + +begin + # Get real or double version of descriptor. + ct = CT_D(a_ct) + + ltm = CT_LTM(ct) + ltv = CT_LTV(ct) + + if (CT_TYPE(ct) == LNR) { + # Simple linear, nonrotated transformation. + x2 = Memd[ltm ] * x1 + Memd[ltv ] + y2 = Memd[ltm+3] * y1 + Memd[ltv+1] + } else if (CT_TYPE(ct) == LRO) { + # Linear, rotated transformation. + p1[1] = x1; p1[2] = y1 + x2 = Memd[ltm ] * p1[1] + Memd[ltm+1] * p1[2] + Memd[ltv ] + y2 = Memd[ltm+2] * p1[1] + Memd[ltm+3] * p1[2] + Memd[ltv+1] + } else { + # General case involving one or more functional terms. + p1[1] = x1; p1[2] = y1 + call mw_ctrand (a_ct, p1, p2, 2) + x2 = p2[1]; y2 = p2[2] + } +end diff --git a/sys/mwcs/gen/mwc2tranr.x b/sys/mwcs/gen/mwc2tranr.x new file mode 100644 index 00000000..ef5b5ef7 --- /dev/null +++ b/sys/mwcs/gen/mwc2tranr.x @@ -0,0 +1,38 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_C2TRAN -- Optimized 2D coordinate transformation. + +procedure mw_c2tranr (a_ct, x1,y1, x2,y2) + +pointer a_ct #I pointer to CTRAN descriptor +real x1,y1 #I coordinates in input system +real x2,y2 #O coordinates in output system + +pointer ct, ltm, ltv +real p1[2], p2[2] + +begin + # Get real or double version of descriptor. + ct = CT_R(a_ct) + + ltm = CT_LTM(ct) + ltv = CT_LTV(ct) + + if (CT_TYPE(ct) == LNR) { + # Simple linear, nonrotated transformation. + x2 = Memr[ltm ] * x1 + Memr[ltv ] + y2 = Memr[ltm+3] * y1 + Memr[ltv+1] + } else if (CT_TYPE(ct) == LRO) { + # Linear, rotated transformation. + p1[1] = x1; p1[2] = y1 + x2 = Memr[ltm ] * p1[1] + Memr[ltm+1] * p1[2] + Memr[ltv ] + y2 = Memr[ltm+2] * p1[1] + Memr[ltm+3] * p1[2] + Memr[ltv+1] + } else { + # General case involving one or more functional terms. + p1[1] = x1; p1[2] = y1 + call mw_ctranr (a_ct, p1, p2, 2) + x2 = p2[1]; y2 = p2[2] + } +end diff --git a/sys/mwcs/gen/mwctrand.x b/sys/mwcs/gen/mwctrand.x new file mode 100644 index 00000000..70e575cc --- /dev/null +++ b/sys/mwcs/gen/mwctrand.x @@ -0,0 +1,97 @@ +include "../mwcs.h" + +# MW_CTRAN -- Transform a single N-dimensional point, using the optimized +# transformation set up by a prior call to MW_SCTRAN. + +procedure mw_ctrand (a_ct, p1, p2, ndim) + +pointer a_ct #I pointer to CTRAN descriptor +double p1[ndim] #I coordinates of point in input system +double p2[ndim] #O coordinates of point in output system +int ndim #I dimensionality of point + +int naxes, i, j +pointer ct, fc, ltm, ltv, d_ct +double v1[MAX_DIM], v2[MAX_DIM], iv[MAX_DIM], ov[MAX_DIM] +errchk zcall3 + +begin + # Get real or double version of descriptor. + ct = CT_D(a_ct) + + ltm = CT_LTM(ct) + ltv = CT_LTV(ct) + + # Specially optimized cases. + if (CT_TYPE(ct) == LNR) { + # Simple linear, nonrotated transformation. + do i = 1, ndim + p2[i] = Memd[ltm+(i-1)*(ndim+1)] * p1[i] + Memd[ltv+i-1] + return + } else if (CT_TYPE(ct) == LRO) { + # Simple linear, rotated transformation. + call mw_ltrand (p1, p2, Memd[ltm], Memd[ltv], ndim) + return + } + + # If we get here the transformation involves a call to one or more + # WCS functions. In this general case, the transformation consists + # of zero or more calls to WCS functions to transform the input + # world coordinates to the linear input system, followed by a general + # linear transformation to the linear output system, followed by zero + # or more calls to WCS functions to do the forward transformation + # to generate the final output world coordinates. The WCS function + # calls are always evaluated in double precision. + + # Make zero or more WCS function calls for the different axes of the + # input system (inverse transform). + + call achtdd (p1, iv, ndim) + do j = 1, CT_NCALLI(ct) { + # Get pointer to function call descriptor. + fc = CT_FCI(ct,j) + naxes = FC_NAXES(fc) + + # Extract the coordinate vector for the function call. + do i = 1, naxes + v1[i] = p1[FC_AXIS(fc,i)] + + # Call the WCS function. + call zcall3 (FC_FCN(fc), fc, v1, v2) + + # Edit the vector IV, replacing the entries associated with + # the WCS function by the transformed values. + + do i = 1, naxes + iv[FC_AXIS(fc,i)] = v2[i] + } + + # Apply the general linear transformation. We may as well do this in + # double since we already have to use double for the function calls. + + d_ct = CT_D(a_ct) + call mw_ltrand (iv, ov, Memd[CT_LTM(d_ct)], Memd[CT_LTV(d_ct)], ndim) + + # Make zero or more WCS function calls for the different axes of the + # output system (forward transform to final world system). + + call achtdd (ov, p2, ndim) + do j = 1, CT_NCALLO(ct) { + # Get pointer to function call descriptor. + fc = CT_FCO(ct,j) + naxes = FC_NAXES(fc) + + # Extract the coordinate vector for the function call. + do i = 1, naxes + v1[i] = ov[FC_AXIS(fc,i)] + + # Call the WCS function. + call zcall3 (FC_FCN(fc), fc, v1, v2) + + # Edit the final output vector, replacing the entries for the + # function axes by their transformed values. + + do i = 1, naxes + p2[FC_AXIS(fc,i)] = v2[i] + } +end diff --git a/sys/mwcs/gen/mwctranr.x b/sys/mwcs/gen/mwctranr.x new file mode 100644 index 00000000..0574a563 --- /dev/null +++ b/sys/mwcs/gen/mwctranr.x @@ -0,0 +1,97 @@ +include "../mwcs.h" + +# MW_CTRAN -- Transform a single N-dimensional point, using the optimized +# transformation set up by a prior call to MW_SCTRAN. + +procedure mw_ctranr (a_ct, p1, p2, ndim) + +pointer a_ct #I pointer to CTRAN descriptor +real p1[ndim] #I coordinates of point in input system +real p2[ndim] #O coordinates of point in output system +int ndim #I dimensionality of point + +int naxes, i, j +pointer ct, fc, ltm, ltv, d_ct +double v1[MAX_DIM], v2[MAX_DIM], iv[MAX_DIM], ov[MAX_DIM] +errchk zcall3 + +begin + # Get real or double version of descriptor. + ct = CT_R(a_ct) + + ltm = CT_LTM(ct) + ltv = CT_LTV(ct) + + # Specially optimized cases. + if (CT_TYPE(ct) == LNR) { + # Simple linear, nonrotated transformation. + do i = 1, ndim + p2[i] = Memr[ltm+(i-1)*(ndim+1)] * p1[i] + Memr[ltv+i-1] + return + } else if (CT_TYPE(ct) == LRO) { + # Simple linear, rotated transformation. + call mw_ltranr (p1, p2, Memr[ltm], Memr[ltv], ndim) + return + } + + # If we get here the transformation involves a call to one or more + # WCS functions. In this general case, the transformation consists + # of zero or more calls to WCS functions to transform the input + # world coordinates to the linear input system, followed by a general + # linear transformation to the linear output system, followed by zero + # or more calls to WCS functions to do the forward transformation + # to generate the final output world coordinates. The WCS function + # calls are always evaluated in double precision. + + # Make zero or more WCS function calls for the different axes of the + # input system (inverse transform). + + call achtrd (p1, iv, ndim) + do j = 1, CT_NCALLI(ct) { + # Get pointer to function call descriptor. + fc = CT_FCI(ct,j) + naxes = FC_NAXES(fc) + + # Extract the coordinate vector for the function call. + do i = 1, naxes + v1[i] = p1[FC_AXIS(fc,i)] + + # Call the WCS function. + call zcall3 (FC_FCN(fc), fc, v1, v2) + + # Edit the vector IV, replacing the entries associated with + # the WCS function by the transformed values. + + do i = 1, naxes + iv[FC_AXIS(fc,i)] = v2[i] + } + + # Apply the general linear transformation. We may as well do this in + # double since we already have to use double for the function calls. + + d_ct = CT_D(a_ct) + call mw_ltrand (iv, ov, Memd[CT_LTM(d_ct)], Memd[CT_LTV(d_ct)], ndim) + + # Make zero or more WCS function calls for the different axes of the + # output system (forward transform to final world system). + + call achtdr (ov, p2, ndim) + do j = 1, CT_NCALLO(ct) { + # Get pointer to function call descriptor. + fc = CT_FCO(ct,j) + naxes = FC_NAXES(fc) + + # Extract the coordinate vector for the function call. + do i = 1, naxes + v1[i] = ov[FC_AXIS(fc,i)] + + # Call the WCS function. + call zcall3 (FC_FCN(fc), fc, v1, v2) + + # Edit the final output vector, replacing the entries for the + # function axes by their transformed values. + + do i = 1, naxes + p2[FC_AXIS(fc,i)] = v2[i] + } +end diff --git a/sys/mwcs/gen/mwgctrand.x b/sys/mwcs/gen/mwgctrand.x new file mode 100644 index 00000000..cfdca886 --- /dev/null +++ b/sys/mwcs/gen/mwgctrand.x @@ -0,0 +1,44 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_GCTRAN -- Get a coordinate transformation compiled in a previous call +# to mw_sctran. When the transformation is compiled, it is specified by +# naming the input and output systems and the axes over which the transform +# is to be performed. Rather than return this information, which the +# application already knows, we return the actual transform, i.e., the +# linear transformation matrix and translation vector comprising the linear +# portion of the transform, and axis class arrays for the input and output +# systems defining the axis types. If the axis types are all zero, there +# are no WCS function calls for any axis in either system, and the +# transformation is completely linear (hence computable by the application +# if desired, e.g., with mw_ltr). + +int procedure mw_gctrand (a_ct, o_ltm, o_ltv, axtype1, axtype2, maxdim) + +pointer a_ct #I pointer to CTRAN descriptor +double o_ltm[ARB] #O linear tranformation matrix +double o_ltv[ARB] #O translation matrix +int axtype1[ARB] #O axis types for input system +int axtype2[ARB] #O axis types for output system +int maxdim #I how much stuff to return + +pointer ct +int pdim, ndim, i, j + +begin + ct = CT_D(a_ct) + pdim = CT_NDIM(ct) + ndim = min (pdim, maxdim) + + # Output the goods. + do j = 1, ndim { + axtype1[j] = WCS_AXCLASS(CT_WCSI(ct),j) + axtype2[j] = WCS_AXCLASS(CT_WCSO(ct),j) + o_ltv[j] = Memd[CT_LTV(ct)+(j-1)] + do i = 1, ndim + o_ltm[(j-1)*ndim+i] = Memd[CT_LTM(ct)+(j-1)*pdim+(i-1)] + } + + return (pdim) +end diff --git a/sys/mwcs/gen/mwgctranr.x b/sys/mwcs/gen/mwgctranr.x new file mode 100644 index 00000000..7825c6df --- /dev/null +++ b/sys/mwcs/gen/mwgctranr.x @@ -0,0 +1,44 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_GCTRAN -- Get a coordinate transformation compiled in a previous call +# to mw_sctran. When the transformation is compiled, it is specified by +# naming the input and output systems and the axes over which the transform +# is to be performed. Rather than return this information, which the +# application already knows, we return the actual transform, i.e., the +# linear transformation matrix and translation vector comprising the linear +# portion of the transform, and axis class arrays for the input and output +# systems defining the axis types. If the axis types are all zero, there +# are no WCS function calls for any axis in either system, and the +# transformation is completely linear (hence computable by the application +# if desired, e.g., with mw_ltr). + +int procedure mw_gctranr (a_ct, o_ltm, o_ltv, axtype1, axtype2, maxdim) + +pointer a_ct #I pointer to CTRAN descriptor +real o_ltm[ARB] #O linear tranformation matrix +real o_ltv[ARB] #O translation matrix +int axtype1[ARB] #O axis types for input system +int axtype2[ARB] #O axis types for output system +int maxdim #I how much stuff to return + +pointer ct +int pdim, ndim, i, j + +begin + ct = CT_R(a_ct) + pdim = CT_NDIM(ct) + ndim = min (pdim, maxdim) + + # Output the goods. + do j = 1, ndim { + axtype1[j] = WCS_AXCLASS(CT_WCSI(ct),j) + axtype2[j] = WCS_AXCLASS(CT_WCSO(ct),j) + o_ltv[j] = Memr[CT_LTV(ct)+(j-1)] + do i = 1, ndim + o_ltm[(j-1)*ndim+i] = Memr[CT_LTM(ct)+(j-1)*pdim+(i-1)] + } + + return (pdim) +end diff --git a/sys/mwcs/gen/mwltrand.x b/sys/mwcs/gen/mwltrand.x new file mode 100644 index 00000000..d35670c7 --- /dev/null +++ b/sys/mwcs/gen/mwltrand.x @@ -0,0 +1,26 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_LTRAN -- Perform a general N-dimensional linear transformation, i.e., +# matrix multiply and translation. + +procedure mw_ltrand (p1, p2, ltm, ltv, ndim) + +double p1[ndim] #I input point +double p2[ndim] #O transformed output point +double ltm[ndim,ndim] #I linear transformation matrix +double ltv[ndim] #I linear translation vector +int ndim #I dimension of system + +int i, j +double p3[MAX_DIM] + +begin + call amovd (p1, p3, ndim) + do j = 1, ndim { + p2[j] = ltv[j] + do i = 1, ndim + p2[j] = p2[j] + ltm[i,j] * p3[i] + } +end diff --git a/sys/mwcs/gen/mwltranr.x b/sys/mwcs/gen/mwltranr.x new file mode 100644 index 00000000..9cafe4d2 --- /dev/null +++ b/sys/mwcs/gen/mwltranr.x @@ -0,0 +1,26 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_LTRAN -- Perform a general N-dimensional linear transformation, i.e., +# matrix multiply and translation. + +procedure mw_ltranr (p1, p2, ltm, ltv, ndim) + +real p1[ndim] #I input point +real p2[ndim] #O transformed output point +real ltm[ndim,ndim] #I linear transformation matrix +real ltv[ndim] #I linear translation vector +int ndim #I dimension of system + +int i, j +real p3[MAX_DIM] + +begin + call amovr (p1, p3, ndim) + do j = 1, ndim { + p2[j] = ltv[j] + do i = 1, ndim + p2[j] = p2[j] + ltm[i,j] * p3[i] + } +end diff --git a/sys/mwcs/gen/mwmmuld.x b/sys/mwcs/gen/mwmmuld.x new file mode 100644 index 00000000..ae35f082 --- /dev/null +++ b/sys/mwcs/gen/mwmmuld.x @@ -0,0 +1,21 @@ +# MW_MMUL -- Matrix multiply. + +procedure mw_mmuld (a, b, c, ndim) + +double a[ndim,ndim] #I left input matrix +double b[ndim,ndim] #I right input matrix +double c[ndim,ndim] #O output matrix +int ndim #I dimensionality of system + +int i, j, k +double v + +begin + do j = 1, ndim + do i = 1, ndim { + v = 0 + do k = 1, ndim + v = v + a[k,j] * b[i,k] + c[i,j] = v + } +end diff --git a/sys/mwcs/gen/mwmmulr.x b/sys/mwcs/gen/mwmmulr.x new file mode 100644 index 00000000..83e14d2c --- /dev/null +++ b/sys/mwcs/gen/mwmmulr.x @@ -0,0 +1,21 @@ +# MW_MMUL -- Matrix multiply. + +procedure mw_mmulr (a, b, c, ndim) + +real a[ndim,ndim] #I left input matrix +real b[ndim,ndim] #I right input matrix +real c[ndim,ndim] #O output matrix +int ndim #I dimensionality of system + +int i, j, k +real v + +begin + do j = 1, ndim + do i = 1, ndim { + v = 0 + do k = 1, ndim + v = v + a[k,j] * b[i,k] + c[i,j] = v + } +end diff --git a/sys/mwcs/gen/mwv1trand.x b/sys/mwcs/gen/mwv1trand.x new file mode 100644 index 00000000..3c3ac124 --- /dev/null +++ b/sys/mwcs/gen/mwv1trand.x @@ -0,0 +1,32 @@ +include "../mwcs.h" + +# MW_V1TRAN -- Optimized 1D coordinate transformation for an array of points. + +procedure mw_v1trand (a_ct, x1, x2, npts) + +pointer a_ct #I pointer to CTRAN descriptor +double x1[ARB] #I coordinates in input system +double x2[ARB] #O coordinates in output system +int npts + +int i +pointer ct +double scale, offset +errchk mw_ctrand + +begin + # Get real or double version of descriptor. + ct = CT_D(a_ct) + + scale = Memd[CT_LTM(ct)] + offset = Memd[CT_LTV(ct)] + + # Perform the transformation; case LNR is a simple linear transform. + if (CT_TYPE(ct) == LNR) { + do i = 1, npts + x2[i] = scale * x1[i] + offset + } else { + do i = 1, npts + call mw_ctrand (a_ct, x1[i], x2[i], 1) + } +end diff --git a/sys/mwcs/gen/mwv1tranr.x b/sys/mwcs/gen/mwv1tranr.x new file mode 100644 index 00000000..045f6a33 --- /dev/null +++ b/sys/mwcs/gen/mwv1tranr.x @@ -0,0 +1,32 @@ +include "../mwcs.h" + +# MW_V1TRAN -- Optimized 1D coordinate transformation for an array of points. + +procedure mw_v1tranr (a_ct, x1, x2, npts) + +pointer a_ct #I pointer to CTRAN descriptor +real x1[ARB] #I coordinates in input system +real x2[ARB] #O coordinates in output system +int npts + +int i +pointer ct +real scale, offset +errchk mw_ctranr + +begin + # Get real or double version of descriptor. + ct = CT_R(a_ct) + + scale = Memr[CT_LTM(ct)] + offset = Memr[CT_LTV(ct)] + + # Perform the transformation; case LNR is a simple linear transform. + if (CT_TYPE(ct) == LNR) { + do i = 1, npts + x2[i] = scale * x1[i] + offset + } else { + do i = 1, npts + call mw_ctranr (a_ct, x1[i], x2[i], 1) + } +end diff --git a/sys/mwcs/gen/mwv2trand.x b/sys/mwcs/gen/mwv2trand.x new file mode 100644 index 00000000..3a1cf329 --- /dev/null +++ b/sys/mwcs/gen/mwv2trand.x @@ -0,0 +1,49 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_V2TRAN -- Optimized 2D coordinate transformation for an array of points. + +procedure mw_v2trand (a_ct, x1,y1, x2,y2, npts) + +pointer a_ct #I pointer to CTRAN descriptor +double x1[ARB],y1[ARB] #I coordinates in input system +double x2[ARB],y2[ARB] #O coordinates in output system +int npts + +int i +pointer ct, ltm, ltv +double p1[2], p2[2] +errchk mw_ctrand + +begin + # Get real or double version of descriptor. + ct = CT_D(a_ct) + + ltm = CT_LTM(ct) + ltv = CT_LTV(ct) + + if (CT_TYPE(ct) == LNR) { + # Simple linear, nonrotated transformation. + do i = 1, npts { + x2[i] = Memd[ltm ] * x1[i] + Memd[ltv ] + y2[i] = Memd[ltm+3] * y1[i] + Memd[ltv+1] + } + } else if (CT_TYPE(ct) == LRO) { + # Linear, rotated transformation. + do i = 1, npts { + p1[1] = x1[i]; p1[2] = y1[i] + x2[i] = Memd[ltm ] * p1[1] + Memd[ltm+1] * p2[1] + + Memd[ltv ] + y2[i] = Memd[ltm+2] * p1[1] + Memd[ltm+3] * p2[1] + + Memd[ltv+1] + } + } else { + # General case involving one or more functional terms. + do i = 1, npts { + p1[1] = x1[i]; p1[2] = y1[i] + call mw_ctrand (a_ct, p1, p2, 2) + x2[i] = p2[1]; y2[i] = p2[2] + } + } +end diff --git a/sys/mwcs/gen/mwv2tranr.x b/sys/mwcs/gen/mwv2tranr.x new file mode 100644 index 00000000..dc2fe58f --- /dev/null +++ b/sys/mwcs/gen/mwv2tranr.x @@ -0,0 +1,49 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../mwcs.h" + +# MW_V2TRAN -- Optimized 2D coordinate transformation for an array of points. + +procedure mw_v2tranr (a_ct, x1,y1, x2,y2, npts) + +pointer a_ct #I pointer to CTRAN descriptor +real x1[ARB],y1[ARB] #I coordinates in input system +real x2[ARB],y2[ARB] #O coordinates in output system +int npts + +int i +pointer ct, ltm, ltv +real p1[2], p2[2] +errchk mw_ctranr + +begin + # Get real or double version of descriptor. + ct = CT_R(a_ct) + + ltm = CT_LTM(ct) + ltv = CT_LTV(ct) + + if (CT_TYPE(ct) == LNR) { + # Simple linear, nonrotated transformation. + do i = 1, npts { + x2[i] = Memr[ltm ] * x1[i] + Memr[ltv ] + y2[i] = Memr[ltm+3] * y1[i] + Memr[ltv+1] + } + } else if (CT_TYPE(ct) == LRO) { + # Linear, rotated transformation. + do i = 1, npts { + p1[1] = x1[i]; p1[2] = y1[i] + x2[i] = Memr[ltm ] * p1[1] + Memr[ltm+1] * p2[1] + + Memr[ltv ] + y2[i] = Memr[ltm+2] * p1[1] + Memr[ltm+3] * p2[1] + + Memr[ltv+1] + } + } else { + # General case involving one or more functional terms. + do i = 1, npts { + p1[1] = x1[i]; p1[2] = y1[i] + call mw_ctranr (a_ct, p1, p2, 2) + x2[i] = p2[1]; y2[i] = p2[2] + } + } +end diff --git a/sys/mwcs/gen/mwvmuld.x b/sys/mwcs/gen/mwvmuld.x new file mode 100644 index 00000000..0af8dfa7 --- /dev/null +++ b/sys/mwcs/gen/mwvmuld.x @@ -0,0 +1,20 @@ +# MW_VMUL -- Vector multiply. + +procedure mw_vmuld (a, b, c, ndim) + +double a[ndim,ndim] #I input matrix +double b[ndim] #I input vector +double c[ndim] #O output vector +int ndim #I system dimension + +int i, j +double v + +begin + do j = 1, ndim { + v = 0 + do i = 1, ndim + v = v + a[i,j] * b[i] + c[j] = v + } +end diff --git a/sys/mwcs/gen/mwvmulr.x b/sys/mwcs/gen/mwvmulr.x new file mode 100644 index 00000000..54a0776e --- /dev/null +++ b/sys/mwcs/gen/mwvmulr.x @@ -0,0 +1,20 @@ +# MW_VMUL -- Vector multiply. + +procedure mw_vmulr (a, b, c, ndim) + +real a[ndim,ndim] #I input matrix +real b[ndim] #I input vector +real c[ndim] #O output vector +int ndim #I system dimension + +int i, j +real v + +begin + do j = 1, ndim { + v = 0 + do i = 1, ndim + v = v + a[i,j] * b[i] + c[j] = v + } +end diff --git a/sys/mwcs/gen/mwvtrand.x b/sys/mwcs/gen/mwvtrand.x new file mode 100644 index 00000000..1a1cb662 --- /dev/null +++ b/sys/mwcs/gen/mwvtrand.x @@ -0,0 +1,18 @@ +# MW_VTRAN -- Transform an array of N-dimensional points, expressed as a +# 2D vector where v[1,i] is point I of vector V. + +procedure mw_vtrand (ct, v1, v2, ndim, npts) + +pointer ct #I pointer to CTRAN descriptor +double v1[ndim,npts] #I points to be transformed +double v2[ndim,npts] #O vector to get the transformed points +int ndim #I dimensionality of each point +int npts #I number of points + +int i +errchk mw_ctrand + +begin + do i = 1, npts + call mw_ctrand (ct, v1[1,i], v2[1,i], ndim) +end diff --git a/sys/mwcs/gen/mwvtranr.x b/sys/mwcs/gen/mwvtranr.x new file mode 100644 index 00000000..ca705c8b --- /dev/null +++ b/sys/mwcs/gen/mwvtranr.x @@ -0,0 +1,18 @@ +# MW_VTRAN -- Transform an array of N-dimensional points, expressed as a +# 2D vector where v[1,i] is point I of vector V. + +procedure mw_vtranr (ct, v1, v2, ndim, npts) + +pointer ct #I pointer to CTRAN descriptor +real v1[ndim,npts] #I points to be transformed +real v2[ndim,npts] #O vector to get the transformed points +int ndim #I dimensionality of each point +int npts #I number of points + +int i +errchk mw_ctranr + +begin + do i = 1, npts + call mw_ctranr (ct, v1[1,i], v2[1,i], ndim) +end -- cgit