diff options
Diffstat (limited to 'math/slalib/el2ue.f')
| -rw-r--r-- | math/slalib/el2ue.f | 329 |
1 files changed, 329 insertions, 0 deletions
diff --git a/math/slalib/el2ue.f b/math/slalib/el2ue.f new file mode 100644 index 00000000..4edb620a --- /dev/null +++ b/math/slalib/el2ue.f @@ -0,0 +1,329 @@ + SUBROUTINE slELUE (DATE, JFORM, EPOCH, ORBINC, ANODE, + : PERIH, AORQ, E, AORL, DM, + : U, JSTAT) +*+ +* - - - - - - +* E L U E +* - - - - - - +* +* Transform conventional osculating orbital elements into "universal" +* form. +* +* Given: +* DATE d epoch (TT MJD) of osculation (Note 3) +* JFORM i choice of element set (1-3, Note 6) +* EPOCH d epoch (TT MJD) of the elements +* ORBINC d inclination (radians) +* ANODE d longitude of the ascending node (radians) +* PERIH d longitude or argument of perihelion (radians) +* AORQ d mean distance or perihelion distance (AU) +* E d eccentricity +* AORL d mean anomaly or longitude (radians, JFORM=1,2 only) +* DM d daily motion (radians, JFORM=1 only) +* +* Returned: +* U d(13) universal orbital elements (Note 1) +* +* (1) combined mass (M+m) +* (2) total energy of the orbit (alpha) +* (3) reference (osculating) epoch (t0) +* (4-6) position at reference epoch (r0) +* (7-9) velocity at reference epoch (v0) +* (10) heliocentric distance at reference epoch +* (11) r0.v0 +* (12) date (t) +* (13) universal eccentric anomaly (psi) of date, approx +* +* JSTAT i status: 0 = OK +* -1 = illegal JFORM +* -2 = illegal E +* -3 = illegal AORQ +* -4 = illegal DM +* -5 = numerical error +* +* Called: slUEPV, slPVUE +* +* Notes +* +* 1 The "universal" elements are those which define the orbit for the +* purposes of the method of universal variables (see reference). +* They consist of the combined mass of the two bodies, an epoch, +* and the position and velocity vectors (arbitrary reference frame) +* at that epoch. The parameter set used here includes also various +* quantities that can, in fact, be derived from the other +* information. This approach is taken to avoiding unnecessary +* computation and loss of accuracy. The supplementary quantities +* are (i) alpha, which is proportional to the total energy of the +* orbit, (ii) the heliocentric distance at epoch, (iii) the +* outwards component of the velocity at the given epoch, (iv) an +* estimate of psi, the "universal eccentric anomaly" at a given +* date and (v) that date. +* +* 2 The companion routine is slUEPV. This takes the set of numbers +* that the present routine outputs and uses them to derive the +* object's position and velocity. A single prediction requires one +* call to the present routine followed by one call to slUEPV; +* for convenience, the two calls are packaged as the routine +* slPLNE. Multiple predictions may be made by again calling the +* present routine once, but then calling slUEPV multiple times, +* which is faster than multiple calls to slPLNE. +* +* 3 DATE is the epoch of osculation. It is in the TT timescale +* (formerly Ephemeris Time, ET) and is a Modified Julian Date +* (JD-2400000.5). +* +* 4 The supplied orbital elements are with respect to the J2000 +* ecliptic and equinox. The position and velocity parameters +* returned in the array U are with respect to the mean equator and +* equinox of epoch J2000, and are for the perihelion prior to the +* specified epoch. +* +* 5 The universal elements returned in the array U are in canonical +* units (solar masses, AU and canonical days). +* +* 6 Three different element-format options are available: +* +* Option JFORM=1, suitable for the major planets: +* +* EPOCH = epoch of elements (TT MJD) +* ORBINC = inclination i (radians) +* ANODE = longitude of the ascending node, big omega (radians) +* PERIH = longitude of perihelion, curly pi (radians) +* AORQ = mean distance, a (AU) +* E = eccentricity, e (range 0 to <1) +* AORL = mean longitude L (radians) +* DM = daily motion (radians) +* +* Option JFORM=2, suitable for minor planets: +* +* EPOCH = epoch of elements (TT MJD) +* ORBINC = inclination i (radians) +* ANODE = longitude of the ascending node, big omega (radians) +* PERIH = argument of perihelion, little omega (radians) +* AORQ = mean distance, a (AU) +* E = eccentricity, e (range 0 to <1) +* AORL = mean anomaly M (radians) +* +* Option JFORM=3, suitable for comets: +* +* EPOCH = epoch of perihelion (TT MJD) +* ORBINC = inclination i (radians) +* ANODE = longitude of the ascending node, big omega (radians) +* PERIH = argument of perihelion, little omega (radians) +* AORQ = perihelion distance, q (AU) +* E = eccentricity, e (range 0 to 10) +* +* 7 Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are +* not accessed. +* +* 8 The algorithm was originally adapted from the EPHSLA program of +* D.H.P.Jones (private communication, 1996). The method is based +* on Stumpff's Universal Variables. +* +* Reference: Everhart & Pitkin, Am.J.Phys. 51, 712 (1983). +* +* Last revision: 8 September 2005 +* +* Copyright P.T.Wallace. All rights reserved. +* +* License: +* This program is free software; you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation; either version 2 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 General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program (see SLA_CONDITIONS); if not, write to the +* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +* Boston, MA 02110-1301 USA +* +* Copyright (C) 1995 Association of Universities for Research in Astronomy Inc. +*- + + IMPLICIT NONE + + DOUBLE PRECISION DATE + INTEGER JFORM + DOUBLE PRECISION EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM,U(13) + INTEGER JSTAT + +* Gaussian gravitational constant (exact) + DOUBLE PRECISION GCON + PARAMETER (GCON=0.01720209895D0) + +* Sin and cos of J2000 mean obliquity (IAU 1976) + DOUBLE PRECISION SE,CE + PARAMETER (SE=0.3977771559319137D0, + : CE=0.9174820620691818D0) + + INTEGER J + + DOUBLE PRECISION PHT,ARGPH,Q,W,CM,ALPHA,PHS,SW,CW,SI,CI,SO,CO, + : X,Y,Z,PX,PY,PZ,VX,VY,VZ,DT,FC,FP,PSI, + : UL(13),PV(6) + + + +* Validate arguments. + IF (JFORM.LT.1.OR.JFORM.GT.3) THEN + JSTAT = -1 + GO TO 9999 + END IF + IF (E.LT.0D0.OR.E.GT.10D0.OR.(E.GE.1D0.AND.JFORM.NE.3)) THEN + JSTAT = -2 + GO TO 9999 + END IF + IF (AORQ.LE.0D0) THEN + JSTAT = -3 + GO TO 9999 + END IF + IF (JFORM.EQ.1.AND.DM.LE.0D0) THEN + JSTAT = -4 + GO TO 9999 + END IF + +* +* Transform elements into standard form: +* +* PHT = epoch of perihelion passage +* ARGPH = argument of perihelion (little omega) +* Q = perihelion distance (q) +* CM = combined mass, M+m (mu) + + IF (JFORM.EQ.1) THEN + +* Major planet. + PHT = EPOCH-(AORL-PERIH)/DM + ARGPH = PERIH-ANODE + Q = AORQ*(1D0-E) + W = DM/GCON + CM = W*W*AORQ*AORQ*AORQ + + ELSE IF (JFORM.EQ.2) THEN + +* Minor planet. + PHT = EPOCH-AORL*SQRT(AORQ*AORQ*AORQ)/GCON + ARGPH = PERIH + Q = AORQ*(1D0-E) + CM = 1D0 + + ELSE + +* Comet. + PHT = EPOCH + ARGPH = PERIH + Q = AORQ + CM = 1D0 + + END IF + +* The universal variable alpha. This is proportional to the total +* energy of the orbit: -ve for an ellipse, zero for a parabola, +* +ve for a hyperbola. + + ALPHA = CM*(E-1D0)/Q + +* Speed at perihelion. + + PHS = SQRT(ALPHA+2D0*CM/Q) + +* In a Cartesian coordinate system which has the x-axis pointing +* to perihelion and the z-axis normal to the orbit (such that the +* object orbits counter-clockwise as seen from +ve z), the +* perihelion position and velocity vectors are: +* +* position [Q,0,0] +* velocity [0,PHS,0] +* +* To express the results in J2000 equatorial coordinates we make a +* series of four rotations of the Cartesian axes: +* +* axis Euler angle +* +* 1 z argument of perihelion (little omega) +* 2 x inclination (i) +* 3 z longitude of the ascending node (big omega) +* 4 x J2000 obliquity (epsilon) +* +* In each case the rotation is clockwise as seen from the +ve end of +* the axis concerned. + +* Functions of the Euler angles. + SW = SIN(ARGPH) + CW = COS(ARGPH) + SI = SIN(ORBINC) + CI = COS(ORBINC) + SO = SIN(ANODE) + CO = COS(ANODE) + +* Position at perihelion (AU). + X = Q*CW + Y = Q*SW + Z = Y*SI + Y = Y*CI + PX = X*CO-Y*SO + Y = X*SO+Y*CO + PY = Y*CE-Z*SE + PZ = Y*SE+Z*CE + +* Velocity at perihelion (AU per canonical day). + X = -PHS*SW + Y = PHS*CW + Z = Y*SI + Y = Y*CI + VX = X*CO-Y*SO + Y = X*SO+Y*CO + VY = Y*CE-Z*SE + VZ = Y*SE+Z*CE + +* Time from perihelion to date (in Canonical Days: a canonical day +* is 58.1324409... days, defined as 1/GCON). + + DT = (DATE-PHT)*GCON + +* First approximation to the Universal Eccentric Anomaly, PSI, +* based on the circle (FC) and parabola (FP) values. + + FC = DT/Q + W = (3D0*DT+SQRT(9D0*DT*DT+8D0*Q*Q*Q))**(1D0/3D0) + FP = W-2D0*Q/W + PSI = (1D0-E)*FC+E*FP + +* Assemble local copy of element set. + UL(1) = CM + UL(2) = ALPHA + UL(3) = PHT + UL(4) = PX + UL(5) = PY + UL(6) = PZ + UL(7) = VX + UL(8) = VY + UL(9) = VZ + UL(10) = Q + UL(11) = 0D0 + UL(12) = DATE + UL(13) = PSI + +* Predict position+velocity at epoch of osculation. + CALL slUEPV(DATE,UL,PV,J) + IF (J.NE.0) GO TO 9010 + +* Convert back to universal elements. + CALL slPVUE(PV,DATE,CM-1D0,U,J) + IF (J.NE.0) GO TO 9010 + +* OK exit. + JSTAT = 0 + GO TO 9999 + +* Quasi-impossible numerical errors. + 9010 CONTINUE + JSTAT = -5 + + 9999 CONTINUE + END |