/* *+ * Name: * palEl2ue * Purpose: * Transform conventional elements into "universal" form * Language: * Starlink ANSI C * Type of Module: * Library routine * Invocation: * void palEl2ue ( double date, int jform, double epoch, double orbinc, * double anode, double perih, double aorq, double e, * double aorl, double dm, double u[13], int *jstat ); * Arguments: * date = double (Given) * Epoch (TT MJD) of osculation (Note 3) * jform = int (Given) * Element set actually returned (1-3; Note 6) * epoch = double (Given) * Epoch of elements (TT MJD) * orbinc = double (Given) * inclination (radians) * anode = double (Given) * longitude of the ascending node (radians) * perih = double (Given) * longitude or argument of perihelion (radians) * aorq = double (Given) * mean distance or perihelion distance (AU) * e = double (Given) * eccentricity * aorl = double (Given) * mean anomaly or longitude (radians, JFORM=1,2 only) * dm = double (Given) * daily motion (radians, JFORM=1 only) * u = double [13] (Returned) * Universal orbital elements (Note 1) * - (0) combined mass (M+m) * - (1) total energy of the orbit (alpha) * - (2) reference (osculating) epoch (t0) * - (3-5) position at reference epoch (r0) * - (6-8) velocity at reference epoch (v0) * - (9) heliocentric distance at reference epoch * - (10) r0.v0 * - (11) date (t) * - (12) universal eccentric anomaly (psi) of date, approx * jstat = int * (Returned) * status: 0 = OK * - -1 = illegal JFORM * - -2 = illegal E * - -3 = illegal AORQ * - -4 = illegal DM * - -5 = numerical error * Description: * Transform conventional osculating elements into "universal" form. * Authors: * PTW: Pat Wallace (STFC) * TIMJ: Tim Jenness (JAC, Hawaii) * {enter_new_authors_here} * Notes: * - 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. * - The companion routine is palUe2pv. 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 palUe2pv; * for convenience, the two calls are packaged as the routine * palPlanel. Multiple predictions may be made by again calling the * present routine once, but then calling palUe2pv multiple times, * which is faster than multiple calls to palPlanel. * - 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). * - 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. * - The universal elements returned in the array U are in canonical * units (solar masses, AU and canonical days). * - 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) * * - Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are * not accessed. * - 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. * * See Also: * Everhart & Pitkin, Am.J.Phys. 51, 712 (1983). * History: * 2012-03-12 (TIMJ): * Initial version taken directly from SLA/F. * Adapted with permission from the Fortran SLALIB library. * {enter_further_changes_here} * Copyright: * Copyright (C) 2005 Patrick T. Wallace * Copyright (C) 2012 Science and Technology Facilities Council. * All Rights Reserved. * Licence: * 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 3 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; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, * MA 02110-1301, USA. * Bugs: * {note_any_bugs_here} *- */ #include #include "pal.h" #include "palmac.h" void palEl2ue ( double date, int jform, double epoch, double orbinc, double anode, double perih, double aorq, double e, double aorl, double dm, double u[13], int *jstat ) { /* Sin and cos of J2000 mean obliquity (IAU 1976) */ const double SE=0.3977771559319137; const double CE=0.9174820620691818; int J; double 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 < 1 || jform > 3) { *jstat = -1; return; } if (e < 0.0 || e > 10.0 || (e >= 1.0 && jform != 3)) { *jstat = -2; return; } if (aorq <= 0.0) { *jstat = -3; return; } if (jform == 1 && dm <= 0.0) { *jstat = -4; return; } /* * 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 == 1) { /* Major planet. */ PHT = epoch-(aorl-perih)/dm; ARGPH = perih-anode; Q = aorq*(1.0-e); W = dm/PAL__GCON; CM = W*W*aorq*aorq*aorq; } else if (jform == 2) { /* Minor planet. */ PHT = epoch-aorl*sqrt(aorq*aorq*aorq)/PAL__GCON; ARGPH = perih; Q = aorq*(1.0-e); CM = 1.0; } else { /* Comet. */ PHT = epoch; ARGPH = perih; Q = aorq; CM = 1.0; } /* 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-1.0)/Q; /* Speed at perihelion. */ PHS = sqrt(ALPHA+2.0*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/PAL__GCON). */ DT = (date-PHT)*PAL__GCON; /* First approximation to the Universal Eccentric Anomaly, PSI, * based on the circle (FC) and parabola (FP) values. */ FC = DT/Q; W = pow(3.0*DT+sqrt(9.0*DT*DT+8.0*Q*Q*Q), 1.0/3.0); FP = W-2.0*Q/W; PSI = (1.0-e)*FC+e*FP; /* Assemble local copy of element set. */ UL[0] = CM; UL[1] = ALPHA; UL[2] = PHT; UL[3] = PX; UL[4] = PY; UL[5] = PZ; UL[6] = VX; UL[7] = VY; UL[8] = VZ; UL[9] = Q; UL[10] = 0.0; UL[11] = date; UL[12] = PSI; /* Predict position+velocity at epoch of osculation. */ palUe2pv( date, UL, PV, &J ); if (J != 0) { *jstat = -5; return; } /* Convert back to universal elements. */ palPv2ue( PV, date, CM-1.0, u, &J ); if (J != 0) { *jstat = -5; return; } /* OK exit. */ *jstat = 0; }