source: branches/FACT++_lidctrl_new_eth/pal/palPertel.c@ 19814

Last change on this file since 19814 was 18347, checked in by tbretz, 9 years ago
File size: 7.5 KB
Line 
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
171void 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}
Note: See TracBrowser for help on using the repository browser.