home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / AA_51.ZIP / RPLANET.C < prev    next >
C/C++ Source or Header  |  1993-02-09  |  3KB  |  123 lines
  1.  
  2. /* The following program reduces the heliocentric equatorial
  3.  * rectangular coordinates of the earth and object that
  4.  * were computed by kepler() and produces apparent geocentric
  5.  * right ascension and declination.
  6.  */
  7.  
  8. #include "kep.h"
  9.  
  10. int reduce( elemnt, q, e )
  11. struct orbit *elemnt;    /* orbital elements of q */
  12. double q[], e[];    /* heliocentric coordinates */
  13. {
  14. double p[3], temp[3], polar[3];
  15. double a, b, s;
  16. int i;
  17. double sqrt(), asin(), log();
  18.  
  19. /* Save the geometric coordinates at TDT
  20.  */
  21. for( i=0; i<3; i++ )
  22.     temp[i] = q[i];
  23.  
  24. /* Display ecliptic longitude and latitude
  25.  */
  26. if( prtflg )
  27.     lonlat( q, TDT, polar );
  28.  
  29. /* Adjust for light time (planetary aberration)
  30. */
  31. lightt( elemnt, q, e );
  32.  
  33. /* Find Euclidean vectors between earth, object, and the sun
  34.  */
  35. for( i=0; i<3; i++ )
  36.     p[i] = q[i] - e[i];
  37.  
  38. angles( p, q, e );
  39.  
  40. if( prtflg )
  41.     {
  42.     a = 0.0;
  43.     for( i=0; i<3; i++ )
  44.         {
  45.         b = temp[i] - e[i];
  46.         a += b * b;
  47.         }
  48.     a = sqrt(a);
  49.     printf( "true geocentric distance %.7f au    ", a ); /* was EO */
  50.     printf( "equatorial diameter %.2f\"\n", 2.0*elemnt->sdiam/EO );
  51.  
  52.  
  53. /* Calculate visual magnitude.
  54.  * "Visual" refers to the spectrum of visible light.
  55.  * Phase = 0.5(1+pq) = geometric fraction of disc illuminated.
  56.  * where pq = cos( sun-object-earth angle )
  57.  * The magnitude is
  58.  *    V(1,0) + 2.5 log10( SE^2 SO^2 / Phase)
  59.  * where V(1,0) = elemnt->mag is the magnitude at 1au from
  60.  * both earth and sun and 100% illumination.
  61.  */
  62.     a = 0.5 * (1.0 + pq);
  63. /* Fudge the phase for light leakage in magnitude estimation.
  64.  * Note this phase term estimate does not reflect reality well.
  65.  * Calculated magnitudes of Mercury and Venus are inaccurate.
  66.  */
  67.     b = 0.5 * (1.01 + 0.99*pq);
  68.     s = elemnt->mag + 2.1715 * log( EO*SO ) - 1.085*log(b);
  69.     printf( "approx. visual magnitude %.1f, phase %.3f\n", s, a );
  70.     }
  71.  
  72. /* Find unit vector from earth in direction of object
  73.  */
  74. for( i=0; i<3; i++ )
  75.     {
  76.     p[i] /= EO;
  77.     temp[i] = p[i];
  78.     }
  79.  
  80. if( prtflg )
  81.     {
  82. /* Report astrometric position
  83.  */
  84.     showrd( "Astrometric J2000.0:", p, polar );
  85.  
  86. /* Also in 1950 coordinates
  87.  */
  88.     precess( temp, B1950, -1 );
  89.     showrd( "Astrometric B1950.0:", temp, polar );
  90.     }
  91.  
  92. /* Correct position for light deflection
  93.  */
  94. relativity( p, q, e );
  95.  
  96. /* Correct for annual aberration
  97.  */
  98. annuab( p );
  99.  
  100. /* Precession of the equinox and ecliptic
  101.  * from J2000.0 to ephemeris date
  102.  */
  103. precess( p, TDT, -1 );
  104.  
  105. /* Ajust for nutation
  106.  * at current ecliptic.
  107.  */
  108. epsiln( TDT );
  109. nutate( TDT, p );
  110.  
  111.  
  112. /* Display the final apparent R.A. and Dec.
  113.  * for equinox of date.
  114.  */
  115. showrd( "    Apparent:", p, polar );
  116.  
  117. /* Go do topocentric reductions.
  118.  */
  119. polar[2] = EO;
  120. altaz( polar, UT );
  121. return(0);
  122. }
  123.