| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palEl2ue
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Transform conventional elements into "universal" form
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palEl2ue ( double date, int jform, double epoch, double orbinc,
|
|---|
| 17 | * double anode, double perih, double aorq, double e,
|
|---|
| 18 | * double aorl, double dm, double u[13], int *jstat );
|
|---|
| 19 |
|
|---|
| 20 | * Arguments:
|
|---|
| 21 | * date = double (Given)
|
|---|
| 22 | * Epoch (TT MJD) of osculation (Note 3)
|
|---|
| 23 | * jform = int (Given)
|
|---|
| 24 | * Element set actually returned (1-3; Note 6)
|
|---|
| 25 | * epoch = double (Given)
|
|---|
| 26 | * Epoch of elements (TT MJD)
|
|---|
| 27 | * orbinc = double (Given)
|
|---|
| 28 | * inclination (radians)
|
|---|
| 29 | * anode = double (Given)
|
|---|
| 30 | * longitude of the ascending node (radians)
|
|---|
| 31 | * perih = double (Given)
|
|---|
| 32 | * longitude or argument of perihelion (radians)
|
|---|
| 33 | * aorq = double (Given)
|
|---|
| 34 | * mean distance or perihelion distance (AU)
|
|---|
| 35 | * e = double (Given)
|
|---|
| 36 | * eccentricity
|
|---|
| 37 | * aorl = double (Given)
|
|---|
| 38 | * mean anomaly or longitude (radians, JFORM=1,2 only)
|
|---|
| 39 | * dm = double (Given)
|
|---|
| 40 | * daily motion (radians, JFORM=1 only)
|
|---|
| 41 | * u = double [13] (Returned)
|
|---|
| 42 | * Universal orbital elements (Note 1)
|
|---|
| 43 | * - (0) combined mass (M+m)
|
|---|
| 44 | * - (1) total energy of the orbit (alpha)
|
|---|
| 45 | * - (2) reference (osculating) epoch (t0)
|
|---|
| 46 | * - (3-5) position at reference epoch (r0)
|
|---|
| 47 | * - (6-8) velocity at reference epoch (v0)
|
|---|
| 48 | * - (9) heliocentric distance at reference epoch
|
|---|
| 49 | * - (10) r0.v0
|
|---|
| 50 | * - (11) date (t)
|
|---|
| 51 | * - (12) universal eccentric anomaly (psi) of date, approx
|
|---|
| 52 | * jstat = int * (Returned)
|
|---|
| 53 | * status: 0 = OK
|
|---|
| 54 | * - -1 = illegal JFORM
|
|---|
| 55 | * - -2 = illegal E
|
|---|
| 56 | * - -3 = illegal AORQ
|
|---|
| 57 | * - -4 = illegal DM
|
|---|
| 58 | * - -5 = numerical error
|
|---|
| 59 |
|
|---|
| 60 | * Description:
|
|---|
| 61 | * Transform conventional osculating elements into "universal" form.
|
|---|
| 62 |
|
|---|
| 63 | * Authors:
|
|---|
| 64 | * PTW: Pat Wallace (STFC)
|
|---|
| 65 | * TIMJ: Tim Jenness (JAC, Hawaii)
|
|---|
| 66 | * {enter_new_authors_here}
|
|---|
| 67 |
|
|---|
| 68 | * Notes:
|
|---|
| 69 | * - The "universal" elements are those which define the orbit for the
|
|---|
| 70 | * purposes of the method of universal variables (see reference).
|
|---|
| 71 | * They consist of the combined mass of the two bodies, an epoch,
|
|---|
| 72 | * and the position and velocity vectors (arbitrary reference frame)
|
|---|
| 73 | * at that epoch. The parameter set used here includes also various
|
|---|
| 74 | * quantities that can, in fact, be derived from the other
|
|---|
| 75 | * information. This approach is taken to avoiding unnecessary
|
|---|
| 76 | * computation and loss of accuracy. The supplementary quantities
|
|---|
| 77 | * are (i) alpha, which is proportional to the total energy of the
|
|---|
| 78 | * orbit, (ii) the heliocentric distance at epoch, (iii) the
|
|---|
| 79 | * outwards component of the velocity at the given epoch, (iv) an
|
|---|
| 80 | * estimate of psi, the "universal eccentric anomaly" at a given
|
|---|
| 81 | * date and (v) that date.
|
|---|
| 82 | * - The companion routine is palUe2pv. This takes the set of numbers
|
|---|
| 83 | * that the present routine outputs and uses them to derive the
|
|---|
| 84 | * object's position and velocity. A single prediction requires one
|
|---|
| 85 | * call to the present routine followed by one call to palUe2pv;
|
|---|
| 86 | * for convenience, the two calls are packaged as the routine
|
|---|
| 87 | * palPlanel. Multiple predictions may be made by again calling the
|
|---|
| 88 | * present routine once, but then calling palUe2pv multiple times,
|
|---|
| 89 | * which is faster than multiple calls to palPlanel.
|
|---|
| 90 | * - DATE is the epoch of osculation. It is in the TT timescale
|
|---|
| 91 | * (formerly Ephemeris Time, ET) and is a Modified Julian Date
|
|---|
| 92 | * (JD-2400000.5).
|
|---|
| 93 | * - The supplied orbital elements are with respect to the J2000
|
|---|
| 94 | * ecliptic and equinox. The position and velocity parameters
|
|---|
| 95 | * returned in the array U are with respect to the mean equator and
|
|---|
| 96 | * equinox of epoch J2000, and are for the perihelion prior to the
|
|---|
| 97 | * specified epoch.
|
|---|
| 98 | * - The universal elements returned in the array U are in canonical
|
|---|
| 99 | * units (solar masses, AU and canonical days).
|
|---|
| 100 | * - Three different element-format options are available:
|
|---|
| 101 | *
|
|---|
| 102 | * Option JFORM=1, suitable for the major planets:
|
|---|
| 103 | *
|
|---|
| 104 | * EPOCH = epoch of elements (TT MJD)
|
|---|
| 105 | * ORBINC = inclination i (radians)
|
|---|
| 106 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 107 | * PERIH = longitude of perihelion, curly pi (radians)
|
|---|
| 108 | * AORQ = mean distance, a (AU)
|
|---|
| 109 | * E = eccentricity, e (range 0 to <1)
|
|---|
| 110 | * AORL = mean longitude L (radians)
|
|---|
| 111 | * DM = daily motion (radians)
|
|---|
| 112 | *
|
|---|
| 113 | * Option JFORM=2, suitable for minor planets:
|
|---|
| 114 | *
|
|---|
| 115 | * EPOCH = epoch of elements (TT MJD)
|
|---|
| 116 | * ORBINC = inclination i (radians)
|
|---|
| 117 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 118 | * PERIH = argument of perihelion, little omega (radians)
|
|---|
| 119 | * AORQ = mean distance, a (AU)
|
|---|
| 120 | * E = eccentricity, e (range 0 to <1)
|
|---|
| 121 | * AORL = mean anomaly M (radians)
|
|---|
| 122 | *
|
|---|
| 123 | * Option JFORM=3, suitable for comets:
|
|---|
| 124 | *
|
|---|
| 125 | * EPOCH = epoch of perihelion (TT MJD)
|
|---|
| 126 | * ORBINC = inclination i (radians)
|
|---|
| 127 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 128 | * PERIH = argument of perihelion, little omega (radians)
|
|---|
| 129 | * AORQ = perihelion distance, q (AU)
|
|---|
| 130 | * E = eccentricity, e (range 0 to 10)
|
|---|
| 131 | *
|
|---|
| 132 | * - Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are
|
|---|
| 133 | * not accessed.
|
|---|
| 134 | * - The algorithm was originally adapted from the EPHSLA program of
|
|---|
| 135 | * D.H.P.Jones (private communication, 1996). The method is based
|
|---|
| 136 | * on Stumpff's Universal Variables.
|
|---|
| 137 | *
|
|---|
| 138 | * See Also:
|
|---|
| 139 | * Everhart & Pitkin, Am.J.Phys. 51, 712 (1983).
|
|---|
| 140 |
|
|---|
| 141 | * History:
|
|---|
| 142 | * 2012-03-12 (TIMJ):
|
|---|
| 143 | * Initial version taken directly from SLA/F.
|
|---|
| 144 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 145 | * {enter_further_changes_here}
|
|---|
| 146 |
|
|---|
| 147 | * Copyright:
|
|---|
| 148 | * Copyright (C) 2005 Patrick T. Wallace
|
|---|
| 149 | * Copyright (C) 2012 Science and Technology Facilities Council.
|
|---|
| 150 | * All Rights Reserved.
|
|---|
| 151 |
|
|---|
| 152 | * Licence:
|
|---|
| 153 | * This program is free software; you can redistribute it and/or
|
|---|
| 154 | * modify it under the terms of the GNU General Public License as
|
|---|
| 155 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 156 | * the License, or (at your option) any later version.
|
|---|
| 157 | *
|
|---|
| 158 | * This program is distributed in the hope that it will be
|
|---|
| 159 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 160 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 161 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 162 | *
|
|---|
| 163 | * You should have received a copy of the GNU General Public License
|
|---|
| 164 | * along with this program; if not, write to the Free Software
|
|---|
| 165 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 166 | * MA 02110-1301, USA.
|
|---|
| 167 |
|
|---|
| 168 | * Bugs:
|
|---|
| 169 | * {note_any_bugs_here}
|
|---|
| 170 | *-
|
|---|
| 171 | */
|
|---|
| 172 |
|
|---|
| 173 | #include <math.h>
|
|---|
| 174 |
|
|---|
| 175 | #include "pal.h"
|
|---|
| 176 | #include "palmac.h"
|
|---|
| 177 |
|
|---|
| 178 | void palEl2ue ( double date, int jform, double epoch, double orbinc,
|
|---|
| 179 | double anode, double perih, double aorq, double e,
|
|---|
| 180 | double aorl, double dm, double u[13], int *jstat ) {
|
|---|
| 181 |
|
|---|
| 182 | /* Sin and cos of J2000 mean obliquity (IAU 1976) */
|
|---|
| 183 | const double SE=0.3977771559319137;
|
|---|
| 184 | const double CE=0.9174820620691818;
|
|---|
| 185 |
|
|---|
| 186 | int J;
|
|---|
| 187 |
|
|---|
| 188 | double PHT,ARGPH,Q,W,CM,ALPHA,PHS,SW,CW,SI,CI,SO,CO,
|
|---|
| 189 | X,Y,Z,PX,PY,PZ,VX,VY,VZ,DT,FC,FP,PSI,
|
|---|
| 190 | UL[13],PV[6];
|
|---|
| 191 |
|
|---|
| 192 | /* Validate arguments. */
|
|---|
| 193 | if (jform < 1 || jform > 3) {
|
|---|
| 194 | *jstat = -1;
|
|---|
| 195 | return;
|
|---|
| 196 | }
|
|---|
| 197 | if (e < 0.0 || e > 10.0 || (e >= 1.0 && jform != 3)) {
|
|---|
| 198 | *jstat = -2;
|
|---|
| 199 | return;
|
|---|
| 200 | }
|
|---|
| 201 | if (aorq <= 0.0) {
|
|---|
| 202 | *jstat = -3;
|
|---|
| 203 | return;
|
|---|
| 204 | }
|
|---|
| 205 | if (jform == 1 && dm <= 0.0) {
|
|---|
| 206 | *jstat = -4;
|
|---|
| 207 | return;
|
|---|
| 208 | }
|
|---|
| 209 |
|
|---|
| 210 | /*
|
|---|
| 211 | * Transform elements into standard form:
|
|---|
| 212 | *
|
|---|
| 213 | * PHT = epoch of perihelion passage
|
|---|
| 214 | * ARGPH = argument of perihelion (little omega)
|
|---|
| 215 | * Q = perihelion distance (q)
|
|---|
| 216 | * CM = combined mass, M+m (mu)
|
|---|
| 217 | */
|
|---|
| 218 |
|
|---|
| 219 | if (jform == 1) {
|
|---|
| 220 |
|
|---|
| 221 | /* Major planet. */
|
|---|
| 222 | PHT = epoch-(aorl-perih)/dm;
|
|---|
| 223 | ARGPH = perih-anode;
|
|---|
| 224 | Q = aorq*(1.0-e);
|
|---|
| 225 | W = dm/PAL__GCON;
|
|---|
| 226 | CM = W*W*aorq*aorq*aorq;
|
|---|
| 227 |
|
|---|
| 228 | } else if (jform == 2) {
|
|---|
| 229 |
|
|---|
| 230 | /* Minor planet. */
|
|---|
| 231 | PHT = epoch-aorl*sqrt(aorq*aorq*aorq)/PAL__GCON;
|
|---|
| 232 | ARGPH = perih;
|
|---|
| 233 | Q = aorq*(1.0-e);
|
|---|
| 234 | CM = 1.0;
|
|---|
| 235 |
|
|---|
| 236 | } else {
|
|---|
| 237 |
|
|---|
| 238 | /* Comet. */
|
|---|
| 239 | PHT = epoch;
|
|---|
| 240 | ARGPH = perih;
|
|---|
| 241 | Q = aorq;
|
|---|
| 242 | CM = 1.0;
|
|---|
| 243 |
|
|---|
| 244 | }
|
|---|
| 245 |
|
|---|
| 246 | /* The universal variable alpha. This is proportional to the total
|
|---|
| 247 | * energy of the orbit: -ve for an ellipse, zero for a parabola,
|
|---|
| 248 | * +ve for a hyperbola. */
|
|---|
| 249 |
|
|---|
| 250 | ALPHA = CM*(e-1.0)/Q;
|
|---|
| 251 |
|
|---|
| 252 | /* Speed at perihelion. */
|
|---|
| 253 |
|
|---|
| 254 | PHS = sqrt(ALPHA+2.0*CM/Q);
|
|---|
| 255 |
|
|---|
| 256 | /* In a Cartesian coordinate system which has the x-axis pointing
|
|---|
| 257 | * to perihelion and the z-axis normal to the orbit (such that the
|
|---|
| 258 | * object orbits counter-clockwise as seen from +ve z), the
|
|---|
| 259 | * perihelion position and velocity vectors are:
|
|---|
| 260 | *
|
|---|
| 261 | * position [Q,0,0]
|
|---|
| 262 | * velocity [0,PHS,0]
|
|---|
| 263 | *
|
|---|
| 264 | * To express the results in J2000 equatorial coordinates we make a
|
|---|
| 265 | * series of four rotations of the Cartesian axes:
|
|---|
| 266 | *
|
|---|
| 267 | * axis Euler angle
|
|---|
| 268 | *
|
|---|
| 269 | * 1 z argument of perihelion (little omega)
|
|---|
| 270 | * 2 x inclination (i)
|
|---|
| 271 | * 3 z longitude of the ascending node (big omega)
|
|---|
| 272 | * 4 x J2000 obliquity (epsilon)
|
|---|
| 273 | *
|
|---|
| 274 | * In each case the rotation is clockwise as seen from the +ve end of
|
|---|
| 275 | * the axis concerned.
|
|---|
| 276 | */
|
|---|
| 277 |
|
|---|
| 278 | /* Functions of the Euler angles. */
|
|---|
| 279 | SW = sin(ARGPH);
|
|---|
| 280 | CW = cos(ARGPH);
|
|---|
| 281 | SI = sin(orbinc);
|
|---|
| 282 | CI = cos(orbinc);
|
|---|
| 283 | SO = sin(anode);
|
|---|
| 284 | CO = cos(anode);
|
|---|
| 285 |
|
|---|
| 286 | /* Position at perihelion (AU). */
|
|---|
| 287 | X = Q*CW;
|
|---|
| 288 | Y = Q*SW;
|
|---|
| 289 | Z = Y*SI;
|
|---|
| 290 | Y = Y*CI;
|
|---|
| 291 | PX = X*CO-Y*SO;
|
|---|
| 292 | Y = X*SO+Y*CO;
|
|---|
| 293 | PY = Y*CE-Z*SE;
|
|---|
| 294 | PZ = Y*SE+Z*CE;
|
|---|
| 295 |
|
|---|
| 296 | /* Velocity at perihelion (AU per canonical day). */
|
|---|
| 297 | X = -PHS*SW;
|
|---|
| 298 | Y = PHS*CW;
|
|---|
| 299 | Z = Y*SI;
|
|---|
| 300 | Y = Y*CI;
|
|---|
| 301 | VX = X*CO-Y*SO;
|
|---|
| 302 | Y = X*SO+Y*CO;
|
|---|
| 303 | VY = Y*CE-Z*SE;
|
|---|
| 304 | VZ = Y*SE+Z*CE;
|
|---|
| 305 |
|
|---|
| 306 | /* Time from perihelion to date (in Canonical Days: a canonical day
|
|---|
| 307 | * is 58.1324409... days, defined as 1/PAL__GCON). */
|
|---|
| 308 |
|
|---|
| 309 | DT = (date-PHT)*PAL__GCON;
|
|---|
| 310 |
|
|---|
| 311 | /* First approximation to the Universal Eccentric Anomaly, PSI,
|
|---|
| 312 | * based on the circle (FC) and parabola (FP) values. */
|
|---|
| 313 |
|
|---|
| 314 | FC = DT/Q;
|
|---|
| 315 | W = pow(3.0*DT+sqrt(9.0*DT*DT+8.0*Q*Q*Q), 1.0/3.0);
|
|---|
| 316 | FP = W-2.0*Q/W;
|
|---|
| 317 | PSI = (1.0-e)*FC+e*FP;
|
|---|
| 318 |
|
|---|
| 319 | /* Assemble local copy of element set. */
|
|---|
| 320 | UL[0] = CM;
|
|---|
| 321 | UL[1] = ALPHA;
|
|---|
| 322 | UL[2] = PHT;
|
|---|
| 323 | UL[3] = PX;
|
|---|
| 324 | UL[4] = PY;
|
|---|
| 325 | UL[5] = PZ;
|
|---|
| 326 | UL[6] = VX;
|
|---|
| 327 | UL[7] = VY;
|
|---|
| 328 | UL[8] = VZ;
|
|---|
| 329 | UL[9] = Q;
|
|---|
| 330 | UL[10] = 0.0;
|
|---|
| 331 | UL[11] = date;
|
|---|
| 332 | UL[12] = PSI;
|
|---|
| 333 |
|
|---|
| 334 | /* Predict position+velocity at epoch of osculation. */
|
|---|
| 335 | palUe2pv( date, UL, PV, &J );
|
|---|
| 336 | if (J != 0) {
|
|---|
| 337 | *jstat = -5;
|
|---|
| 338 | return;
|
|---|
| 339 | }
|
|---|
| 340 |
|
|---|
| 341 | /* Convert back to universal elements. */
|
|---|
| 342 | palPv2ue( PV, date, CM-1.0, u, &J );
|
|---|
| 343 | if (J != 0) {
|
|---|
| 344 | *jstat = -5;
|
|---|
| 345 | return;
|
|---|
| 346 | }
|
|---|
| 347 |
|
|---|
| 348 | /* OK exit. */
|
|---|
| 349 | *jstat = 0;
|
|---|
| 350 |
|
|---|
| 351 | }
|
|---|
