Merge commit '23c7930d27fe11c4655e1291a07a821dbbaba78d' as 'funtools'

1 files changed, 364 insertions, 0 deletions

diff --git a/funtools/wcs/slasubs.c b/funtools/wcs/slasubs.c
new file mode 100644
index 0000000..74ddb88
--- /dev/null
+++ b/funtools/wcs/slasubs.c
@@ -0,0 +1,364 @@
+/* File slasubs.c
+ *** Starlink subroutines by Patrick Wallace used by wcscon.c subroutines
+ *** April 13, 1998
+ */
+
+#include <math.h>
+#include <string.h>
+
+/*  slaDcs2c (a, b, v): Spherical coordinates to direction cosines.
+ *  slaDcc2s (v, a, b):  Direction cosines to spherical coordinates.
+ *  slaDmxv (dm, va, vb): vector vb = matrix dm * vector va
+ *  slaImxv (rm, va, vb): vector vb = (inverse of matrix rm) * vector va
+ *  slaDranrm (angle):  Normalize angle into range 0-2 pi.
+ *  slaDrange (angle):  Normalize angle into range +/- pi.
+ *  slaDeuler (order, phi, theta, psi, rmat)
+ *	      Form a rotation matrix from the Euler angles - three successive
+ *	      rotations about specified Cartesian axes.
+ */
+
+void
+slaDcs2c (a, b, v)
+
+double a;	/* Right ascension in radians */
+double b;	/* Declination in radians */
+double *v;	/* x,y,z unit vector (returned) */
+
+/*
+**  slaDcs2c: Spherical coordinates to direction cosines.
+**
+**  The spherical coordinates are longitude (+ve anticlockwise
+**  looking from the +ve latitude pole) and latitude.  The
+**  Cartesian coordinates are right handed, with the x axis
+**  at zero longitude and latitude, and the z axis at the
+**  +ve latitude pole.
+**
+**  P.T.Wallace   Starlink   31 October 1993
+*/
+{
+    double cosb;
+ 
+    cosb = cos ( b );
+    v[0] = cos ( a ) * cosb;
+    v[1] = sin ( a ) * cosb;
+    v[2] = sin ( b );
+}
+
+
+void
+slaDcc2s (v, a, b)
+
+double *v;	/* x,y,z vector */
+double *a;	/* Right ascension in radians */
+double *b;	/* Declination in radians */
+
+/*
+**  slaDcc2s:
+**  Direction cosines to spherical coordinates.
+**
+**  Returned:
+**     *a,*b  double      spherical coordinates in radians
+**
+**  The spherical coordinates are longitude (+ve anticlockwise
+**  looking from the +ve latitude pole) and latitude.  The
+**  Cartesian coordinates are right handed, with the x axis
+**  at zero longitude and latitude, and the z axis at the
+**  +ve latitude pole.
+**
+**  If v is null, zero a and b are returned.
+**  At either pole, zero a is returned.
+**
+**  P.T.Wallace   Starlink   31 October 1993
+*/
+{
+    double x, y, z, r;
+ 
+    x = v[0];
+    y = v[1];
+    z = v[2];
+    r = sqrt ( x * x + y * y );
+ 
+    *a = ( r != 0.0 ) ? atan2 ( y, x ) : 0.0;
+    *b = ( z != 0.0 ) ? atan2 ( z, r ) : 0.0;
+}
+
+
+void
+slaDmxv (dm, va, vb)
+
+double (*dm)[3];	/* 3x3 Matrix */
+double *va;		/* Vector */
+double *vb;		/* Result vector (returned) */
+
+/*
+**  slaDmxv:
+**  Performs the 3-d forward unitary transformation:
+**     vector vb = matrix dm * vector va
+**
+**  P.T.Wallace   Starlink   31 October 1993
+*/
+{
+    int i, j;
+    double w, vw[3];
+ 
+    /* Matrix dm * vector va -> vector vw */
+    for ( j = 0; j < 3; j++ ) {
+	w = 0.0;
+	for ( i = 0; i < 3; i++ ) {
+	    w += dm[j][i] * va[i];
+	    }
+	vw[j] = w;
+	}
+ 
+    /* Vector vw -> vector vb */
+    for ( j = 0; j < 3; j++ ) {
+	vb[j] = vw[j];
+	}
+}
+
+
+void slaDimxv (dm, va, vb)
+     double (*dm)[3];
+     double *va;
+     double *vb;
+/*
+**  - - - - - - - - -
+**   s l a D i m x v
+**  - - - - - - - - -
+**
+**  Performs the 3-d backward unitary transformation:
+**
+**     vector vb = (inverse of matrix dm) * vector va
+**
+**  (double precision)
+**
+**  (n.b.  The matrix must be unitary, as this routine assumes that
+**   the inverse and transpose are identical)
+**
+**
+**  Given:
+**     dm       double[3][3]   matrix
+**     va       double[3]      vector
+**
+**  Returned:
+**     vb       double[3]      result vector
+**
+**  P.T.Wallace   Starlink   31 October 1993
+*/
+{
+  long i, j;
+  double w, vw[3];
+ 
+/* Inverse of matrix dm * vector va -> vector vw */
+   for ( j = 0; j < 3; j++ ) {
+      w = 0.0;
+      for ( i = 0; i < 3; i++ ) {
+         w += dm[i][j] * va[i];
+      }
+      vw[j] = w;
+   }
+ 
+/* Vector vw -> vector vb */
+   for ( j = 0; j < 3; j++ ) {
+     vb[j] = vw[j];
+   }
+}
+
+ 
+/* 2pi */
+#define D2PI 6.2831853071795864769252867665590057683943387987502
+
+/* pi */
+#define DPI 3.1415926535897932384626433832795028841971693993751
+
+double slaDranrm (angle)
+
+double angle;	/* angle in radians */
+
+/*
+**  slaDranrm:
+**  Normalize angle into range 0-2 pi.
+**  The result is angle expressed in the range 0-2 pi (double).
+**  Defined in slamac.h:  D2PI
+**
+**  P.T.Wallace   Starlink   30 October 1993
+*/
+{
+    double w;
+ 
+    w = fmod ( angle, D2PI );
+    return ( w >= 0.0 ) ? w : w + D2PI;
+}
+
+#ifndef dsign
+#define dsign(A,B) ((B)<0.0?-(A):(A))
+#endif
+
+double
+slaDrange (angle)
+     double angle;
+/*
+**  - - - - - - - - - -
+**   s l a D r a n g e
+**  - - - - - - - - - -
+**
+**  Normalize angle into range +/- pi.
+**
+**  (double precision)
+**
+**  Given:
+**     angle     double      the angle in radians
+**
+**  The result is angle expressed in the +/- pi (double precision).
+**
+**  Defined in slamac.h:  DPI, D2PI
+**
+**  P.T.Wallace   Starlink   31 October 1993
+*/
+{
+  double w;
+ 
+  w = fmod ( angle, D2PI );
+  return ( fabs ( w ) < DPI ) ? w : w - dsign ( D2PI, angle );
+}
+
+
+void
+slaDeuler (order, phi, theta, psi, rmat)
+
+char *order;		/* specifies about which axes the rotations occur */
+double phi;		/* 1st rotation (radians) */
+double theta;		/* 2nd rotation (radians) */
+double psi;		/* 3rd rotation (radians) */
+double (*rmat)[3];	/* 3x3 Rotation matrix (returned) */
+
+/*
+**  slaDeuler:
+**  Form a rotation matrix from the Euler angles - three successive
+**  rotations about specified Cartesian axes.
+**
+**  A rotation is positive when the reference frame rotates
+**  anticlockwise as seen looking towards the origin from the
+**  positive region of the specified axis.
+**
+**  The characters of order define which axes the three successive
+**  rotations are about.  A typical value is 'zxz', indicating that
+**  rmat is to become the direction cosine matrix corresponding to
+**  rotations of the reference frame through phi radians about the
+**  old z-axis, followed by theta radians about the resulting x-axis,
+**  then psi radians about the resulting z-axis.
+**
+**  The axis names can be any of the following, in any order or
+**  combination:  x, y, z, uppercase or lowercase, 1, 2, 3.  Normal
+**  axis labelling/numbering conventions apply;  the xyz (=123)
+**  triad is right-handed.  Thus, the 'zxz' example given above
+**  could be written 'zxz' or '313' (or even 'zxz' or '3xz').  Order
+**  is terminated by length or by the first unrecognised character.
+**
+**  Fewer than three rotations are acceptable, in which case the later
+**  angle arguments are ignored.  Zero rotations produces a unit rmat.
+**
+**  P.T.Wallace   Starlink   17 November 1993
+*/
+{
+   int j, i, l, n, k;
+   double result[3][3], rotn[3][3], angle, s, c , w, wm[3][3];
+   char axis;
+ 
+/* Initialize result matrix */
+   for ( j = 0; j < 3; j++ ) {
+      for ( i = 0; i < 3; i++ ) {
+         result[i][j] = ( i == j ) ? 1.0 : 0.0;
+      }
+   }
+ 
+/* Establish length of axis string */
+   l = strlen ( order );
+ 
+/* Look at each character of axis string until finished */
+   for ( n = 0; n < 3; n++ ) {
+      if ( n <= l ) {
+ 
+      /* Initialize rotation matrix for the current rotation */
+         for ( j = 0; j < 3; j++ ) {
+            for ( i = 0; i < 3; i++ ) {
+               rotn[i][j] = ( i == j ) ? 1.0 : 0.0;
+            }
+         }
+ 
+      /* Pick up the appropriate Euler angle and take sine & cosine */
+         switch ( n ) {
+         case 0 :
+           angle = phi;
+           break;
+         case 1 :
+           angle = theta;
+           break;
+         case 2 :
+           angle = psi;
+           break;
+         }
+         s = sin ( angle );
+         c = cos ( angle );
+ 
+      /* Identify the axis */
+         axis =  order[n];
+         if ( ( axis == 'X' ) || ( axis == 'x' ) || ( axis == '1' ) ) {
+ 
+         /* Matrix for x-rotation */
+            rotn[1][1] = c;
+            rotn[1][2] = s;
+            rotn[2][1] = -s;
+            rotn[2][2] = c;
+         }
+         else if ( ( axis == 'Y' ) || ( axis == 'y' ) || ( axis == '2' ) ) {
+ 
+         /* Matrix for y-rotation */
+            rotn[0][0] = c;
+            rotn[0][2] = -s;
+            rotn[2][0] = s;
+            rotn[2][2] = c;
+         }
+         else if ( ( axis == 'Z' ) || ( axis == 'z' ) || ( axis == '3' ) ) {
+ 
+         /* Matrix for z-rotation */
+            rotn[0][0] = c;
+            rotn[0][1] = s;
+            rotn[1][0] = -s;
+            rotn[1][1] = c;
+         } else {
+ 
+         /* Unrecognized character - fake end of string */
+            l = 0;
+         }
+ 
+      /* Apply the current rotation (matrix rotn x matrix result) */
+         for ( i = 0; i < 3; i++ ) {
+            for ( j = 0; j < 3; j++ ) {
+               w = 0.0;
+               for ( k = 0; k < 3; k++ ) {
+                  w += rotn[i][k] * result[k][j];
+               }
+               wm[i][j] = w;
+            }
+         }
+         for ( j = 0; j < 3; j++ ) {
+            for ( i= 0; i < 3; i++ ) {
+               result[i][j] = wm[i][j];
+            }
+         }
+      }
+   }
+ 
+/* Copy the result */
+   for ( j = 0; j < 3; j++ ) {
+      for ( i = 0; i < 3; i++ ) {
+         rmat[i][j] = result[i][j];
+      }
+   }
+}
+/*
+ * Nov  4 1996	New file
+ *
+ * Apr 13 1998	Add list of subroutines to start of file
+ */