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
palPertel(int jform,double date0,double date1,double epoch0,double orbi0,double anode0,double perih0,double aorq0,double e0,double am0,double * epoch1,double * orbi1,double * anode1,double * perih1,double * aorq1,double * e1,double * am1,int * jstat)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 }
213