1 /* conics.f -- translated by f2c (version 19980913).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* $Procedure CONICS ( Determine state from conic elements ) */
conics_(doublereal * elts,doublereal * et,doublereal * state)9 /* Subroutine */ int conics_(doublereal *elts, doublereal *et, doublereal *
10 state)
11 {
12 /* System generated locals */
13 doublereal d__1;
14
15 /* Builtin functions */
16 double cos(doublereal), sin(doublereal), sqrt(doublereal), d_mod(
17 doublereal *, doublereal *);
18
19 /* Local variables */
20 doublereal cnci, argp, snci, cosi, sini, cosn, sinn;
21 extern /* Subroutine */ int vscl_(doublereal *, doublereal *, doublereal *
22 );
23 doublereal cosw, sinw, n, v;
24 extern /* Subroutine */ int chkin_(char *, ftnlen);
25 doublereal lnode;
26 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen);
27 doublereal m0;
28 extern doublereal twopi_(void);
29 doublereal t0;
30 extern /* Subroutine */ int prop2b_(doublereal *, doublereal *,
31 doublereal *, doublereal *);
32 doublereal dt, rp, mu, basisp[3], period, basisq[3];
33 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
34 ftnlen);
35 doublereal pstate[6], ainvrs;
36 extern /* Subroutine */ int setmsg_(char *, ftnlen);
37 extern logical return_(void);
38 doublereal ecc, inc;
39
40 /* $ Abstract */
41
42 /* Determine the state (position, velocity) of an orbiting body */
43 /* from a set of elliptic, hyperbolic, or parabolic orbital */
44 /* elements. */
45
46 /* $ Disclaimer */
47
48 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
49 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
50 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
51 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
52 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
53 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
54 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
55 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
56 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
57 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
58
59 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
60 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
61 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
62 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
63 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
64 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
65
66 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
67 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
68 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
69 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
70
71 /* $ Required_Reading */
72
73 /* None. */
74
75 /* $ Keywords */
76
77 /* CONIC */
78 /* EPHEMERIS */
79
80 /* $ Declarations */
81 /* $ Brief_I/O */
82
83 /* VARIABLE I/O DESCRIPTION */
84 /* -------- --- -------------------------------------------------- */
85 /* ELTS I Conic elements. */
86 /* ET I Input time. */
87 /* STATE O State of orbiting body at ET. */
88
89 /* $ Detailed_Input */
90
91 /* ELTS are conic elements describing the orbit of a body */
92 /* around a primary. The elements are, in order: */
93
94 /* RP Perifocal distance. */
95 /* ECC Eccentricity. */
96 /* INC Inclination. */
97 /* LNODE Longitude of the ascending node. */
98 /* ARGP Argument of periapse. */
99 /* M0 Mean anomaly at epoch. */
100 /* T0 Epoch. */
101 /* MU Gravitational parameter. */
102
103 /* Units are km, rad, rad/sec, km**3/sec**2. The epoch */
104 /* is given in ephemeris seconds past J2000. The same */
105 /* elements are used to describe all three types */
106 /* (elliptic, hyperbolic, and parabolic) of conic orbit. */
107
108 /* ET is the time at which the state of the orbiting body */
109 /* is to be determined, in ephemeris seconds J2000. */
110
111 /* $ Detailed_Output */
112
113 /* STATE is the state (position and velocity) of the body at */
114 /* time ET. Components are x, y, z, dx/dt, dy/dt, dz/dt. */
115
116 /* $ Parameters */
117
118 /* None. */
119
120 /* $ Exceptions */
121
122 /* 1) If the eccentricity supplied is less than 0, the error */
123 /* 'SPICE(BADECCENTRICITY)' is signalled. */
124
125 /* 2) If a non-positive periapse distance is supplied, the error */
126 /* 'SPICE(BADPERIAPSEVALUE)' is signalled. */
127
128 /* 3) If a non-positive value for the attracting mass is supplied, */
129 /* the error 'SPICE(BADGM)', is signalled. */
130
131 /* 4) Errors such as an out of bounds value for ET are diagnosed */
132 /* by routines called by this routine. */
133
134 /* $ Files */
135
136 /* None. */
137
138 /* $ Particulars */
139
140 /* None. */
141
142 /* $ Examples */
143
144 /* Let VINIT contain the initial state of a spacecraft relative to */
145 /* the center of a planet at epoch ET, and let GM be the gravitation */
146 /* parameter of the planet. The call */
147
148 /* CALL OSCELT ( VINIT, ET, GM, ELTS ) */
149
150 /* produces a set of osculating elements describing the nominal */
151 /* orbit that the spacecraft would follow in the absence of all */
152 /* other bodies in the solar system and non-gravitational forces */
153 /* on the spacecraft. */
154
155 /* Now let STATE contain the state of the same spacecraft at some */
156 /* other epoch, LATER. The difference between this state and the */
157 /* state predicted by the nominal orbit at the same epoch can be */
158 /* computed as follows. */
159
160 /* CALL CONICS ( ELTS, LATER, NOMINAL ) */
161 /* CALL VSUBG ( NOMINAL, STATE, 6, DIFF ) */
162
163 /* WRITE (*,*) 'Perturbation in x, dx/dt = ', DIFF(1), DIFF(4) */
164 /* WRITE (*,*) ' y, dy/dt = ', DIFF(2), DIFF(5) */
165 /* WRITE (*,*) ' z, dz/dt = ', DIFF(3), DIFF(6) */
166
167 /* $ Restrictions */
168
169 /* None. */
170
171 /* $ Literature_References */
172
173 /* [1] Roger Bate, Fundamentals of Astrodynamics, Dover, 1971. */
174
175 /* $ Author_and_Institution */
176
177 /* I.M. Underwood (JPL) */
178 /* W.L. Taber (JPL) */
179
180 /* $ Version */
181
182 /* - SPICELIB Version 4.0.0, 26-MAR-1998 (WLT) */
183
184 /* There was a coding error in the computation of the mean */
185 /* anomaly in the parabolic case. This problem has been */
186 /* corrected. */
187
188 /* - SPICELIB Version 3.0.1, 15-OCT-1996 (WLT) */
189
190 /* Corrected a typo in the description of the units associated */
191 /* with the input elements. */
192
193 /* - SPICELIB Version 3.0.0, 12-NOV-1992 (WLT) */
194
195 /* The routine was re-written to make use of NAIF's universal */
196 /* variables formulation for state propagation (PROP2B). As */
197 /* a result, several problems were simultaneously corrected. */
198
199 /* A major bug was fixed that caused improper state evaluations */
200 /* for ET's that precede the epoch of the elements in the */
201 /* elliptic case. */
202
203 /* A danger of non-convergence in the solution of Kepler's */
204 /* equation has been eliminated. */
205
206 /* In addition to this reformulation of CONICS checks were */
207 /* installed that ensure the elements supplied are physically */
208 /* meaningful. Eccentricity must be non-negative. The */
209 /* distance at periapse and central mass must be positive. If */
210 /* not errors are signalled. */
211
212 /* - SPICELIB Version 2.0.1, 10-MAR-1992 (WLT) */
213
214 /* Comment section for permuted index source lines was added */
215 /* following the header. */
216
217 /* - SPICELIB Version 2.0.0, 19-APR-1991 (WLT) */
218
219 /* An error in the hyperbolic state generation was corrected. */
220
221 /* - SPICELIB Version 1.0.0, 31-JAN-1990 (IMU) */
222
223 /* -& */
224 /* $ Index_Entries */
225
226 /* state from conic elements */
227
228 /* -& */
229 /* $ Revisions */
230
231 /* - SPICELIB Version 3.0.1, 15-OCT-1996 (WLT) */
232
233 /* Corrected a typo in the description of the units associated */
234 /* with the input elements. */
235
236 /* - SPICELIB Version 3.0.0, 12-NOV-1992 (WLT) */
237
238 /* The routine was re-written to make use of NAIF's universal */
239 /* variables formulation for state propagation (PROP2B). As */
240 /* a result, several problems were simultaneously corrected. */
241
242 /* A major bug was fixed that caused improper state evaluations */
243 /* for ET's that precede the epoch of the elements in the */
244 /* elliptic case. */
245
246 /* A danger of non-convergence in the solution of Kepler's */
247 /* equation has been eliminated. */
248
249 /* In addition to this reformulation of CONICS checks were */
250 /* installed that ensure the elements supplied are physically */
251 /* meaningful. Eccentricity must be non-negative. The */
252 /* distance at periapse and central mass must be positive. If */
253 /* not errors are signalled. */
254
255 /* These changes were prompted by the discovery that the old */
256 /* formulation had a severe bug for elliptic orbits and epochs */
257 /* prior to the epoch of the input elements, and by the discovery */
258 /* that the time of flight routines had problems with convergence. */
259
260 /* - SPICELIB Version 2.0.0, 19-APR-1991 (WLT) */
261
262 /* The original version of the routine had a bug in that */
263 /* it attempted to restrict the hyperbolic anomaly to */
264 /* the interval 0 to 2*PI. This has been fixed. */
265
266 /* - Beta Version 1.0.1, 27-JAN-1989 (IMU) */
267
268 /* Examples section completed. */
269
270 /* -& */
271
272 /* SPICELIB functions */
273
274
275 /* Local variables */
276
277
278 /* The only real work required by this routine is the construction */
279 /* of a preliminary state vector from the input elements. Once this */
280 /* is in hand, we can simply let the routine PROP2B do the real */
281 /* work, free from the instabilities inherent in the classical */
282 /* elements formulation of two-body motion. */
283
284 /* To do this we shall construct a basis of vectors that lie in the */
285 /* plane of the orbit. The first vector P shall point towards the */
286 /* position of the orbiting body at periapse. The second */
287 /* vector Q shall point along the velocity vector of the body at */
288 /* periapse. */
289
290 /* The only other consideration is determining an epoch, TP, of */
291 /* this state and the delta time ET - TP. */
292
293
294 /* Standard SPICE error handling. */
295
296 if (return_()) {
297 return 0;
298 } else {
299 chkin_("CONICS", (ftnlen)6);
300 }
301
302 /* Unpack the element vector. */
303
304 rp = elts[0];
305 ecc = elts[1];
306 inc = elts[2];
307 lnode = elts[3];
308 argp = elts[4];
309 m0 = elts[5];
310 t0 = elts[6];
311 mu = elts[7];
312
313 /* Handle all of the exceptions first. */
314
315 if (ecc < 0.) {
316 setmsg_("The eccentricity supplied was negative. Only positive value"
317 "s are meaningful. The value was #", (ftnlen)93);
318 errdp_("#", &ecc, (ftnlen)1);
319 sigerr_("SPICE(BADECCENTRICITY)", (ftnlen)22);
320 chkout_("CONICS", (ftnlen)6);
321 return 0;
322 }
323 if (rp <= 0.) {
324 setmsg_("The value of periapse range supplied was non-positive. Onl"
325 "y positive values are allowed. The value supplied was #. ", (
326 ftnlen)117);
327 errdp_("#", &rp, (ftnlen)1);
328 sigerr_("SPICE(BADPERIAPSEVALUE)", (ftnlen)23);
329 chkout_("CONICS", (ftnlen)6);
330 return 0;
331 }
332 if (mu <= 0.) {
333 setmsg_("The value of GM supplied was non-positive. Only positive v"
334 "alues are allowed. The value supplied was #. ", (ftnlen)105);
335 errdp_("#", &mu, (ftnlen)1);
336 sigerr_("SPICE(BADGM)", (ftnlen)12);
337 chkout_("CONICS", (ftnlen)6);
338 return 0;
339 }
340
341 /* First construct the orthonormal basis vectors that span the orbit */
342 /* plane. */
343
344 cosi = cos(inc);
345 sini = sin(inc);
346 cosn = cos(lnode);
347 sinn = sin(lnode);
348 cosw = cos(argp);
349 sinw = sin(argp);
350 snci = sinn * cosi;
351 cnci = cosn * cosi;
352 basisp[0] = cosn * cosw - snci * sinw;
353 basisp[1] = sinn * cosw + cnci * sinw;
354 basisp[2] = sini * sinw;
355 basisq[0] = -cosn * sinw - snci * cosw;
356 basisq[1] = -sinn * sinw + cnci * cosw;
357 basisq[2] = sini * cosw;
358
359 /* Next construct the state at periapse. */
360
361 /* The position at periapse is just BASISP scaled by the distance */
362 /* at periapse. */
363
364 /* The velocity must be constructed so that we can get an orbit */
365 /* of this shape. Recall that the magnitude of the specific angular */
366 /* momentum vector is given by DSQRT ( MU*RP*(1+ECC) ) */
367 /* The velocity will be given by V * BASISQ. But we must have the */
368 /* magnitude of the cross product of position and velocity be */
369 /* equal to DSQRT ( MU*RP*(1+ECC) ). So we must have */
370
371 /* RP*V = DSQRT( MU*RP*(1+ECC) ) */
372
373 /* so that: */
374
375 v = sqrt(mu * (ecc + 1.) / rp);
376 vscl_(&rp, basisp, pstate);
377 vscl_(&v, basisq, &pstate[3]);
378
379 /* Finally compute DT the elapsed time since the epoch of periapse. */
380 /* Ellipses first, since they are the most common. */
381
382 if (ecc < 1.) {
383
384 /* Recall that: */
385
386 /* N ( mean motion ) is given by DSQRT( MU / A**3 ). */
387 /* But since, A = RP / ( 1 - ECC ) ... */
388
389 ainvrs = (1. - ecc) / rp;
390 n = sqrt(mu * ainvrs) * ainvrs;
391 period = twopi_() / n;
392
393 /* In general the mean anomaly is given by */
394
395 /* M = (T - TP) * N */
396
397 /* Where TP is the time of periapse passage. M0 is the mean */
398 /* anomaly at time T0 so that */
399 /* Thus */
400
401 /* M0 = ( T0 - TP ) * N */
402
403 /* So TP = T0-M0/N hence the time since periapse at time ET */
404 /* is given by ET - T0 + M0/N. Finally, since elliptic orbits are */
405 /* periodic, we can mod this value by the period of the orbit. */
406
407 d__1 = *et - t0 + m0 / n;
408 dt = d_mod(&d__1, &period);
409
410 /* Hyperbolas next. */
411
412 } else if (ecc > 1.) {
413
414 /* Again, recall that: */
415
416 /* N ( mean motion ) is given by DSQRT( MU / |A**3| ). */
417 /* But since, |A| = RP / ( ECC - 1 ) ... */
418
419 ainvrs = (ecc - 1.) / rp;
420 n = sqrt(mu * ainvrs) * ainvrs;
421 dt = *et - t0 + m0 / n;
422
423 /* Finally, parabolas. */
424
425 } else {
426 n = sqrt(mu / (rp * 2.)) / rp;
427 dt = *et - t0 + m0 / n;
428 }
429
430 /* Now let PROP2B do the work of propagating the state. */
431
432 prop2b_(&mu, pstate, &dt, state);
433 chkout_("CONICS", (ftnlen)6);
434 return 0;
435 } /* conics_ */
436
437