| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palPertel
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Update elements by applying planetary perturbations
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palPertel (int jform, double date0, double date1,
|
|---|
| 17 | * double epoch0, double orbi0, double anode0,
|
|---|
| 18 | * double perih0, double aorq0, double e0, double am0,
|
|---|
| 19 | * double *epoch1, double *orbi1, double *anode1,
|
|---|
| 20 | * double *perih1, double *aorq1, double *e1, double *am1,
|
|---|
| 21 | * int *jstat );
|
|---|
| 22 |
|
|---|
| 23 | * Arguments:
|
|---|
| 24 | * jform = int (Given)
|
|---|
| 25 | * Element set actually returned (1-3; Note 6)
|
|---|
| 26 | * date0 = double (Given)
|
|---|
| 27 | * Date of osculation (TT MJD) for the given elements.
|
|---|
| 28 | * date1 = double (Given)
|
|---|
| 29 | * Date of osculation (TT MJD) for the updated elements.
|
|---|
| 30 | * epoch0 = double (Given)
|
|---|
| 31 | * Epoch of elements (TT MJD)
|
|---|
| 32 | * orbi0 = double (Given)
|
|---|
| 33 | * inclination (radians)
|
|---|
| 34 | * anode0 = double (Given)
|
|---|
| 35 | * longitude of the ascending node (radians)
|
|---|
| 36 | * perih0 = double (Given)
|
|---|
| 37 | * longitude or argument of perihelion (radians)
|
|---|
| 38 | * aorq0 = double (Given)
|
|---|
| 39 | * mean distance or perihelion distance (AU)
|
|---|
| 40 | * e0 = double (Given)
|
|---|
| 41 | * eccentricity
|
|---|
| 42 | * am0 = double (Given)
|
|---|
| 43 | * mean anomaly (radians, JFORM=2 only)
|
|---|
| 44 | * epoch1 = double * (Returned)
|
|---|
| 45 | * Epoch of elements (TT MJD)
|
|---|
| 46 | * orbi1 = double * (Returned)
|
|---|
| 47 | * inclination (radians)
|
|---|
| 48 | * anode1 = double * (Returned)
|
|---|
| 49 | * longitude of the ascending node (radians)
|
|---|
| 50 | * perih1 = double * (Returned)
|
|---|
| 51 | * longitude or argument of perihelion (radians)
|
|---|
| 52 | * aorq1 = double * (Returned)
|
|---|
| 53 | * mean distance or perihelion distance (AU)
|
|---|
| 54 | * e1 = double * (Returned)
|
|---|
| 55 | * eccentricity
|
|---|
| 56 | * am1 = double * (Returned)
|
|---|
| 57 | * mean anomaly (radians, JFORM=2 only)
|
|---|
| 58 | * jstat = int * (Returned)
|
|---|
| 59 | * status:
|
|---|
| 60 | * - +102 = warning, distant epoch
|
|---|
| 61 | * - +101 = warning, large timespan ( > 100 years)
|
|---|
| 62 | * - +1 to +10 = coincident with planet (Note 6)
|
|---|
| 63 | * - 0 = OK
|
|---|
| 64 | * - -1 = illegal JFORM
|
|---|
| 65 | * - -2 = illegal E0
|
|---|
| 66 | * - -3 = illegal AORQ0
|
|---|
| 67 | * - -4 = internal error
|
|---|
| 68 | * - -5 = numerical error
|
|---|
| 69 |
|
|---|
| 70 | * Description:
|
|---|
| 71 | * Update the osculating orbital elements of an asteroid or comet by
|
|---|
| 72 | * applying planetary perturbations.
|
|---|
| 73 |
|
|---|
| 74 | * Authors:
|
|---|
| 75 | * PTW: Pat Wallace (STFC)
|
|---|
| 76 | * TIMJ: Tim Jenness (JAC, Hawaii)
|
|---|
| 77 | * {enter_new_authors_here}
|
|---|
| 78 |
|
|---|
| 79 | * Notes:
|
|---|
| 80 | * - Two different element-format options are available:
|
|---|
| 81 | *
|
|---|
| 82 | * Option JFORM=2, suitable for minor planets:
|
|---|
| 83 | *
|
|---|
| 84 | * EPOCH = epoch of elements (TT MJD)
|
|---|
| 85 | * ORBI = inclination i (radians)
|
|---|
| 86 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 87 | * PERIH = argument of perihelion, little omega (radians)
|
|---|
| 88 | * AORQ = mean distance, a (AU)
|
|---|
| 89 | * E = eccentricity, e
|
|---|
| 90 | * AM = mean anomaly M (radians)
|
|---|
| 91 | *
|
|---|
| 92 | * Option JFORM=3, suitable for comets:
|
|---|
| 93 | *
|
|---|
| 94 | * EPOCH = epoch of perihelion (TT MJD)
|
|---|
| 95 | * ORBI = inclination i (radians)
|
|---|
| 96 | * ANODE = longitude of the ascending node, big omega (radians)
|
|---|
| 97 | * PERIH = argument of perihelion, little omega (radians)
|
|---|
| 98 | * AORQ = perihelion distance, q (AU)
|
|---|
| 99 | * E = eccentricity, e
|
|---|
| 100 | *
|
|---|
| 101 | * - DATE0, DATE1, EPOCH0 and EPOCH1 are all instants of time in
|
|---|
| 102 | * the TT timescale (formerly Ephemeris Time, ET), expressed
|
|---|
| 103 | * as Modified Julian Dates (JD-2400000.5).
|
|---|
| 104 | *
|
|---|
| 105 | * DATE0 is the instant at which the given (i.e. unperturbed)
|
|---|
| 106 | * osculating elements are correct.
|
|---|
| 107 | *
|
|---|
| 108 | * DATE1 is the specified instant at which the updated osculating
|
|---|
| 109 | * elements are correct.
|
|---|
| 110 | *
|
|---|
| 111 | * EPOCH0 and EPOCH1 will be the same as DATE0 and DATE1
|
|---|
| 112 | * (respectively) for the JFORM=2 case, normally used for minor
|
|---|
| 113 | * planets. For the JFORM=3 case, the two epochs will refer to
|
|---|
| 114 | * perihelion passage and so will not, in general, be the same as
|
|---|
| 115 | * DATE0 and/or DATE1 though they may be similar to one another.
|
|---|
| 116 | * - The elements are with respect to the J2000 ecliptic and equinox.
|
|---|
| 117 | * - Unused elements (AM0 and AM1 for JFORM=3) are not accessed.
|
|---|
| 118 | * - See the palPertue routine for details of the algorithm used.
|
|---|
| 119 | * - This routine is not intended to be used for major planets, which
|
|---|
| 120 | * is why JFORM=1 is not available and why there is no opportunity
|
|---|
| 121 | * to specify either the longitude of perihelion or the daily
|
|---|
| 122 | * motion. However, if JFORM=2 elements are somehow obtained for a
|
|---|
| 123 | * major planet and supplied to the routine, sensible results will,
|
|---|
| 124 | * in fact, be produced. This happens because the palPertue routine
|
|---|
| 125 | * that is called to perform the calculations checks the separation
|
|---|
| 126 | * between the body and each of the planets and interprets a
|
|---|
| 127 | * suspiciously small value (0.001 AU) as an attempt to apply it to
|
|---|
| 128 | * the planet concerned. If this condition is detected, the
|
|---|
| 129 | * contribution from that planet is ignored, and the status is set to
|
|---|
| 130 | * the planet number (1-10 = Mercury, Venus, EMB, Mars, Jupiter,
|
|---|
| 131 | * Saturn, Uranus, Neptune, Earth, Moon) as a warning.
|
|---|
| 132 | *
|
|---|
| 133 | * See Also:
|
|---|
| 134 | * - Sterne, Theodore E., "An Introduction to Celestial Mechanics",
|
|---|
| 135 | * Interscience Publishers Inc., 1960. Section 6.7, p199.
|
|---|
| 136 |
|
|---|
| 137 | * History:
|
|---|
| 138 | * 2012-03-12 (TIMJ):
|
|---|
| 139 | * Initial version direct conversion of SLA/F.
|
|---|
| 140 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 141 | * {enter_further_changes_here}
|
|---|
| 142 |
|
|---|
| 143 | * Copyright:
|
|---|
| 144 | * Copyright (C) 2004 Patrick T. Wallace
|
|---|
| 145 | * Copyright (C) 2012 Science and Technology Facilities Council.
|
|---|
| 146 | * All Rights Reserved.
|
|---|
| 147 |
|
|---|
| 148 | * Licence:
|
|---|
| 149 | * This program is free software; you can redistribute it and/or
|
|---|
| 150 | * modify it under the terms of the GNU General Public License as
|
|---|
| 151 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 152 | * the License, or (at your option) any later version.
|
|---|
| 153 | *
|
|---|
| 154 | * This program is distributed in the hope that it will be
|
|---|
| 155 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 156 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 157 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 158 | *
|
|---|
| 159 | * You should have received a copy of the GNU General Public License
|
|---|
| 160 | * along with this program; if not, write to the Free Software
|
|---|
| 161 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 162 | * MA 02110-1301, USA.
|
|---|
| 163 |
|
|---|
| 164 | * Bugs:
|
|---|
| 165 | * {note_any_bugs_here}
|
|---|
| 166 | *-
|
|---|
| 167 | */
|
|---|
| 168 |
|
|---|
| 169 | #include "pal.h"
|
|---|
| 170 |
|
|---|
| 171 | void palPertel (int jform, double date0, double date1,
|
|---|
| 172 | double epoch0, double orbi0, double anode0,
|
|---|
| 173 | double perih0, double aorq0, double e0, double am0,
|
|---|
| 174 | double *epoch1, double *orbi1, double *anode1,
|
|---|
| 175 | double *perih1, double *aorq1, double *e1, double *am1,
|
|---|
| 176 | int *jstat ) {
|
|---|
| 177 |
|
|---|
| 178 | double u[13], dm;
|
|---|
| 179 | int j, jf;
|
|---|
| 180 |
|
|---|
| 181 | /* Check that the elements are either minor-planet or comet format. */
|
|---|
| 182 | if (jform < 2 || jform > 3) {
|
|---|
| 183 | *jstat = -1;
|
|---|
| 184 | return;
|
|---|
| 185 | } else {
|
|---|
| 186 |
|
|---|
| 187 | /* Provisionally set the status to OK. */
|
|---|
| 188 | *jstat = 0;
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 | /* Transform the elements from conventional to universal form. */
|
|---|
| 192 | palEl2ue(date0,jform,epoch0,orbi0,anode0,perih0,
|
|---|
| 193 | aorq0,e0,am0,0.0,u,&j);
|
|---|
| 194 | if (j != 0) {
|
|---|
| 195 | *jstat = j;
|
|---|
| 196 | return;
|
|---|
| 197 | }
|
|---|
| 198 |
|
|---|
| 199 | /* Update the universal elements. */
|
|---|
| 200 | palPertue(date1,u,&j);
|
|---|
| 201 | if (j > 0) {
|
|---|
| 202 | *jstat = j;
|
|---|
| 203 | } else if (j < 0) {
|
|---|
| 204 | *jstat = -5;
|
|---|
| 205 | return;
|
|---|
| 206 | }
|
|---|
| 207 |
|
|---|
| 208 | /* Transform from universal to conventional elements. */
|
|---|
| 209 | palUe2el(u, jform, &jf, epoch1, orbi1, anode1, perih1,
|
|---|
| 210 | aorq1, e1, am1, &dm, &j);
|
|---|
| 211 | if (jf != jform || j != 0) *jstat = -5;
|
|---|
| 212 | }
|
|---|
