diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2016-10-27 18:59:29 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2016-10-27 18:59:29 (GMT) |
commit | d4d595fa7fb12903db9227d33d48b2b00120dbd1 (patch) | |
tree | 7d18365de0d6d1b29399b6a17c7eb01c2eb3ed49 /tksao/wcssubs/slasubs.c | |
parent | 949f96e29bfe0bd8710d775ce220e597064e2589 (diff) | |
download | blt-d4d595fa7fb12903db9227d33d48b2b00120dbd1.zip blt-d4d595fa7fb12903db9227d33d48b2b00120dbd1.tar.gz blt-d4d595fa7fb12903db9227d33d48b2b00120dbd1.tar.bz2 |
Initial commit
Diffstat (limited to 'tksao/wcssubs/slasubs.c')
-rw-r--r-- | tksao/wcssubs/slasubs.c | 364 |
1 files changed, 364 insertions, 0 deletions
diff --git a/tksao/wcssubs/slasubs.c b/tksao/wcssubs/slasubs.c new file mode 100644 index 0000000..74ddb88 --- /dev/null +++ b/tksao/wcssubs/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 + */ |