home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / AA_51.ZIP / FK4FK5.C < prev    next >
C/C++ Source or Header  |  1993-02-09  |  3KB  |  122 lines
  1.  
  2.  
  3. /* Convert FK4 B1950.0 catalogue coordinates
  4.  * to FK5 J2000.0 coordinates.
  5.  * AA page B58.
  6.  */
  7. #include "kep.h"
  8.  
  9. /* Factors to eliminate E terms of aberration
  10.  */
  11. static double A[3] = {-1.62557e-6, -3.1919e-7, - 1.3843e-7};
  12. static double Ad[3] = {1.244e-3, -1.579e-3, -6.60e-4};
  13.  
  14. /* Transformation matrix for unit direction vector,
  15.  * and motion vector in arc seconds per century
  16.  */
  17. static double Mat[36] = {
  18. 0.9999256782, -0.0111820611, -4.8579477e-3,
  19.     2.42395018e-6, -2.710663e-8, -1.177656e-8,
  20. 0.0111820610, 0.9999374784, -2.71765e-5,
  21.     2.710663e-8, 2.42397878e-6, -6.587e-11,
  22. 4.8579479e-3, -2.71474e-5, 0.9999881997,
  23.     1.177656e-8, -6.582e-11, 2.42410173e-6,
  24. -5.51e-4, -0.238565, 0.435739,
  25.     0.99994704, -0.01118251, -4.85767e-3,
  26. 0.238514, -2.667e-3, -8.541e-3,
  27.     0.01118251, 0.99995883, -2.718e-5,
  28. -0.435623, 0.012254, 2.117e-3,
  29.     4.85767e-3, -2.714e-5, 1.00000956
  30. };
  31.  
  32.  
  33. int fk4fk5( p, m, el )
  34. double p[], m[];
  35. struct star *el;
  36. {
  37. double a, b, c;
  38. double *u, *v;
  39. double R[6];
  40. int i, j;
  41.  
  42. printf( "Converting to FK5 system\n" );
  43.  
  44. /* Note the direction vector and motion vector
  45.  * are already supplied by rstar.c.
  46.  */
  47. a = 0.0;
  48. b = 0.0;
  49. for( i=0; i<3; i++ )
  50.     {
  51.     m[i] *= RTS; /* motion must be in arc seconds per century */
  52.     a += A[i] * p[i];
  53.     b += Ad[i] * p[i];
  54.     }
  55. /* Remove E terms of aberration from FK4
  56.  */
  57. for( i=0; i<3; i++ )
  58.     {
  59.     R[i] = p[i] - A[i] + a * p[i];
  60.     R[i+3] = m[i] - Ad[i] + b * p[i];
  61.     }
  62.  
  63. /* Perform matrix multiplication
  64.  */
  65. v = &Mat[0];
  66. for( i=0; i<6; i++ )
  67.     {
  68.     a = 0.0;
  69.     u = &R[0];
  70.     for( j=0; j<6; j++ )
  71.         a += *u++ * *v++;
  72.     if( i < 3 )
  73.         p[i] = a;
  74.     else
  75.         m[i-3] = a;
  76.     }
  77.  
  78. /* Transform the answers into J2000 catalogue entries
  79.  * in radian measure.
  80.  */
  81.  
  82. b = p[0]*p[0] + p[1]*p[1];
  83. a = b + p[2]*p[2];
  84. c = a;
  85. a = sqrt(a);
  86.  
  87. el->ra = zatan2( p[0], p[1] );
  88. el->dec = asin( p[2]/a );
  89.  
  90. /* Note motion converted back to radians per (Julian) century */
  91. el->mura = (p[0]*m[1] - p[1]*m[0])/(RTS*b);
  92. el->mudec = (m[2]*b - p[2]*(p[0]*m[0] + p[1]*m[1]) )/(RTS*c*sqrt(b));
  93.  
  94. if( el->px > 0.0 )
  95.     {
  96.     c = 0.0;
  97.     for( i=0; i<3; i++ )
  98.         c += p[i] * m[i];
  99.  
  100. /* divide by RTS to deconvert m (and therefore c)
  101.  * from arc seconds back to radians
  102.  */
  103.     el->v = c/(21.094952663 * el->px * RTS * a);
  104.     }
  105. el->px = el->px/a;    /* a is dimensionless */
  106. el->epoch = J2000;
  107.  
  108. #if DEBUG
  109. /* Display the computed J2000 catalogue entries
  110.  */
  111. hms( el->ra );
  112. dms( el->dec );
  113.  
  114. u = (double *)&el->px;
  115. for( i=0; i<3; i++ )
  116.     printf( " %.4lf ", *u++ * RTS );
  117. printf( "\n" );
  118. #endif
  119. return(0);
  120. }
  121.  
  122.