| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palUe2pv
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Heliocentric position and velocity of a planet, asteroid or comet, from universal elements
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palUe2pv( double date, double u[13], double pv[6], int *jstat );
|
|---|
| 17 |
|
|---|
| 18 | * Arguments:
|
|---|
| 19 | * date = double (Given)
|
|---|
| 20 | * TT Modified Julian date (JD-2400000.5).
|
|---|
| 21 | * u = double [13] (Given & Returned)
|
|---|
| 22 | * Universal orbital elements (updated, see note 1)
|
|---|
| 23 | * given (0) combined mass (M+m)
|
|---|
| 24 | * " (1) total energy of the orbit (alpha)
|
|---|
| 25 | * " (2) reference (osculating) epoch (t0)
|
|---|
| 26 | * " (3-5) position at reference epoch (r0)
|
|---|
| 27 | * " (6-8) velocity at reference epoch (v0)
|
|---|
| 28 | * " (9) heliocentric distance at reference epoch
|
|---|
| 29 | * " (10) r0.v0
|
|---|
| 30 | * returned (11) date (t)
|
|---|
| 31 | * " (12) universal eccentric anomaly (psi) of date
|
|---|
| 32 | * pv = double [6] (Returned)
|
|---|
| 33 | * Position (AU) and velocity (AU/s)
|
|---|
| 34 | * jstat = int * (Returned)
|
|---|
| 35 | * status: 0 = OK
|
|---|
| 36 | * -1 = radius vector zero
|
|---|
| 37 | * -2 = failed to converge
|
|---|
| 38 |
|
|---|
| 39 | * Description:
|
|---|
| 40 | * Heliocentric position and velocity of a planet, asteroid or comet,
|
|---|
| 41 | * starting from orbital elements in the "universal variables" form.
|
|---|
| 42 |
|
|---|
| 43 | * Authors:
|
|---|
| 44 | * PTW: Pat Wallace (STFC)
|
|---|
| 45 | * TIMJ: Tim Jenness (JAC, Hawaii)
|
|---|
| 46 | * {enter_new_authors_here}
|
|---|
| 47 |
|
|---|
| 48 | * Notes:
|
|---|
| 49 | * - The "universal" elements are those which define the orbit for the
|
|---|
| 50 | * purposes of the method of universal variables (see reference).
|
|---|
| 51 | * They consist of the combined mass of the two bodies, an epoch,
|
|---|
| 52 | * and the position and velocity vectors (arbitrary reference frame)
|
|---|
| 53 | * at that epoch. The parameter set used here includes also various
|
|---|
| 54 | * quantities that can, in fact, be derived from the other
|
|---|
| 55 | * information. This approach is taken to avoiding unnecessary
|
|---|
| 56 | * computation and loss of accuracy. The supplementary quantities
|
|---|
| 57 | * are (i) alpha, which is proportional to the total energy of the
|
|---|
| 58 | * orbit, (ii) the heliocentric distance at epoch, (iii) the
|
|---|
| 59 | * outwards component of the velocity at the given epoch, (iv) an
|
|---|
| 60 | * estimate of psi, the "universal eccentric anomaly" at a given
|
|---|
| 61 | * date and (v) that date.
|
|---|
| 62 | * - The companion routine is palEl2ue. This takes the conventional
|
|---|
| 63 | * orbital elements and transforms them into the set of numbers
|
|---|
| 64 | * needed by the present routine. A single prediction requires one
|
|---|
| 65 | * one call to palEl2ue followed by one call to the present routine;
|
|---|
| 66 | * for convenience, the two calls are packaged as the routine
|
|---|
| 67 | * palPlanel. Multiple predictions may be made by again
|
|---|
| 68 | * calling palEl2ue once, but then calling the present routine
|
|---|
| 69 | * multiple times, which is faster than multiple calls to palPlanel.
|
|---|
| 70 | * - It is not obligatory to use palEl2ue to obtain the parameters.
|
|---|
| 71 | * However, it should be noted that because palEl2ue performs its
|
|---|
| 72 | * own validation, no checks on the contents of the array U are made
|
|---|
| 73 | * by the present routine.
|
|---|
| 74 | * - DATE is the instant for which the prediction is required. It is
|
|---|
| 75 | * in the TT timescale (formerly Ephemeris Time, ET) and is a
|
|---|
| 76 | * Modified Julian Date (JD-2400000.5).
|
|---|
| 77 | * - The universal elements supplied in the array U are in canonical
|
|---|
| 78 | * units (solar masses, AU and canonical days). The position and
|
|---|
| 79 | * velocity are not sensitive to the choice of reference frame. The
|
|---|
| 80 | * palEl2ue routine in fact produces coordinates with respect to the
|
|---|
| 81 | * J2000 equator and equinox.
|
|---|
| 82 | * - The algorithm was originally adapted from the EPHSLA program of
|
|---|
| 83 | * D.H.P.Jones (private communication, 1996). The method is based
|
|---|
| 84 | * on Stumpff's Universal Variables.
|
|---|
| 85 | * - Reference: Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
|
|---|
| 86 |
|
|---|
| 87 | * History:
|
|---|
| 88 | * 2012-03-09 (TIMJ):
|
|---|
| 89 | * Initial version cloned from SLA/F.
|
|---|
| 90 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 91 | * {enter_further_changes_here}
|
|---|
| 92 |
|
|---|
| 93 | * Copyright:
|
|---|
| 94 | * Copyright (C) 2005 Rutherford Appleton Laboratory
|
|---|
| 95 | * Copyright (C) 2012 Science and Technology Facilities Council.
|
|---|
| 96 | * All Rights Reserved.
|
|---|
| 97 |
|
|---|
| 98 | * Licence:
|
|---|
| 99 | * This program is free software; you can redistribute it and/or
|
|---|
| 100 | * modify it under the terms of the GNU General Public License as
|
|---|
| 101 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 102 | * the License, or (at your option) any later version.
|
|---|
| 103 | *
|
|---|
| 104 | * This program is distributed in the hope that it will be
|
|---|
| 105 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 106 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 107 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 108 | *
|
|---|
| 109 | * You should have received a copy of the GNU General Public License
|
|---|
| 110 | * along with this program; if not, write to the Free Software
|
|---|
| 111 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 112 | * MA 02110-1301, USA.
|
|---|
| 113 |
|
|---|
| 114 | * Bugs:
|
|---|
| 115 | * {note_any_bugs_here}
|
|---|
| 116 | *-
|
|---|
| 117 | */
|
|---|
| 118 |
|
|---|
| 119 | #include <math.h>
|
|---|
| 120 |
|
|---|
| 121 | #include "pal.h"
|
|---|
| 122 | #include "palmac.h"
|
|---|
| 123 |
|
|---|
| 124 | void palUe2pv( double date, double u[13], double pv[6], int *jstat ) {
|
|---|
| 125 |
|
|---|
| 126 | /* Canonical days to seconds */
|
|---|
| 127 | const double CD2S = PAL__GCON / PAL__SPD;
|
|---|
| 128 |
|
|---|
| 129 | /* Test value for solution and maximum number of iterations */
|
|---|
| 130 | const double TEST = 1e-13;
|
|---|
| 131 | const int NITMAX = 25;
|
|---|
| 132 |
|
|---|
| 133 | int I, NIT, N;
|
|---|
| 134 | double CM,ALPHA,T0,P0[3],V0[3],R0,SIGMA0,T,PSI,DT,W,
|
|---|
| 135 | TOL,PSJ,PSJ2,BETA,S0,S1,S2,S3,
|
|---|
| 136 | FF,R,F,G,FD,GD;
|
|---|
| 137 |
|
|---|
| 138 | double PLAST = 0.0;
|
|---|
| 139 | double FLAST = 0.0;
|
|---|
| 140 |
|
|---|
| 141 | /* Unpack the parameters. */
|
|---|
| 142 | CM = u[0];
|
|---|
| 143 | ALPHA = u[1];
|
|---|
| 144 | T0 = u[2];
|
|---|
| 145 | for (I=0; I<3; I++) {
|
|---|
| 146 | P0[I] = u[I+3];
|
|---|
| 147 | V0[I] = u[I+6];
|
|---|
| 148 | }
|
|---|
| 149 | R0 = u[9];
|
|---|
| 150 | SIGMA0 = u[10];
|
|---|
| 151 | T = u[11];
|
|---|
| 152 | PSI = u[12];
|
|---|
| 153 |
|
|---|
| 154 | /* Approximately update the universal eccentric anomaly. */
|
|---|
| 155 | PSI = PSI+(date-T)*PAL__GCON/R0;
|
|---|
| 156 |
|
|---|
| 157 | /* Time from reference epoch to date (in Canonical Days: a canonical
|
|---|
| 158 | * day is 58.1324409... days, defined as 1/PAL__GCON). */
|
|---|
| 159 | DT = (date-T0)*PAL__GCON;
|
|---|
| 160 |
|
|---|
| 161 | /* Refine the universal eccentric anomaly, psi. */
|
|---|
| 162 | NIT = 1;
|
|---|
| 163 | W = 1.0;
|
|---|
| 164 | TOL = 0.0;
|
|---|
| 165 | while (fabs(W) >= TOL) {
|
|---|
| 166 |
|
|---|
| 167 | /* Form half angles until BETA small enough. */
|
|---|
| 168 | N = 0;
|
|---|
| 169 | PSJ = PSI;
|
|---|
| 170 | PSJ2 = PSJ*PSJ;
|
|---|
| 171 | BETA = ALPHA*PSJ2;
|
|---|
| 172 | while (fabs(BETA) > 0.7) {
|
|---|
| 173 | N = N+1;
|
|---|
| 174 | BETA = BETA/4.0;
|
|---|
| 175 | PSJ = PSJ/2.0;
|
|---|
| 176 | PSJ2 = PSJ2/4.0;
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 | /* Calculate Universal Variables S0,S1,S2,S3 by nested series. */
|
|---|
| 180 | S3 = PSJ*PSJ2*((((((BETA/210.0+1.0)
|
|---|
| 181 | *BETA/156.0+1.0)
|
|---|
| 182 | *BETA/110.0+1.0)
|
|---|
| 183 | *BETA/72.0+1.0)
|
|---|
| 184 | *BETA/42.0+1.0)
|
|---|
| 185 | *BETA/20.0+1.0)/6.0;
|
|---|
| 186 | S2 = PSJ2*((((((BETA/182.0+1.0)
|
|---|
| 187 | *BETA/132.0+1.0)
|
|---|
| 188 | *BETA/90.0+1.0)
|
|---|
| 189 | *BETA/56.0+1.0)
|
|---|
| 190 | *BETA/30.0+1.0)
|
|---|
| 191 | *BETA/12.0+1.0)/2.0;
|
|---|
| 192 | S1 = PSJ+ALPHA*S3;
|
|---|
| 193 | S0 = 1.0+ALPHA*S2;
|
|---|
| 194 |
|
|---|
| 195 | /* Undo the angle-halving. */
|
|---|
| 196 | TOL = TEST;
|
|---|
| 197 | while (N > 0) {
|
|---|
| 198 | S3 = 2.0*(S0*S3+PSJ*S2);
|
|---|
| 199 | S2 = 2.0*S1*S1;
|
|---|
| 200 | S1 = 2.0*S0*S1;
|
|---|
| 201 | S0 = 2.0*S0*S0-1.0;
|
|---|
| 202 | PSJ = PSJ+PSJ;
|
|---|
| 203 | TOL += TOL;
|
|---|
| 204 | N--;
|
|---|
| 205 | }
|
|---|
| 206 |
|
|---|
| 207 | /* Values of F and F' corresponding to the current value of psi. */
|
|---|
| 208 | FF = R0*S1+SIGMA0*S2+CM*S3-DT;
|
|---|
| 209 | R = R0*S0+SIGMA0*S1+CM*S2;
|
|---|
| 210 |
|
|---|
| 211 | /* If first iteration, create dummy "last F". */
|
|---|
| 212 | if ( NIT == 1) FLAST = FF;
|
|---|
| 213 |
|
|---|
| 214 | /* Check for sign change. */
|
|---|
| 215 | if ( FF*FLAST < 0.0 ) {
|
|---|
| 216 |
|
|---|
| 217 | /* Sign change: get psi adjustment using secant method. */
|
|---|
| 218 | W = FF*(PLAST-PSI)/(FLAST-FF);
|
|---|
| 219 | } else {
|
|---|
| 220 |
|
|---|
| 221 | /* No sign change: use Newton-Raphson method instead. */
|
|---|
| 222 | if (R == 0.0) {
|
|---|
| 223 | /* Null radius vector */
|
|---|
| 224 | *jstat = -1;
|
|---|
| 225 | return;
|
|---|
| 226 | }
|
|---|
| 227 | W = FF/R;
|
|---|
| 228 | }
|
|---|
| 229 |
|
|---|
| 230 | /* Save the last psi and F values. */
|
|---|
| 231 | PLAST = PSI;
|
|---|
| 232 | FLAST = FF;
|
|---|
| 233 |
|
|---|
| 234 | /* Apply the Newton-Raphson or secant adjustment to psi. */
|
|---|
| 235 | PSI = PSI-W;
|
|---|
| 236 |
|
|---|
| 237 | /* Next iteration, unless too many already. */
|
|---|
| 238 | if (NIT > NITMAX) {
|
|---|
| 239 | *jstat = -2; /* Failed to converge */
|
|---|
| 240 | return;
|
|---|
| 241 | }
|
|---|
| 242 | NIT++;
|
|---|
| 243 | }
|
|---|
| 244 |
|
|---|
| 245 | /* Project the position and velocity vectors (scaling velocity to AU/s). */
|
|---|
| 246 | W = CM*S2;
|
|---|
| 247 | F = 1.0-W/R0;
|
|---|
| 248 | G = DT-CM*S3;
|
|---|
| 249 | FD = -CM*S1/(R0*R);
|
|---|
| 250 | GD = 1.0-W/R;
|
|---|
| 251 | for (I=0; I<3; I++) {
|
|---|
| 252 | pv[I] = P0[I]*F+V0[I]*G;
|
|---|
| 253 | pv[I+3] = CD2S*(P0[I]*FD+V0[I]*GD);
|
|---|
| 254 | }
|
|---|
| 255 |
|
|---|
| 256 | /* Update the parameters to allow speedy prediction of PSI next time. */
|
|---|
| 257 | u[11] = date;
|
|---|
| 258 | u[12] = PSI;
|
|---|
| 259 |
|
|---|
| 260 | /* OK exit. */
|
|---|
| 261 | *jstat = 0;
|
|---|
| 262 | return;
|
|---|
| 263 | }
|
|---|
