view release on metacpan or  search on metacpan

palsrc/palFk524.c  view on Meta::CPAN

*  History:
*     2012-02-13 (DSB):
*        Initial version with documentation taken from Fortran SLA
*        Adapted with permission from the Fortran SLALIB library.
*     {enter_further_changes_here}

*  Copyright:
*     Copyright (C) 1995 Rutherford Appleton Laboratory
*     Copyright (C) 2012 Science and Technology Facilities Council.
*     All Rights Reserved.

*  Licence:
*     This program is free software: you can redistribute it and/or
*     modify it under the terms of the GNU Lesser General Public
*     License as published by the Free Software Foundation, either
*     version 3 of the License, or (at your option) any later
*     version.
*
*     This program is distributed in the hope that it will be useful,
*     but WITHOUT ANY WARRANTY; without even the implied warranty of
*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*     GNU Lesser General Public License for more details.
*
*     You should have received a copy of the GNU Lesser General
*     License along with this program.  If not, see
*     <http://www.gnu.org/licenses/>.

*  Bugs:
*     {note_any_bugs_here}
*-
*/

#include "pal.h"
#include "palmac.h"
#include "math.h"

void palFk524( double r2000, double d2000, double dr2000, double dd2000,
               double p2000, double v2000, double *r1950, double *d1950,
               double *dr1950, double *dd1950, double *p1950, double *v1950 ){

/* Local Variables; */
   double r, d, ur, ud, px, rv;
   double sr, cr, sd, cd, x, y, z, w;
   double v1[ 6 ], v2[ 6 ];
   double xd, yd, zd;
   double rxyz, wd, rxysq, rxy;
   int i, j;

/* Small number to avoid arithmetic problems. */
   static const double tiny = 1.0E-30;

/* Canonical constants (see references). Constant vector and matrix. */
   double a[ 6 ] = { -1.62557E-6,  -0.31919E-6, -0.13843E-6,
                     +1.245E-3,    -1.580E-3,   -0.659E-3 };
   double emi[ 6 ][ 6 ] = {
                 { 0.9999256795,      0.0111814828,      0.0048590039,
                  -0.00000242389840, -0.00000002710544, -0.00000001177742},
                 {-0.0111814828,      0.9999374849,     -0.0000271771,
                   0.00000002710544, -0.00000242392702,  0.00000000006585 },
                 {-0.0048590040,     -0.0000271557,      0.9999881946,
                   0.00000001177742,  0.00000000006585, -0.00000242404995 },
                 {-0.000551,          0.238509,         -0.435614,
                   0.99990432,        0.01118145,        0.00485852 },
                 {-0.238560,         -0.002667,          0.012254,
                  -0.01118145,        0.99991613,       -0.00002717},
                 { 0.435730,         -0.008541,          0.002117,
                  -0.00485852,       -0.00002716,        0.99996684 } };

/* Pick up J2000 data (units radians and arcsec/JC). */
   r = r2000;
   d = d2000;
   ur = dr2000*PAL__PMF;
   ud = dd2000*PAL__PMF;
   px = p2000;
   rv = v2000;

/* Spherical to Cartesian. */
   sr = sin( r );
   cr = cos( r );
   sd = sin( d );
   cd = cos( d );

   x = cr*cd;
   y = sr*cd;
   z =    sd;

   w = PAL__VF*rv*px;

   v1[ 0 ] = x;
   v1[ 1 ] = y;
   v1[ 2 ] = z;

   v1[ 3 ] = -ur*y - cr*sd*ud + w*x;
   v1[ 4 ] =  ur*x - sr*sd*ud + w*y;
   v1[ 5 ] =            cd*ud + w*z;

/* Convert position+velocity vector to BN system. */
   for( i = 0; i < 6; i++ ) {
      w = 0.0;
      for( j = 0; j < 6; j++ ) {
         w += emi[ i ][ j ]*v1[ j ];
      }
      v2[ i ] = w;
   }

/* Position vector components and magnitude. */
   x = v2[ 0 ];
   y = v2[ 1 ];
   z = v2[ 2 ];
   rxyz = sqrt( x*x + y*y + z*z );

/* Apply E-terms to position. */
   w = x*a[ 0 ] + y*a[ 1 ] + z*a[ 2 ];
   x += a[ 0 ]*rxyz - w*x;
   y += a[ 1 ]*rxyz - w*y;
   z += a[ 2 ]*rxyz - w*z;

/* Recompute magnitude. */
   rxyz = sqrt( x*x + y*y + z*z );

/* Apply E-terms to both position and velocity. */