/* elltof.f -- translated by f2c (version 19980913). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* $Procedure ELLTOF ( Elliptic time of flight ) */ /* Subroutine */ int elltof_(doublereal *ma, doublereal *ecc, doublereal *e) { /* System generated locals */ doublereal d__1, d__2; /* Builtin functions */ double sin(doublereal), sqrt(doublereal), cos(doublereal); /* Local variables */ doublereal a, b, m; integer n; doublereal q, r__, y; extern /* Subroutine */ int chkin_(char *, ftnlen); extern doublereal dcbrt_(doublereal *); doublereal deriv, m0; extern doublereal twopi_(void); doublereal deriv2, fn, change; extern doublereal pi_(void), halfpi_(void); doublereal qr, mprime; extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *, ftnlen); extern logical return_(void); /* $ Abstract */ /* Solve the time of flight equation MA = E - e sin(E) for the */ /* elliptic eccentric anomaly E, given mean anomaly the MA and */ /* the eccentricity ECC. */ /* $ Copyright */ /* Copyright (1995), California Institute of Technology. */ /* U.S. Government sponsorship acknowledged. */ /* $ Required_Reading */ /* None. */ /* $ Keywords */ /* CONIC */ /* $ Declarations */ /* $ Brief_I/O */ /* VARIABLE I/O DESCRIPTION */ /* -------- --- -------------------------------------------------- */ /* MA I Mean anomaly at epoch. */ /* ECC I Eccentricity. */ /* E O Elliptic eccentric anomaly. */ /* $ Detailed_Input */ /* MA is the elliptic mean anomaly of an orbiting body at */ /* some epoch t, */ /* 3 1/2 */ /* MA = (t-T)(mu/a ) */ /* where T is the time of periapsis passage, a is */ /* the semi-major axis of the orbit, and mu is the */ /* gravitational parameter of the primary body. */ /* ECC is the eccentricity of the orbit. */ /* $ Detailed_Output */ /* E is the corresponding eccentric anomaly. This is the */ /* solution to the time of flight equation */ /* MA = E - e sin(E) */ /* $ Parameters */ /* None. */ /* $ Exceptions */ /* 1) If the eccentricity (ECC) is outside the range [0,1), */ /* the error 'SPICE(WRONGCONIC)' is signalled. */ /* $ Files */ /* None. */ /* $ Particulars */ /* Iterate to solve */ /* f(E,MA,e) = E - e sin(E) - MA = 0 */ /* $ Examples */ /* ELLTOF, HYPTOF, and PARTOF are used by CONICS. */ /* $ Restrictions */ /* None. */ /* $ Literature_References */ /* [1] Roger Bate, Fundamentals of Astrodynamics, Dover, 1971. */ /* [2] Ed Ng, "A General Algorithm for the Solution of Kepler's */ /* Equation for Elliptic Orbits", Cel. Mech. 20, 243, 1979. */ /* $ Author_and_Institution */ /* K.R. Gehringer (JPL) */ /* I.M. Underwood (JPL) */ /* $ Version */ /* - SPICELIB Version 1.1.0, 29-FEB-1996 (KRG) */ /* The declaration for the SPICELIB function PI is now */ /* preceded by an EXTERNAL statement declaring PI to be an */ /* external function. This removes a conflict with any */ /* compilers that have a PI intrinsic function. */ /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */ /* Comment section for permuted index source lines was added */ /* following the header. */ /* - SPICELIB Version 1.0.0, 31-JAN-1990 (IMU) */ /* -& */ /* $ Index_Entries */ /* elliptic time of flight */ /* -& */ /* $ Revisions */ /* - Beta Version 1.1.0, 8-JAN-1989 (IMU) */ /* The routine now verifies that the eccentricity is in the */ /* proper range---[0,1)---before proceeding. */ /* -& */ /* SPICELIB functions */ /* Local parameters */ /* Local variables */ /* Standard SPICE error handling. */ if (return_()) { return 0; } else { chkin_("ELLTOF", (ftnlen)6); } if (*ecc < 0. || *ecc >= 1.) { sigerr_("SPICE(WRONGCONIC)", (ftnlen)17); chkout_("ELLTOF", (ftnlen)6); return 0; } /* For reasons of numerical stability, we would like to restrict */ /* our solution to the interval [0,pi]. Because E, M, and sin E */ /* are always positive or negative together, we can pretend that M */ /* is positive and adjust the sign of the result. And for M, E > pi, */ /* we can define */ /* M = 2n pi + M' and E = 2n pi + E' */ /* where M' and E' are in the interval [-pi,pi]. Solving for E' */ /* gives us E. */ /* So, we begin by reducing the input mean anomaly to [0,pi]. */ m = abs(*ma); if (m > pi_()) { n = (integer) ((m - pi_()) / twopi_()) + 1; mprime = m - n * twopi_(); } else { n = 0; mprime = m; } m = abs(mprime); /* The convergence of the iterative scheme below depends on a good */ /* initial estimate for E. */ /* For small eccentricity, the initial estimate E = M is sufficient. */ /* However, as the eccentricity increases, so does the number of */ /* iterations required for convergence. For sufficiently large */ /* eccentricity, this estimate leads to divergence. */ /* Ng [2] notes that the function y(M,e) */ /* E - M */ /* ------- = sin(e y + M) */ /* e */ /* increases and decreases monotonically when M is in the ranges */ /* [0,M0] and [m0,pi], respectively. */ /* When M0 < M < pi, where M0 = (pi/2) - e, the cubic */ /* - - */ /* pi - M 2 pi - M pi - M */ /* B(M,e) = 1 - (1 - -------) (1 + 2 ------- - -------) */ /* pi - M0 pi - M0 1 + e */ /* provides a good initial estimate of y for all values of e. */ m0 = halfpi_() - *ecc; if (m >= m0) { a = pi_() - m; b = pi_() - m0; /* Computing 2nd power */ d__1 = 1. - a / b; y = 1. - d__1 * d__1 * (a * 2. / b + 1. - a / (*ecc + 1.)); *e = *ecc * sin(*ecc * y + m) + m; /* The situation is a little more troublesome, however, when M < M0. */ /* For small eccentricity, the cubic */ /* 2 */ /* A(M,e) = 1 - (1 - M/M0) (1 + 2M/M0 - M/(1-e) ) */ /* gives a reasonable first estimate of y. However, as e -> 1, */ /* successive approximations of the form */ /* k k */ /* C (M,e) = 1 - (-1) (1 - M/M0) */ /* k */ /* are used, where k = 4 for e > 0.7, and k = 8 for e > 0.85. */ /* For high eccentricity (e > 0.96) and low mean anomaly (M < 0.05), */ /* these successive approximations eventually fail. Fortunately, in */ /* just these cases, the cubic */ /* 3 2 1/3 3 2 1/3 */ /* D(M,e) = [r + (q + r )] + [r - (q + r )] */ /* where */ /* r = 3M/e, q = (2/e)(1 - e) */ /* provides a reasonable estimate of E directly. */ } else if (*ecc <= .7) { /* Computing 2nd power */ d__1 = 1. - m / m0; y = 1. - d__1 * d__1 * (m * 2. / m0 + 1. - m / (1. - *ecc)); *e = *ecc * sin(*ecc * y + m) + m; } else if (*ecc <= .85) { /* Computing 4th power */ d__1 = 1. - m / m0, d__1 *= d__1; y = 1. - d__1 * d__1; *e = *ecc * sin(*ecc * y + m) + m; } else if (*ecc <= .96 || m > .05) { /* Computing 8th power */ d__1 = 1. - m / m0, d__1 *= d__1, d__1 *= d__1; y = 1. - d__1 * d__1; *e = *ecc * sin(*ecc * y + m) + m; } else { q = 2. / *ecc * (1. - *ecc); r__ = m / *ecc * 3.; /* Computing 3rd power */ d__1 = q; /* Computing 2nd power */ d__2 = r__; qr = sqrt(d__1 * (d__1 * d__1) + d__2 * d__2); d__1 = r__ + qr; d__2 = r__ - qr; *e = dcbrt_(&d__1) + dcbrt_(&d__2); } /* Use the Newton second-order method, */ /* 2 */ /* E = E - (f/f')*(1 + f*f''/2f' ) */ /* i+1 i */ /* where */ /* f = E - e sin(E) - M */ /* f' = 1 - e cos(E) */ /* f'' = e sin(E) */ change = 1.; while(abs(change) > 1e-15) { fn = *e - *ecc * sin(*e) - m; deriv = 1. - *ecc * cos(*e); deriv2 = *ecc * sin(*e); /* Computing 2nd power */ d__1 = deriv; change = fn / deriv * (fn * deriv2 / (d__1 * d__1 * 2.) + 1.); *e -= change; } /* "Unwrap" E' into the actual value of E. */ if (mprime < 0.) { *e = -(*e); } if (n > 0) { *e += n * twopi_(); } if (*ma < 0.) { *e = -(*e); } chkout_("ELLTOF", (ftnlen)6); return 0; } /* elltof_ */