| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palPv2el
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Position velocity to heliocentirc osculating elements
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palPv2el ( const double pv[6], double date, double pmass, int jformr,
|
|---|
| 17 | * int *jform, double *epoch, double *orbinc,
|
|---|
| 18 | * double *anode, double *perih, double *aorq, double *e,
|
|---|
| 19 | * double *aorl, double *dm, int *jstat );
|
|---|
| 20 |
|
|---|
| 21 | * Arguments:
|
|---|
| 22 | * pv = const double [6] (Given)
|
|---|
| 23 | * Heliocentric x,y,z,xdot,ydot,zdot of date,
|
|---|
| 24 | * J2000 equatorial triad (AU,AU/s; Note 1)
|
|---|
| 25 | * date = double (Given)
|
|---|
| 26 | * Date (TT Modified Julian Date = JD-2400000.5)
|
|---|
| 27 | * pmass = double (Given)
|
|---|
| 28 | * Mass of the planet (Sun=1; Note 2)
|
|---|
| 29 | * jformr = int (Given)
|
|---|
| 30 | * Requested element set (1-3; Note 3)
|
|---|
| 31 | * jform = int * (Returned)
|
|---|
| 32 | * Element set actually returned (1-3; Note 4)
|
|---|
| 33 | * epoch = double * (Returned)
|
|---|
| 34 | * Epoch of elements (TT MJD)
|
|---|
| 35 | * orbinc = double * (Returned)
|
|---|
| 36 | * inclination (radians)
|
|---|
| 37 | * anode = double * (Returned)
|
|---|
| 38 | * longitude of the ascending node (radians)
|
|---|
| 39 | * perih = double * (Returned)
|
|---|
| 40 | * longitude or argument of perihelion (radians)
|
|---|
| 41 | * aorq = double * (Returned)
|
|---|
| 42 | * mean distance or perihelion distance (AU)
|
|---|
| 43 | * e = double * (Returned)
|
|---|
| 44 | * eccentricity
|
|---|
| 45 | * aorl = double * (Returned)
|
|---|
| 46 | * mean anomaly or longitude (radians, JFORM=1,2 only)
|
|---|
| 47 | * dm = double * (Returned)
|
|---|
| 48 | * daily motion (radians, JFORM=1 only)
|
|---|
| 49 | * jstat = int * (Returned)
|
|---|
| 50 | * status: 0 = OK
|
|---|
| 51 | * - -1 = illegal PMASS
|
|---|
| 52 | * - -2 = illegal JFORMR
|
|---|
| 53 | * - -3 = position/velocity out of range
|
|---|
| 54 |
|
|---|
| 55 | * Description:
|
|---|
| 56 | * Heliocentric osculating elements obtained from instantaneous position
|
|---|
| 57 | * and velocity.
|
|---|
| 58 |
|
|---|
| 59 | * Authors:
|
|---|
| 60 | * PTW: Pat Wallace (STFC)
|
|---|
| 61 | * TIMJ: Tim Jenness (JAC, Hawaii)
|
|---|
| 62 | * {enter_new_authors_here}
|
|---|
| 63 |
|
|---|
| 64 | * Notes:
|
|---|
| 65 | * - The PV 6-vector is with respect to the mean equator and equinox of
|
|---|
| 66 | * epoch J2000. The orbital elements produced are with respect to
|
|---|
| 67 | * the J2000 ecliptic and mean equinox.
|
|---|
| 68 | * - The mass, PMASS, is important only for the larger planets. For
|
|---|
| 69 | * most purposes (e.g. asteroids) use 0D0. Values less than zero
|
|---|
| 70 | * are illegal.
|
|---|
| 71 | * - Three different element-format options are supported:
|
|---|
| 72 | *
|
|---|
| 73 | * Option JFORM=1, suitable for the major planets:
|
|---|
| 74 | *
|
|---|
| 75 | * EPOCH = epoch of elements (TT MJD)
|
|---|
| 76 | * ORBINC = inclination i (radians)
|
|---|
| 77 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 78 | * PERIH = longitude of perihelion, curly pi (radians)
|
|---|
| 79 | * AORQ = mean distance, a (AU)
|
|---|
| 80 | * E = eccentricity, e
|
|---|
| 81 | * AORL = mean longitude L (radians)
|
|---|
| 82 | * DM = daily motion (radians)
|
|---|
| 83 | *
|
|---|
| 84 | * Option JFORM=2, suitable for minor planets:
|
|---|
| 85 | *
|
|---|
| 86 | * EPOCH = epoch of elements (TT MJD)
|
|---|
| 87 | * ORBINC = inclination i (radians)
|
|---|
| 88 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 89 | * PERIH = argument of perihelion, little omega (radians)
|
|---|
| 90 | * AORQ = mean distance, a (AU)
|
|---|
| 91 | * E = eccentricity, e
|
|---|
| 92 | * AORL = mean anomaly M (radians)
|
|---|
| 93 | *
|
|---|
| 94 | * Option JFORM=3, suitable for comets:
|
|---|
| 95 | *
|
|---|
| 96 | * EPOCH = epoch of perihelion (TT MJD)
|
|---|
| 97 | * ORBINC = inclination i (radians)
|
|---|
| 98 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 99 | * PERIH = argument of perihelion, little omega (radians)
|
|---|
| 100 | * AORQ = perihelion distance, q (AU)
|
|---|
| 101 | * E = eccentricity, e
|
|---|
| 102 | *
|
|---|
| 103 | * - It may not be possible to generate elements in the form
|
|---|
| 104 | * requested through JFORMR. The caller is notified of the form
|
|---|
| 105 | * of elements actually returned by means of the JFORM argument:
|
|---|
| 106 |
|
|---|
| 107 | * JFORMR JFORM meaning
|
|---|
| 108 | *
|
|---|
| 109 | * 1 1 OK - elements are in the requested format
|
|---|
| 110 | * 1 2 never happens
|
|---|
| 111 | * 1 3 orbit not elliptical
|
|---|
| 112 | *
|
|---|
| 113 | * 2 1 never happens
|
|---|
| 114 | * 2 2 OK - elements are in the requested format
|
|---|
| 115 | * 2 3 orbit not elliptical
|
|---|
| 116 | *
|
|---|
| 117 | * 3 1 never happens
|
|---|
| 118 | * 3 2 never happens
|
|---|
| 119 | * 3 3 OK - elements are in the requested format
|
|---|
| 120 | *
|
|---|
| 121 | * - The arguments returned for each value of JFORM (cf Note 5: JFORM
|
|---|
| 122 | * may not be the same as JFORMR) are as follows:
|
|---|
| 123 | *
|
|---|
| 124 | * JFORM 1 2 3
|
|---|
| 125 | * EPOCH t0 t0 T
|
|---|
| 126 | * ORBINC i i i
|
|---|
| 127 | * ANODE Omega Omega Omega
|
|---|
| 128 | * PERIH curly pi omega omega
|
|---|
| 129 | * AORQ a a q
|
|---|
| 130 | * E e e e
|
|---|
| 131 | * AORL L M -
|
|---|
| 132 | * DM n - -
|
|---|
| 133 | *
|
|---|
| 134 | * where:
|
|---|
| 135 | *
|
|---|
| 136 | * t0 is the epoch of the elements (MJD, TT)
|
|---|
| 137 | * T " epoch of perihelion (MJD, TT)
|
|---|
| 138 | * i " inclination (radians)
|
|---|
| 139 | * Omega " longitude of the ascending node (radians)
|
|---|
| 140 | * curly pi " longitude of perihelion (radians)
|
|---|
| 141 | * omega " argument of perihelion (radians)
|
|---|
| 142 | * a " mean distance (AU)
|
|---|
| 143 | * q " perihelion distance (AU)
|
|---|
| 144 | * e " eccentricity
|
|---|
| 145 | * L " longitude (radians, 0-2pi)
|
|---|
| 146 | * M " mean anomaly (radians, 0-2pi)
|
|---|
| 147 | * n " daily motion (radians)
|
|---|
| 148 | * - means no value is set
|
|---|
| 149 | *
|
|---|
| 150 | * - At very small inclinations, the longitude of the ascending node
|
|---|
| 151 | * ANODE becomes indeterminate and under some circumstances may be
|
|---|
| 152 | * set arbitrarily to zero. Similarly, if the orbit is close to
|
|---|
| 153 | * circular, the true anomaly becomes indeterminate and under some
|
|---|
| 154 | * circumstances may be set arbitrarily to zero. In such cases,
|
|---|
| 155 | * the other elements are automatically adjusted to compensate,
|
|---|
| 156 | * and so the elements remain a valid description of the orbit.
|
|---|
| 157 | * - The osculating epoch for the returned elements is the argument
|
|---|
| 158 | * DATE.
|
|---|
| 159 | *
|
|---|
| 160 | * - Reference: Sterne, Theodore E., "An Introduction to Celestial
|
|---|
| 161 | * Mechanics", Interscience Publishers, 1960
|
|---|
| 162 |
|
|---|
| 163 | * History:
|
|---|
| 164 | * 2012-03-09 (TIMJ):
|
|---|
| 165 | * Initial version converted from SLA/F.
|
|---|
| 166 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 167 | * {enter_further_changes_here}
|
|---|
| 168 |
|
|---|
| 169 | * Copyright:
|
|---|
| 170 | * Copyright (C) 2005 Patrick T. Wallace
|
|---|
| 171 | * Copyright (C) 2012 Science and Technology Facilities Council.
|
|---|
| 172 | * All Rights Reserved.
|
|---|
| 173 |
|
|---|
| 174 | * Licence:
|
|---|
| 175 | * This program is free software; you can redistribute it and/or
|
|---|
| 176 | * modify it under the terms of the GNU General Public License as
|
|---|
| 177 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 178 | * the License, or (at your option) any later version.
|
|---|
| 179 | *
|
|---|
| 180 | * This program is distributed in the hope that it will be
|
|---|
| 181 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 182 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 183 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 184 | *
|
|---|
| 185 | * You should have received a copy of the GNU General Public License
|
|---|
| 186 | * along with this program; if not, write to the Free Software
|
|---|
| 187 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 188 | * MA 02110-1301, USA.
|
|---|
| 189 |
|
|---|
| 190 | * Bugs:
|
|---|
| 191 | * {note_any_bugs_here}
|
|---|
| 192 | *-
|
|---|
| 193 | */
|
|---|
| 194 |
|
|---|
| 195 | #include <math.h>
|
|---|
| 196 |
|
|---|
| 197 | #include "pal1sofa.h"
|
|---|
| 198 | #include "pal.h"
|
|---|
| 199 | #include "palmac.h"
|
|---|
| 200 |
|
|---|
| 201 | void palPv2el ( const double pv[6], double date, double pmass, int jformr,
|
|---|
| 202 | int *jform, double *epoch, double *orbinc,
|
|---|
| 203 | double *anode, double *perih, double *aorq, double *e,
|
|---|
| 204 | double *aorl, double *dm, int *jstat ) {
|
|---|
| 205 |
|
|---|
| 206 | /* Sin and cos of J2000 mean obliquity (IAU 1976) */
|
|---|
| 207 | const double SE = 0.3977771559319137;
|
|---|
| 208 | const double CE = 0.9174820620691818;
|
|---|
| 209 |
|
|---|
| 210 | /* Minimum allowed distance (AU) and speed (AU/day) */
|
|---|
| 211 | const double RMIN = 1e-3;
|
|---|
| 212 | const double VMIN = 1e-8;
|
|---|
| 213 |
|
|---|
| 214 | /* How close to unity the eccentricity has to be to call it a parabola */
|
|---|
| 215 | const double PARAB = 1.0e-8;
|
|---|
| 216 |
|
|---|
| 217 | double X,Y,Z,XD,YD,ZD,R,V2,V,RDV,GMU,HX,HY,HZ,
|
|---|
| 218 | HX2PY2,H2,H,OI,BIGOM,AR,ECC,S,C,AT,U,OM,
|
|---|
| 219 | GAR3,EM1,EP1,HAT,SHAT,CHAT,AE,AM,DN,PL,
|
|---|
| 220 | EL,Q,TP,THAT,THHF,F;
|
|---|
| 221 |
|
|---|
| 222 | int JF;
|
|---|
| 223 |
|
|---|
| 224 | /* Validate arguments PMASS and JFORMR.*/
|
|---|
| 225 | if (pmass < 0.0) {
|
|---|
| 226 | *jstat = -1;
|
|---|
| 227 | return;
|
|---|
| 228 | }
|
|---|
| 229 | if (jformr < 1 || jformr > 3) {
|
|---|
| 230 | *jstat = -2;
|
|---|
| 231 | return;
|
|---|
| 232 | }
|
|---|
| 233 |
|
|---|
| 234 | /* Provisionally assume the elements will be in the chosen form. */
|
|---|
| 235 | JF = jformr;
|
|---|
| 236 |
|
|---|
| 237 | /* Rotate the position from equatorial to ecliptic coordinates. */
|
|---|
| 238 | X = pv[0];
|
|---|
| 239 | Y = pv[1]*CE+pv[2]*SE;
|
|---|
| 240 | Z = -pv[1]*SE+pv[2]*CE;
|
|---|
| 241 |
|
|---|
| 242 | /* Rotate the velocity similarly, scaling to AU/day. */
|
|---|
| 243 | XD = PAL__SPD*pv[3];
|
|---|
| 244 | YD = PAL__SPD*(pv[4]*CE+pv[5]*SE);
|
|---|
| 245 | ZD = PAL__SPD*(-pv[4]*SE+pv[5]*CE);
|
|---|
| 246 |
|
|---|
| 247 | /* Distance and speed. */
|
|---|
| 248 | R = sqrt(X*X+Y*Y+Z*Z);
|
|---|
| 249 | V2 = XD*XD+YD*YD+ZD*ZD;
|
|---|
| 250 | V = sqrt(V2);
|
|---|
| 251 |
|
|---|
| 252 | /* Reject unreasonably small values. */
|
|---|
| 253 | if (R < RMIN || V < VMIN) {
|
|---|
| 254 | *jstat = -3;
|
|---|
| 255 | return;
|
|---|
| 256 | }
|
|---|
| 257 |
|
|---|
| 258 | /* R dot V. */
|
|---|
| 259 | RDV = X*XD+Y*YD+Z*ZD;
|
|---|
| 260 |
|
|---|
| 261 | /* Mu. */
|
|---|
| 262 | GMU = (1.0+pmass)*PAL__GCON*PAL__GCON;
|
|---|
| 263 |
|
|---|
| 264 | /* Vector angular momentum per unit reduced mass. */
|
|---|
| 265 | HX = Y*ZD-Z*YD;
|
|---|
| 266 | HY = Z*XD-X*ZD;
|
|---|
| 267 | HZ = X*YD-Y*XD;
|
|---|
| 268 |
|
|---|
| 269 | /* Areal constant. */
|
|---|
| 270 | HX2PY2 = HX*HX+HY*HY;
|
|---|
| 271 | H2 = HX2PY2+HZ*HZ;
|
|---|
| 272 | H = sqrt(H2);
|
|---|
| 273 |
|
|---|
| 274 | /* Inclination. */
|
|---|
| 275 | OI = atan2(sqrt(HX2PY2),HZ);
|
|---|
| 276 |
|
|---|
| 277 | /* Longitude of ascending node. */
|
|---|
| 278 | if (HX != 0.0 || HY != 0.0) {
|
|---|
| 279 | BIGOM = atan2(HX,-HY);
|
|---|
| 280 | } else {
|
|---|
| 281 | BIGOM=0.0;
|
|---|
| 282 | }
|
|---|
| 283 |
|
|---|
| 284 | /* Reciprocal of mean distance etc. */
|
|---|
| 285 | AR = 2.0/R-V2/GMU;
|
|---|
| 286 |
|
|---|
| 287 | /* Eccentricity. */
|
|---|
| 288 | ECC = sqrt(DMAX(1.0-AR*H2/GMU,0.0));
|
|---|
| 289 |
|
|---|
| 290 | /* True anomaly. */
|
|---|
| 291 | S = H*RDV;
|
|---|
| 292 | C = H2-R*GMU;
|
|---|
| 293 | if (S != 0.0 || C != 0.0) {
|
|---|
| 294 | AT = atan2(S,C);
|
|---|
| 295 | } else {
|
|---|
| 296 | AT = 0.0;
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 | /* Argument of the latitude. */
|
|---|
| 300 | S = sin(BIGOM);
|
|---|
| 301 | C = cos(BIGOM);
|
|---|
| 302 | U = atan2((-X*S+Y*C)*cos(OI)+Z*sin(OI),X*C+Y*S);
|
|---|
| 303 |
|
|---|
| 304 | /* Argument of perihelion. */
|
|---|
| 305 | OM = U-AT;
|
|---|
| 306 |
|
|---|
| 307 | /* Capture near-parabolic cases. */
|
|---|
| 308 | if (fabs(ECC-1.0) < PARAB) ECC=1.0;
|
|---|
| 309 |
|
|---|
| 310 | /* Comply with JFORMR = 1 or 2 only if orbit is elliptical. */
|
|---|
| 311 | if (ECC > 1.0) JF=3;
|
|---|
| 312 |
|
|---|
| 313 | /* Functions. */
|
|---|
| 314 | GAR3 = GMU*AR*AR*AR;
|
|---|
| 315 | EM1 = ECC-1.0;
|
|---|
| 316 | EP1 = ECC+1.0;
|
|---|
| 317 | HAT = AT/2.0;
|
|---|
| 318 | SHAT = sin(HAT);
|
|---|
| 319 | CHAT = cos(HAT);
|
|---|
| 320 |
|
|---|
| 321 | /* Variable initializations to avoid compiler warnings. */
|
|---|
| 322 | AM = 0.0;
|
|---|
| 323 | DN = 0.0;
|
|---|
| 324 | PL = 0.0;
|
|---|
| 325 | EL = 0.0;
|
|---|
| 326 | Q = 0.0;
|
|---|
| 327 | TP = 0.0;
|
|---|
| 328 |
|
|---|
| 329 | /* Ellipse? */
|
|---|
| 330 | if (ECC < 1.0 ) {
|
|---|
| 331 |
|
|---|
| 332 | /* Eccentric anomaly. */
|
|---|
| 333 | AE = 2.0*atan2(sqrt(-EM1)*SHAT,sqrt(EP1)*CHAT);
|
|---|
| 334 |
|
|---|
| 335 | /* Mean anomaly. */
|
|---|
| 336 | AM = AE-ECC*sin(AE);
|
|---|
| 337 |
|
|---|
| 338 | /* Daily motion. */
|
|---|
| 339 | DN = sqrt(GAR3);
|
|---|
| 340 | }
|
|---|
| 341 |
|
|---|
| 342 | /* "Major planet" element set? */
|
|---|
| 343 | if (JF == 1) {
|
|---|
| 344 |
|
|---|
| 345 | /* Longitude of perihelion. */
|
|---|
| 346 | PL = BIGOM+OM;
|
|---|
| 347 |
|
|---|
| 348 | /* Longitude at epoch. */
|
|---|
| 349 | EL = PL+AM;
|
|---|
| 350 | }
|
|---|
| 351 |
|
|---|
| 352 | /* "Comet" element set? */
|
|---|
| 353 | if (JF == 3) {
|
|---|
| 354 |
|
|---|
| 355 | /* Perihelion distance. */
|
|---|
| 356 | Q = H2/(GMU*EP1);
|
|---|
| 357 |
|
|---|
| 358 | /* Ellipse, parabola, hyperbola? */
|
|---|
| 359 | if (ECC < 1.0) {
|
|---|
| 360 |
|
|---|
| 361 | /* Ellipse: epoch of perihelion. */
|
|---|
| 362 | TP = date-AM/DN;
|
|---|
| 363 |
|
|---|
| 364 | } else {
|
|---|
| 365 |
|
|---|
| 366 | /* Parabola or hyperbola: evaluate tan ( ( true anomaly ) / 2 ) */
|
|---|
| 367 | THAT = SHAT/CHAT;
|
|---|
| 368 | if (ECC == 1.0) {
|
|---|
| 369 |
|
|---|
| 370 | /* Parabola: epoch of perihelion. */
|
|---|
| 371 | TP = date-THAT*(1.0+THAT*THAT/3.0)*H*H2/(2.0*GMU*GMU);
|
|---|
| 372 |
|
|---|
| 373 | } else {
|
|---|
| 374 |
|
|---|
| 375 | /* Hyperbola: epoch of perihelion. */
|
|---|
| 376 | THHF = sqrt(EM1/EP1)*THAT;
|
|---|
| 377 | F = log(1.0+THHF)-log(1.0-THHF);
|
|---|
| 378 | TP = date-(ECC*sinh(F)-F)/sqrt(-GAR3);
|
|---|
| 379 | }
|
|---|
| 380 | }
|
|---|
| 381 | }
|
|---|
| 382 |
|
|---|
| 383 | /* Return the appropriate set of elements. */
|
|---|
| 384 | *jform = JF;
|
|---|
| 385 | *orbinc = OI;
|
|---|
| 386 | *anode = eraAnp(BIGOM);
|
|---|
| 387 | *e = ECC;
|
|---|
| 388 | if (JF == 1) {
|
|---|
| 389 | *perih = eraAnp(PL);
|
|---|
| 390 | *aorl = eraAnp(EL);
|
|---|
| 391 | *dm = DN;
|
|---|
| 392 | } else {
|
|---|
| 393 | *perih = eraAnp(OM);
|
|---|
| 394 | if (JF == 2) *aorl = eraAnp(AM);
|
|---|
| 395 | }
|
|---|
| 396 | if (JF != 3) {
|
|---|
| 397 | *epoch = date;
|
|---|
| 398 | *aorq = 1.0/AR;
|
|---|
| 399 | } else {
|
|---|
| 400 | *epoch = TP;
|
|---|
| 401 | *aorq = Q;
|
|---|
| 402 | }
|
|---|
| 403 | *jstat = 0;
|
|---|
| 404 |
|
|---|
| 405 | }
|
|---|
