wcscon.c

Convert between various sky coordinate systems
Based on Starlink subroutines by Patrick Wallace

Subroutines

  • fk524 (ra,dec) Convert J2000(FK5) to B1950(FK4) coordinates
  • fk524e (ra, dec, epoch) (more accurate for known position epoch
  • fk524m (ra,dec,rapm,decpm) exact
  • fk425 (ra,dec) Convert B1950(FK4) to J2000(FK5) coordinates
  • fk425e (ra, dec, epoch) (more accurate for known position epoch
  • fk425m (ra, dec, rapm, decpm) exact
  • fk42gal (dtheta,dphi) Convert B1950(FK4) to galactic coordinates
  • fk52gal (dtheta,dphi) Convert J2000(FK5) to galactic coordinates
  • gal2fk4 (dtheta,dphi) Convert galactic coordinates to B1950(FK4)
  • gal2fk5 (dtheta,dphi) Convert galactic coordinates to J2000(FK5)
  • fk5prec (ep0, ep1, ra, dec) Precession ep0 to ep1, FK5 system
  • fk4prec (ep0, ep1, ra, dec) Precession ep0 to ep1, FK4 system

    Convert J2000 FK5 star data to B1950 FK4

    Based on Starlink sla_fk524 by P.T.Wallace

    fk524 (double *ra, double *dec)

    ra
    Right ascension in degrees (J2000 in, B1950 out)
    dec
    Declination in degrees (J2000 in, B1950 out)

    fk524e (double *ra, double *dec, double epoch)

    ra
    Right ascension in degrees (J2000 in, B1950 out)
    dec
    Declination in degrees (J2000 in, B1950 out)
    epoch
    Besselian epoch in years

    fk524m (double *ra, double *dec, double *rapm, double *decpm)

    ra
    Right ascension in radians (J2000 in, B1950 out)
    dec
    Declination in radians (J2000 in, B1950 out)
    rapm
    Proper motion in right ascension (rad/jul.yr. in, rad/trop.yr. out)
    decpm
    Proper motion in declination (rad/jul.yr. in, rad/trop.yr. out)
    These routines convert stars from the new, IAU 1976, FK5, Fricke system, to the old, Bessel-Newcomb, FK4 system, using Yallop's implementation (see ref 2) of a matrix method due to Standish (see ref 3). The numerical values of ref 2 are used canonically.

    Notes

  • The proper motions in ra are dra / dt rather than cos(dec) * dra / dt, and are per year rather than per century.

  • Note that conversion from Julian epoch 2000.0 to Besselian epoch 1950.0 only is provided for. Conversions involving other epochs will require use of the appropriate precession, proper motion, and e-terms routines before and/or after fk524 is called.

  • In the FK4 catalogue the proper motions of stars within 10 degrees of the poles do not embody the differential e-term effect and should, strictly speaking, be handled in a different manner from stars outside these regions. However, given the general lack of homogeneity of the star data available for routine astrometry, the difficulties of handling positions that may have been determined from astrometric fields spanning the polar and non-polar regions, the likelihood that the differential e-terms effect was not taken into account when allowing for proper motion in past astrometry, and the undesirability of a discontinuity in the algorithm, the decision has been made in this routine to include the effect of differential e-terms on the proper motions for all stars, whether polar or not. At epoch 2000, and measuring on the sky rather than in terms of dra, the errors resulting from this simplification are less than 1 milliarcsecond in position and 1 milliarcsecond per century in proper motion.

    References

    1. Mean and apparent place computations in the new IAU System. I. The transformation of astrometric catalog systems to the equinox J2000.0. Smith, C.A.; Kaplan, G.H.; Hughes, J.A.; Seidelmann, P.K.; Yallop, B.D.; Hohenkerk, C.Y. Astronomical Journal vol. 97, Jan. 1989, p. 265-273.
    2. Mean and apparent place computations in the new IAU System. II. Transformation of mean star places from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space. Yallop, B.D.; Hohenkerk, C.Y.; Smith, C.A.; Kaplan, G.H.; Hughes, J.A.; Seidelmann, P.K.; Astronomical Journal vol. 97, Jan. 1989, p. 274-279.
    3. Conversion of positions and proper motions from B1950.0 to the IAU system at J2000.0, Standish, E.M. Astronomy and Astrophysics, vol. 115, no. 1, Nov. 1982, p. 20-22.

    Convert B1950.0 FK4 star data to J2000.0 FK5

    fk425 (double *ra, double *dec)

    ra
    Right ascension in degrees (B1950 in, J2000 out)
    dec
    Declination in degrees (B1950 in, J2000 out)

    fk425e (double *ra, double *dec, double epoch)

    ra
    Right ascension in degrees (B1950 in, J2000 out)
    dec
    Declination in degrees (B1950 in, J2000 out)
    epoch
    Besselian epoch in years

    fk425m (double *ra, double *dec, double *rapm, double *decpm)

    ra
    Right ascension in radians (B1950 in, J2000 out)
    dec
    Declination in radians (B1950 in, J2000 out)
    rapm
    Proper motion in right ascension (rad/trop.yr. in, rad/jul.yr. out)
    decpm
    Proper motion in declination (rad/trop.yr. in, rad/jul.yr. out)
    These routines convert stars from the old, Bessel-Newcomb, FK4 system to the new, IAU 1976, FK5, Fricke system, using Yallop's implementation (see ref 2) of a matrix method due to Standish (see ref 3). The numerical values of ref 2 are used canonically.

    Notes

    1) The proper motions in ra are dra/dt rather than cos(dec)*dra/dt, and are per year rather than per century.

    2) Conversion from besselian epoch 1950.0 to Julian epoch 2000.0 only is provided for. Conversions involving other epochs will require use of the appropriate precession, proper motion, and e-terms routines before and/or after fk425 is called.
    /* l2,b2 system of galactic coordinates * p = 192.25 ra of galactic north pole (mean b1950.0) * q = 62.6 inclination of galactic to mean b1950.0 equator * r = 33 longitude of ascending node * p,q,r are degrees * Equatorial to galactic rotation matrix (The Eulerian angles are p, q, 90-r) +cp.cq.sr-sp.cr +sp.cq.sr+cp.cr -sq.sr -cp.cq.cr-sp.sr -sp.cq.cr+cp.sr +sq.cr cp.sq +sp.sq +cq */ static double bgal[3][3] = { -0.066988739415,-0.872755765852,-0.483538914632, 0.492728466075,-0.450346958020, 0.744584633283, -0.867600811151,-0.188374601723, 0.460199784784}; /*--- Transform b1950.0 'fk4' equatorial coordinates to * IAU 1958 galactic coordinates */ void fk42gal (dtheta,dphi) double *dtheta; /* b1950.0 'fk4' ra in radians Galactic longitude (l2) in radians (returned) */ double *dphi; /* b1950.0 'fk4' dec in radians Galactic latitude (b2) in radians (returned) */ /* Note: The equatorial coordinates are b1950.0 'fk4'. use the routine jpgalj if conversion from j2000.0 coordinates is required. Reference: blaauw et al, MNRAS,121,123 (1960) */ { double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec; double rcon = 0.0174532925199433; void jpcon(),jpcop(); int i; char *eqcoor, *eqstrn(); dra = *dtheta; ddec = *dphi; rra = dra * rcon; rdec = ddec * rcon; /* remove e-terms */ /* call jpabe (rra,rdec,-1,idg) */ /* Spherical to Cartesian */ r = 1.; jpcop (rra,rdec,r,pos); /* rotate to galactic */ for (i = 0; i<3; i++) { pos1[i] = pos[0]*bgal[i][0] + pos[1]*bgal[i][1] + pos[2]*bgal[i][2]; } /* Cartesian to spherical */ jpcon (pos1,&rl,&rb,&r); dl = rl / rcon; db = rb / rcon; *dtheta = dl; *dphi = db; /* Print result if in diagnostic mode */ if (idg) { eqcoor = eqstrn (dra,ddec); printf ("FK42GAL: B1950 RA,Dec= %s\n",eqcoor); printf ("FK42GAL: long = %.5f lat = %.5f\n",dl,db); free (eqcoor); } return; } /*--- Transform IAU 1958 galactic coordinates to B1950.0 'fk4' * equatorial coordinates */ void gal2fk4 (dtheta,dphi) double *dtheta; /* Galactic longitude (l2) in radians B1950 FK4 RA in radians (returned) */ double *dphi; /* Galactic latitude (b2) in radians B1950 FK4 Dec in radians (returned) */ /* Note: The equatorial coordinates are B1950.0 FK4. Use the routine GAL2FK5 if conversion to J2000 coordinates is required. Reference: Blaauw et al, MNRAS,121,123 (1960) */ { double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec; double rcon = 0.0174532925199433; void jpcon(),jpcop(); char *eqcoor, *eqstrn(); int i; /* spherical to cartesian */ dl = *dtheta; db = *dphi; rl = dl * rcon; rb = db * rcon; r = 1.0; jpcop (rl,rb,r,pos); /* rotate to equatorial coordinates */ for (i = 0; i < 3; i++) { pos1[i] = pos[0]*bgal[0][i] + pos[1]*bgal[1][i] + pos[2]*bgal[2][i]; } /* cartesian to spherical */ jpcon (pos1,&rra,&rdec,&r); /* introduce e-terms */ /* jpabe (rra,rdec,-1,idg); */ dra = rra / rcon; ddec = rdec / rcon; *dtheta = dra; *dphi = ddec; /* print result if in diagnostic mode */ if (idg) { printf ("GAL2FK4: long = %.5f lat = %.5f\n",dl,db); eqcoor = eqstrn (dra,ddec); printf ("GAL2FK4: B1950 RA,Dec= %s\n",eqcoor); free (eqcoor); } return; } /* l2,b2 system of galactic coordinates p = 192.25 ra of galactic north pole (mean b1950.0) q = 62.6 inclination of galactic to mean b1950.0 equator r = 33 longitude of ascending node p,q,r are degrees */ /* Equatorial to galactic rotation matrix The eulerian angles are p, q, 90-r +cp.cq.sr-sp.cr +sp.cq.sr+cp.cr -sq.sr -cp.cq.cr-sp.sr -sp.cq.cr+cp.sr +sq.cr +cp.sq +sp.sq +cq */ static double jgal[3][3] = { -0.054875539726,-0.873437108010,-0.483834985808, 0.494109453312,-0.444829589425, 0.746982251810, -0.867666135858,-0.198076386122, 0.455983795705}; /* Transform J2000 equatorial coordinates to IAU 1958 galactic coordinates */ void fk52gal (dtheta,dphi) double *dtheta; /* J2000 right ascension in degrees Galactic longitude (l2) in degrees (returned) */ double *dphi; /* J2000 declination in degrees Galactic latitude (b2) in degrees (returned) */ /* Rotation matrices by P.T.Wallace, Starlink eqgal and galeq, March 1986 */ /* Note: The equatorial coordinates are J2000 FK5. Use the routine GAL2FK4 if conversion from B1950 FK4 coordinates is required. Reference: Blaauw et al, MNRAS,121,123 (1960) */ { double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec; double rcon = 0.0174532925199433; void jpcon(),jpcop(); char *eqcoor, *eqstrn(); int i; /* Spherical to cartesian */ dra = *dtheta; ddec = *dphi; rra = dra * rcon; rdec = ddec * rcon; r = 1.0; (void)jpcop (rra,rdec,r,pos); /* Rotate to galactic */ for (i = 0; i < 3; i++) { pos1[i] = pos[0]*jgal[i][0] + pos[1]*jgal[i][1] + pos[2]*jgal[i][2]; } /* Cartesian to spherical */ jpcon (pos1,&rl,&rb,&r); dl = rl / rcon; db = rb / rcon; *dtheta = dl; *dphi = db; /* Print result if in diagnostic mode */ if (idg) { eqcoor = eqstrn (dra,ddec); printf ("FK52GAL: J2000 RA,Dec= %s\n",eqcoor); printf ("FK52GAL: long = %.5f lat = %.5f\n",dl,db); free (eqcoor); } return; } /*--- Transform IAU 1958 galactic coordinates to J2000 equatorial coordinates */ void gal2fk5 (dtheta,dphi) double *dtheta; /* Galactic longitude (l2) in degrees J2000.0 ra in degrees (returned) */ double *dphi; /* Galactic latitude (b2) in degrees J2000.0 dec in degrees (returned) */ /* Note: The equatorial coordinates are J2000. Use the routine FK42GAL if conversion to J2000 coordinates is required. Reference: Blaauw et al, MNRAS,121,123 (1960) */ { double pos[3],pos1[3],r,dl,db,rl,rb,rra,rdec,dra,ddec; double rcon = 0.0174532925199433; void jpcon(),jpcop(); int i; char *eqcoor, *eqstrn(); /* Spherical to Cartesian */ dl = *dtheta; db = *dphi; rl = dl * rcon; rb = db * rcon; r = 1.0; jpcop (rl,rb,r,pos); /* Rotate to equatorial coordinates */ for (i = 0; i < 3; i++) { pos1[i] = pos[0]*jgal[0][i] + pos[1]*jgal[1][i] + pos[2]*jgal[2][i]; } /* Cartesian to Spherical */ jpcon (pos1,&rra,&rdec,&r); dra = rra / rcon; ddec = rdec / rcon; *dtheta = dra; *dphi = ddec; /* Print result if in diagnostic mode */ if (idg) { printf ("GAL2FK5: long = %.5f lat = %.5f\n",dl,db); eqcoor = eqstrn (dra,ddec); printf ("GAL2FK5: J2000 RA,Dec= %s\n",eqcoor); free (eqcoor); } return; } /* Return string with right ascension in hours and declination in degrees */ char *eqstrn (dra, ddec) double dra; /* Right ascension in degrees */ double ddec; /* Declination in degrees */ { char *eqcoor; /* ASCII character string of position (returned) */ char decp; int rah,irm,decd,decm; double xpos,ypos,xp,yp,ras,decs; /* Right ascension to hours, minutes, and seconds */ xpos = dra / 15.0; rah = (int) xpos; xp = (double) 60.0 * (xpos - (double) rah); irm = (int) xp; ras = (double) 60.0 * (xp - (double) irm); /* Declination to degrees, minutes, seconds */ if (ddec < 0) { ypos = -ddec; decp = '-'; } else { decp = '+'; ypos = ddec; } decd = (int) ypos; yp = (double) 60.0 * (ypos - (double) decd); decm = (int) yp; decs = (double) 60.0 * (yp - (double) decm); eqcoor = malloc (32); (void)sprintf (eqcoor,"%02d:%02d:%06.3f %c%02d:%02d:%05.2f", rah,irm,ras,decp,decd,decm,decs); if (eqcoor[6] == ' ') eqcoor[6] = '0'; if (eqcoor[20] == ' ') eqcoor[20] = '0'; return (eqcoor); } /* Convert geocentric equatorial rectangular coordinates to right ascension and declination, and distance */ void jpcon (pos,rra,rdec,r) double pos[3]; /* x,y,z geocentric equatorial position of object */ double *rra; /* Right ascension in radians */ double *rdec; /* Declination in radians */ double *r; /* Distance to object in same units as pos */ { double x,y,z,rxy,rxy2,z2; x = pos[0]; y = pos[1]; z = pos[2]; *rra = atan2 (y, x); if (*rra < 0.) *rra = *rra + 6.283185307179586; rxy2 = x*x + y*y; rxy = sqrt (rxy2); *rdec = atan2 (z, rxy); z2 = z * z; *r = sqrt (rxy2 + z2); return; } /* Convert right ascension, declination, and distance to geocentric equatorial rectangular coordinates */ void jpcop (rra,rdec,r,pos) double rra; /* Right ascension in radians */ double rdec; /* Declination in radians */ double r; /* Distance to object in same units as pos */ double pos[3]; /* x,y,z geocentric equatorial position of object */ { pos[0] = r * cos (rra) * cos (rdec); pos[1] = r * sin (rra) * cos (rdec); pos[2] = r * sin (rdec); return; } /* The following routines are almost verbatim from Patrick Wallace's SLALIB */ void fk4prec (ep0, ep1, ra, dec) double ep0; /* Starting Besselian epoch */ double ep1; /* Ending Besselian epoch */ double *ra; /* RA in radians mean equator & equinox of epoch ep0 mean equator & equinox of epoch ep1 (returned) */ double *dec; /* Dec in radians mean equator & equinox of epoch ep0 mean equator & equinox of epoch ep1 (returned) */ /* ** slaPreces: ** Precession - FK4 (Bessel-Newcomb, pre-IAU1976) ** ** Note: ** This routine will not correctly convert between the old and ** the new systems - for example conversion from B1950 to J2000. ** For these purposes see fk425, fk524, fk45m and fk54m. ** ** P.T.Wallace Starlink 22 December 1993 */ { double pm[3][3], v1[3], v2[3]; void mprecfk4(), slaDcs2c(), slaDmxv(), slaDcc2s(); double slaDranrm(); /* Generate appropriate precession matrix */ mprecfk4 ( ep0, ep1, pm ); /* Convert RA,Dec to x,y,z */ slaDcs2c ( *ra, *dec, v1 ); /* Precess */ slaDmxv ( pm, v1, v2 ); /* Back to RA,Dec */ slaDcc2s ( v2, ra, dec ); *ra = slaDranrm ( *ra ); } void fk5prec (ep0, ep1, ra, dec) double ep0; /* Starting epoch */ double ep1; /* Ending epoch */ double *ra; /* RA in radians mean equator & equinox of epoch ep0 mean equator & equinox of epoch ep1 (returned) */ double *dec; /* Dec in radians mean equator & equinox of epoch ep0 mean equator & equinox of epoch ep1 (returned) */ /* ** slaPreces: ** Precession - FK5 (Fricke, post-IAU1976) ** ** Note: ** This routine will not correctly convert between the old and ** the new systems - for example conversion from B1950 to J2000. ** For these purposes see fk425, fk524, fk45m and fk54m. ** ** P.T.Wallace Starlink 22 December 1993 */ { double pm[3][3], v1[3], v2[3]; void mprecfk5(), slaDcs2c(), slaDmxv(), slaDcc2s(); double slaDranrm(); /* Generate appropriate precession matrix */ mprecfk5 ( ep0, ep1, pm ); /* Convert RA,Dec to x,y,z */ slaDcs2c ( *ra, *dec, v1 ); /* Precess */ slaDmxv ( pm, v1, v2 ); /* Back to RA,Dec */ slaDcc2s ( v2, ra, dec ); *ra = slaDranrm ( *ra ); } 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 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 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; } /* 2pi */ #define D2PI 6.2831853071795864769252867665590057683943387987502 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; } /* pi/(180*3600): arcseconds to radians */ #define DAS2R 4.8481368110953599358991410235794797595635330237270e-6 void mprecfk4 (bep0, bep1, rmatp) double bep0; /* Beginning Besselian epoch */ double bep1; /* Ending Besselian epoch */ double (*rmatp)[3]; /* 3x3 Precession matrix (returned) */ /* ** slaPrebn: ** Generate the matrix of precession between two epochs, ** using the old, pre-IAU1976, Bessel-Newcomb model, using ** Kinoshita's formulation (double precision) ** ** The matrix is in the sense v(bep1) = rmatp * v(bep0) ** ** Reference: ** Kinoshita, H. (1975) 'Formulas for precession', SAO Special ** Report No. 364, Smithsonian Institution Astrophysical ** Observatory, Cambridge, Massachusetts. ** ** P.T.Wallace Starlink 30 October 1993 */ { double bigt, t, tas2r, w, zeta, z, theta; void slaDeuler(); /* Interval between basic epoch B1850.0 and beginning epoch in TC */ bigt = ( bep0 - 1850.0 ) / 100.0; /* Interval over which precession required, in tropical centuries */ t = ( bep1 - bep0 ) / 100.0; /* Euler angles */ tas2r = t * DAS2R; w = 2303.5548 + ( 1.39720 + 0.000059 * bigt ) * bigt; zeta = (w + ( 0.30242 - 0.000269 * bigt + 0.017996 * t ) * t ) * tas2r; z = (w + ( 1.09478 + 0.000387 * bigt + 0.018324 * t ) * t ) * tas2r; theta = ( 2005.1125 + ( - 0.85294 - 0.000365* bigt ) * bigt + ( - 0.42647 - 0.000365 * bigt - 0.041802 * t ) * t ) * tas2r; /* Rotation matrix */ slaDeuler ( "ZYZ", -zeta, theta, -z, rmatp ); } void mprecfk5 (ep0, ep1, rmatp) double ep0; /* Beginning epoch */ double ep1; /* Ending epoch */ double (*rmatp)[3]; /* 3x3 Precession matrix (returned) */ /* ** slaPrec: ** Form the matrix of precession between two epochs (IAU 1976, FK5). ** Notes: ** 1) The epochs are TDB (loosely ET) Julian epochs. ** 2) The matrix is in the sense v(ep1) = rmatp * v(ep0) . ** ** References: ** Lieske,J.H., 1979. Astron. Astrophys.,73,282. ** equations (6) & (7), p283. ** Kaplan,G.H., 1981. USNO circular no. 163, pa2. ** ** P.T.Wallace Starlink 31 October 1993 */ { double t0, t, tas2r, w, zeta, z, theta; void slaDeuler(); /* Interval between basic epoch J2000.0 and beginning epoch (JC) */ t0 = ( ep0 - 2000.0 ) / 100.0; /* Interval over which precession required (JC) */ t = ( ep1 - ep0 ) / 100.0; /* Euler angles */ tas2r = t * DAS2R; w = 2306.2181 + ( ( 1.39656 - ( 0.000139 * t0 ) ) * t0 ); zeta = (w + ( ( 0.30188 - 0.000344 * t0 ) + 0.017998 * t ) * t ) * tas2r; z = (w + ( ( 1.09468 + 0.000066 * t0 ) + 0.018203 * t ) * t ) * tas2r; theta = ( ( 2004.3109 + ( - 0.85330 - 0.000217 * t0 ) * t0 ) + ( ( -0.42665 - 0.000217 * t0 ) - 0.041833 * t ) * t ) * tas2r; /* Rotation matrix */ slaDeuler ( "ZYZ", -zeta, theta, -z, rmatp ); } 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 6 1995 Include stdlib.h instead of malloc.h * Apr 1 1996 Add arbitrary epoch precession * Apr 26 1996 Add FK4 <-> FK5 subroutines for use when epoch is known * Aug 6 1996 Clean up after lint */