home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / AA_51.ZIP / KEPLER.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  10KB  |  410 lines
  1.  
  2. /* Program to solve Keplerian orbit
  3.  * given orbital parameters and the time.
  4.  * Returns Heliocentric equatorial rectangular coordinates of
  5.  * the object.
  6.  *
  7.  * This program detects several cases of given orbital elements.
  8.  * If a program for perturbations is pointed to, it is called
  9.  * to calculate all the elements.
  10.  * If there is no program, then the mean longitude is calculated
  11.  * from the mean anomaly and daily motion.
  12.  * If the daily motion is not given, it is calculated
  13.  * by Kepler's law.
  14.  * If the eccentricity is given to be 1.0, it means that
  15.  * meandistance is really the perihelion distance, as in a comet
  16.  * specification, and the orbit is parabolic.
  17.  *
  18.  * Reference: Taff, L.G., "Celestial Mechanics, A Computational
  19.  * Guide for the Practitioner."  Wiley, 1985.
  20.  */
  21.  
  22. #include "kep.h"
  23.  
  24. #if !BandS
  25. extern struct orbit earth;    /* orbital elements of the earth */
  26. #endif
  27.  
  28. extern double eps, coseps, sineps; /* obliquity of ecliptic */
  29. double sinh(), cosh(), tanh();
  30. int embofs();
  31.  
  32. int kepler(J, e, rect, polar)
  33. double J, rect[], polar[];
  34. struct orbit *e;
  35. {
  36. double alat, E, M, W, temp;
  37. double epoch, inclination, ascnode, argperih;
  38. double meandistance, dailymotion, eccent, meananomaly;
  39. double r, coso, sino, cosa;
  40. int k;
  41.  
  42. /* Compute orbital elements if a program for doing so
  43.  * is supplied
  44.  */
  45. if( e->oelmnt )
  46.     {
  47.     k = (*(e->oelmnt) )(e,J);
  48.     if( k == -1 )
  49.         goto dobs;
  50.     }
  51. else if( e->celmnt )
  52.     { /* call B & S algorithm */
  53. dobs:
  54.     (*(e->celmnt) )( J, polar );
  55.     polar[0] = modtp( polar[0] );
  56.     E = polar[0]; /* longitude */
  57.     e->L = E;
  58.     W = polar[1]; /* latitude */
  59.     r = polar[2]; /* radius */
  60.     e->r = r;
  61.     e->epoch = J;
  62.     e->equinox = J;
  63.     goto kepdon;
  64.     }
  65.  
  66. /* Decant the parameters from the data structure
  67.  */
  68. epoch = e->epoch;
  69. inclination = e->i;
  70. ascnode = e->W * DTR;
  71. argperih = e->w;
  72. meandistance = e->a; /* semimajor axis */
  73. dailymotion = e->dm;
  74. eccent = e->ecc;
  75. meananomaly = e->M;
  76. /* Check for parabolic orbit. */
  77. if( eccent == 1.0 )
  78.     {
  79. /* meandistance = perihelion distance, q
  80.  * epoch = perihelion passage date
  81.  */
  82.     temp = meandistance * sqrt(meandistance);
  83.     W = (J - epoch ) * 0.0364911624 / temp;
  84. /* The constant above is 3 k / sqrt(2),
  85.  * k = Gaussian gravitational constant = 0.01720209895 . */
  86.     E = 0.0;
  87.     M = 1.0;
  88.     while( fabs(M) > 1.0e-11 )
  89.         {
  90.         temp = E * E;
  91.         temp = (2.0 * E * temp + W)/( 3.0 * (1.0 + temp));
  92.         M = temp - E;
  93.         if( temp != 0.0 )
  94.             M /= temp;
  95.         E = temp;
  96.         }
  97.     r = meandistance * (1.0 + E * E );
  98.     M = atan( E );
  99.     M = 2.0 * M;
  100.     alat = M + DTR*argperih;
  101.     goto parabcon;
  102.     }
  103. if( eccent > 1.0 )
  104.     {
  105. /* The equation of the hyperbola in polar coordinates r, theta
  106.  * is r = a(e^2 - 1)/(1 + e cos(theta))
  107.  * so the perihelion distance q = a(e-1),
  108.  * the "mean distance"  a = q/(e-1).
  109.  */
  110.     meandistance = meandistance/(eccent - 1.0);
  111.     temp = meandistance * sqrt(meandistance);
  112.     W = (J - epoch ) * 0.01720209895 / temp;
  113. /* solve M = -E + e sinh E */
  114.     E = W/(eccent - 1.0);
  115.     M = 1.0;
  116.     while( fabs(M) > 1.0e-11 )
  117.         {
  118.         M = -E + eccent * sinh(E) - W;
  119.         E += M/(1.0 - eccent * cosh(E));
  120.         }
  121.     r = meandistance * (-1.0 + eccent * cosh(E));
  122.     temp = (eccent + 1.0)/(eccent - 1.0);
  123.     M = sqrt(temp) * tanh( 0.5*E );
  124.     M = 2.0 * atan(M);    
  125.     alat = M + DTR*argperih;
  126.     goto parabcon;
  127.     }
  128. /* Calculate the daily motion, if it is not given.
  129.  */
  130. if( dailymotion == 0.0 )
  131.     {
  132. /* The constant is 180 k / pi, k = Gaussian gravitational constant.
  133.  * Assumes object in heliocentric orbit is massless.
  134.  */
  135.     dailymotion = 0.9856076686/(e->a*sqrt(e->a));
  136.     }
  137. dailymotion *= J - epoch;
  138. /* M is proportional to the area swept out by the radius
  139.  * vector of a circular orbit during the time between
  140.  * perihelion passage and Julian date J.
  141.  * It is the mean anomaly at time J.
  142.  */
  143. M = DTR*( meananomaly + dailymotion );
  144. M = modtp(M);
  145. /* If mean longitude was calculated, adjust it also
  146.  * for motion since epoch of elements.
  147.  */
  148. if( e->L )
  149.     {
  150.     e->L += dailymotion;
  151.     e->L = mod360( e->L );
  152.     }
  153.  
  154. /* By Kepler's second law, M must be equal to
  155.  * the area swept out in the same time by an
  156.  * elliptical orbit of same total area.
  157.  * Integrate the ellipse expressed in polar coordinates
  158.  *     r = a(1-e^2)/(1 + e cosW)
  159.  * with respect to the angle W to get an expression for the
  160.  * area swept out by the radius vector.  The area is given
  161.  * by the mean anomaly; the angle is solved numerically.
  162.  * 
  163.  * The answer is obtained in two steps.  We first solve
  164.  * Kepler's equation
  165.  *    M = E - eccent*sin(E)
  166.  * for the eccentric anomaly E.  Then there is a
  167.  * closed form solution for W in terms of E.
  168.  */
  169.  
  170. E = M; /* Initial guess is same as circular orbit. */
  171. temp = 1.0;
  172. do
  173.     {
  174. /* The approximate area swept out in the ellipse */
  175.     temp = E - eccent * sin(E)
  176. /* ...minus the area swept out in the circle */
  177.         - M;
  178. /* ...should be zero.  Use the derivative of the error
  179.  * to converge to solution by Newton's method.
  180.  */
  181.     E -= temp/(1.0 - eccent*cos(E));
  182.     }
  183. while( fabs(temp) > 1.0e-11 );
  184.  
  185. /* The exact formula for the area in the ellipse is
  186.  *    2.0*atan(c2*tan(0.5*W)) - c1*eccent*sin(W)/(1+e*cos(W))
  187.  * where
  188.  *    c1 = sqrt( 1.0 - eccent*eccent )
  189.  *    c2 = sqrt( (1.0-eccent)/(1.0+eccent) ).
  190.  * Substituting the following value of W
  191.  * yields the exact solution.
  192.  */
  193. temp = sqrt( (1.0+eccent)/(1.0-eccent) );
  194. W = 2.0 * atan( temp * tan(0.5*E) );
  195.  
  196. /* The true anomaly.
  197.  */
  198. W = modtp(W);
  199.  
  200. meananomaly *= DTR;
  201. /* Orbital longitude measured from node
  202.  * (argument of latitude)
  203.  */
  204. if( e->L )
  205.     alat = (e->L)*DTR + W - meananomaly - ascnode;
  206. else
  207.     alat = W + DTR*argperih; /* mean longitude not given */
  208.  
  209. /* From the equation of the ellipse, get the
  210.  * radius from central focus to the object.
  211.  */
  212. r = meandistance*(1.0-eccent*eccent)/(1.0+eccent*cos(W));
  213.  
  214. parabcon:
  215.  
  216. /* The heliocentric ecliptic longitude of the object
  217.  * is given by
  218.  *   tan( longitude - ascnode )  =  cos( inclination ) * tan( alat ).
  219.  */
  220. coso = cos( alat );
  221. sino = sin( alat );
  222. inclination *= DTR;
  223. W = sino * cos( inclination );
  224. E = zatan2( coso, W ) + ascnode;
  225.  
  226. /* The ecliptic latitude of the object
  227.  */
  228. W = sino * sin( inclination );
  229. W = asin(W);
  230.  
  231. /* Apply perturbations, if a program is supplied.
  232.  */
  233. if( e->celmnt )
  234.     {
  235.     e->L = E;
  236.     e->r = r;
  237.     (*(e->celmnt) )(e);
  238.     E = e->L;
  239.     r = e->r;
  240.     W += e->plat;
  241.     }
  242.  
  243. /* If earth, Adjust from earth-moon barycenter to earth
  244.  * by AA page E2
  245.  * unless orbital elements are calculated by formula.
  246.  * (The Meeus perturbation formulas include this term for the moon.)
  247.  */
  248. #if !BandS
  249. if( (e == &earth) && (e->oelmnt == 0) )
  250.     {
  251.     embofs( J, &E, &W, &r ); /* see below */
  252.     }
  253. #endif
  254.  
  255. /* Output the polar cooordinates
  256.  */
  257. polar[0] = E; /* longitude */
  258. polar[1] = W; /* latitude */
  259. polar[2] = r; /* radius */
  260.  
  261.  
  262. kepdon:
  263.  
  264.  
  265. /* Convert to rectangular coordinates,
  266.  * using the perturbed latitude.
  267.  */
  268. rect[2] = r * sin(W);
  269. cosa = cos(W);
  270. rect[1] = r * cosa * sin(E);
  271. rect[0] = r * cosa * cos(E);
  272.  
  273. /* Convert from heliocentric ecliptic rectangular
  274.  * to heliocentric equatorial rectangular coordinates
  275.  * by rotating eps radians about the x axis.
  276.  */
  277. epsiln( e->equinox );
  278. W = coseps*rect[1] - sineps*rect[2];
  279. M = sineps*rect[1] + coseps*rect[2];
  280. rect[1] = W;
  281. rect[2] = M;
  282.  
  283. /* Precess the position
  284.  * to ecliptic and equinox of J2000.0
  285.  * if not already there.
  286.  */
  287. precess( rect, e->equinox, 1 );
  288. return(0);
  289. }
  290.  
  291.  
  292. #if !BandS
  293. /* Adjust position from Earth-Moon barycenter to Earth
  294.  *
  295.  * J = Julian day number
  296.  * E = Earth's ecliptic longitude (radians)
  297.  * W = Earth's ecliptic latitude (radians)
  298.  * r = Earth's distance to the Sun (au)
  299.  */
  300. int embofs( J, E, W, r )
  301. double J;
  302. double *E, *W, *r;
  303. {
  304. double T, M, a, L, B, p;
  305. double smp, cmp, s2mp, c2mp, s2d, c2d, sf, cf;
  306. double s2f, sx, cx, xyz[3];
  307.  
  308. /* Short series for position of the Moon
  309.  */
  310. T = (J-J1900)/36525.0;
  311. /* Mean anomaly of moon (MP) */
  312. a = ((1.44e-5*T + 0.009192)*T + 477198.8491)*T + 296.104608;
  313. a *= DTR;
  314. smp = sin(a);
  315. cmp = cos(a);
  316. s2mp = 2.0*smp*cmp;        /* sin(2MP) */
  317. c2mp = cmp*cmp - smp*smp;    /* cos(2MP) */
  318. /* Mean elongation