1 /* oscltx.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 /* Table of constant values */
9
10 static integer c__12 = 12;
11 static doublereal c_b5 = 0.;
12 static doublereal c_b6 = 1e-10;
13 static integer c__1 = 1;
14 static integer c__3 = 3;
15 static doublereal c_b13 = .66666666666666663;
16 static doublereal c_b14 = .33333333333333331;
17 static doublereal c_b15 = 1.5;
18
19 /* $Procedure OSCLTX ( Extended osculating elements from state ) */
oscltx_(doublereal * state,doublereal * et,doublereal * mu,doublereal * elts)20 /* Subroutine */ int oscltx_(doublereal *state, doublereal *et, doublereal *
21 mu, doublereal *elts)
22 {
23 /* Initialized data */
24
25 static logical pass1 = TRUE_;
26 static doublereal limit = 0.;
27
28 /* System generated locals */
29 doublereal d__1;
30
31 /* Builtin functions */
32 double pow_dd(doublereal *, doublereal *);
33
34 /* Local variables */
35 doublereal rmag, vmag, rhat[3], pole[3], xmat[9] /* was [3][3] */;
36 extern doublereal vdot_(doublereal *, doublereal *);
37 extern /* Subroutine */ int vequ_(doublereal *, doublereal *);
38 doublereal a, b, e, range;
39 extern /* Subroutine */ int chkin_(char *, ftnlen);
40 extern doublereal exact_(doublereal *, doublereal *, doublereal *),
41 dpmax_(void);
42 doublereal rperi[3];
43 extern /* Subroutine */ int vlcom_(doublereal *, doublereal *, doublereal
44 *, doublereal *, doublereal *);
45 doublereal c1, c2;
46 extern /* Subroutine */ int ucrss_(doublereal *, doublereal *, doublereal
47 *), unorm_(doublereal *, doublereal *, doublereal *);
48 extern doublereal vnorm_(doublereal *), twopi_(void);
49 extern logical failed_(void);
50 doublereal eccvec[3];
51 extern /* Subroutine */ int cleard_(integer *, doublereal *), recrad_(
52 doublereal *, doublereal *, doublereal *, doublereal *);
53 doublereal nu;
54 logical cmptau;
55 extern /* Subroutine */ int oscelt_(doublereal *, doublereal *,
56 doublereal *, doublereal *);
57 doublereal mucubr;
58 extern /* Subroutine */ int chkout_(char *, ftnlen), twovec_(doublereal *,
59 integer *, doublereal *, integer *, doublereal *);
60 extern logical return_(void);
61 doublereal ecc, dec, vel[3], tau;
62 extern /* Subroutine */ int mxv_(doublereal *, doublereal *, doublereal *)
63 ;
64
65 /* $ Abstract */
66
67 /* Determine the set of osculating conic orbital elements that */
68 /* corresponds to the state (position, velocity) of a body at some */
69 /* epoch. In additional to the classical elements, return the true */
70 /* anomaly, semi-major axis, and period, if applicable. */
71
72 /* $ Disclaimer */
73
74 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
75 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
76 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
77 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
78 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
79 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
80 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
81 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
82 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
83 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
84
85 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
86 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
87 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
88 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
89 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
90 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
91
92 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
93 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
94 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
95 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
96
97 /* $ Required_Reading */
98
99 /* None. */
100
101 /* $ Keywords */
102
103 /* CONIC */
104 /* ELEMENTS */
105 /* EPHEMERIS */
106
107 /* $ Declarations */
108
109 /* Include file oscltx.inc */
110
111 /* $ Abstract */
112
113 /* This file declares public parameters related to */
114 /* the SPICELIB routine OSCLTX. */
115
116 /* $ Version */
117
118 /* - SPICELIB Version 1.0.0, 25-JAN-2017 (NJB) */
119
120 /* -& */
121
122 /* Size of the OSCLTX output argument array ELTS. */
123
124 /* The output array ELTS is intended to contain unused space to hold */
125 /* additional elements that may be added in a later version of this */
126 /* routine. In order to maintain forward compatibility, user */
127 /* applications should declare ELTS as follows: */
128
129 /* DOUBLE PRECISION ELTS( OSCXSZ ) */
130
131
132 /* End of include file oscltx.inc */
133
134 /* $ Brief_I/O */
135
136 /* VARIABLE I/O DESCRIPTION */
137 /* -------- --- -------------------------------------------------- */
138 /* STATE I State of body at epoch of elements. */
139 /* ET I Epoch of elements. */
140 /* MU I Gravitational parameter (GM) of primary body. */
141 /* ELTS O Extended set of classical conic elements. */
142
143 /* $ Detailed_Input */
144
145 /* STATE is the state (position and velocity) of the body */
146 /* at some epoch. Components are x, y, z, dx/dt, dy/dt, */
147 /* dz/dt. STATE must be expressed relative to an */
148 /* inertial reference frame. Units are km and km/sec. */
149
150
151 /* ET is the epoch of the input state, in ephemeris seconds */
152 /* past J2000. */
153
154 /* 3 2 */
155 /* MU is the gravitational parameter (GM, km /sec ) of */
156 /* the primary body. */
157
158 /* $ Detailed_Output */
159
160 /* ELTS are equivalent conic elements describing the orbit */
161 /* of the body around its primary. The elements are, */
162 /* in order: */
163
164 /* RP Perifocal distance. */
165 /* ECC Eccentricity. */
166 /* INC Inclination. */
167 /* LNODE Longitude of the ascending node. */
168 /* ARGP Argument of periapsis. */
169 /* M0 Mean anomaly at epoch. */
170 /* T0 Epoch. */
171 /* MU Gravitational parameter. */
172 /* NU True anomaly at epoch. */
173 /* A Semi-major axis. A is set to zero if */
174 /* it is not computable. */
175 /* TAU Orbital period. Applicable only for */
176 /* elliptical orbits. Set to zero otherwise. */
177
178 /* The epoch of the elements is the epoch of the input */
179 /* state. Units are km, rad, rad/sec. The same elements */
180 /* are used to describe all three types (elliptic, */
181 /* hyperbolic, and parabolic) of conic orbit. */
182
183 /* User applications should declare ELTS using the */
184 /* parameter */
185
186 /* OSCXSZ */
187
188 /* See the Parameters section below. */
189
190
191 /* $ Parameters */
192
193 /* OSCXSZ is the size of the output elements array ELTS. OSCXSZ */
194 /* is declared in the Fortran include file */
195
196 /* oscltx.inc */
197
198 /* The output array ELTS is intended to contain unused */
199 /* space to hold additional elements that may be added */
200 /* in a later version of this routine. In order to */
201 /* maintain forward compatibility, user applications */
202 /* should declare ELTS as follows: */
203
204 /* DOUBLE PRECISION ELTS( OSCXSZ ) */
205
206
207 /* $ Exceptions */
208
209 /* 1) If MU is not positive, the error SPICE(NONPOSITIVEMASS) */
210 /* is signaled. */
211
212 /* 2) If the specific angular momentum vector derived from STATE */
213 /* is the zero vector, the error SPICE(DEGENERATECASE) */
214 /* is signaled. */
215
216 /* 3) If the position or velocity vectors derived from STATE */
217 /* is the zero vector, the error SPICE(DEGENERATECASE) */
218 /* is signaled. */
219
220 /* 4) If the inclination is determined to be zero or 180 degrees, */
221 /* the longitude of the ascending node is set to zero. */
222
223 /* 5) If the eccentricity is determined to be zero, the argument of */
224 /* periapse is set to zero. */
225
226 /* 6) If the eccentricity of the orbit is very close to but not */
227 /* equal to zero, the argument of periapse may not be accurately */
228 /* determined. */
229
230 /* 7) For inclinations near but not equal to 0 or 180 degrees, */
231 /* the longitude of the ascending node may not be determined */
232 /* accurately. The argument of periapse and mean anomaly may */
233 /* also be inaccurate. */
234
235 /* 8) For eccentricities very close to but not equal to 1, the */
236 /* results of this routine are unreliable. */
237
238 /* 9) If the specific angular momentum vector is non-zero but */
239 /* "close" to zero, the results of this routine are unreliable. */
240
241 /* 10) If STATE is expressed relative to a non-inertial reference */
242 /* frame, the resulting elements are invalid. No error checking */
243 /* is done to detect this problem. */
244
245 /* 11) The semi-major axis and period may not be computable for */
246 /* orbits having eccentricity too close to 1. If the semi-major */
247 /* axis is not computable, both it and the period are set to zero. */
248 /* If the period is not computable, it is set to zero. */
249
250 /* $ Files */
251
252 /* None. */
253
254 /* $ Particulars */
255
256 /* This routine returns in the first 8 elements of the array ELTS */
257 /* the outputs computed by OSCELT, and in addition returns in */
258 /* elements 9-11 the quantities: */
259
260 /* ELTS(9) true anomaly at ET, in radians. */
261
262 /* ELTS(10) orbital semi-major axis at ET, in km. Valid */
263 /* if and only if this value is non-zero. */
264
265 /* The semi-major axis won't be computable if the */
266 /* eccentricity of the orbit is too close to 1. */
267 /* In this case A is set to zero. */
268
269 /* ELTS(11) orbital period. If the period is not computable, */
270 /* TAU is set to zero. */
271
272 /* The SPICELIB routine CONICS is an approximate inverse of this */
273 /* routine: CONICS maps a set of osculating elements and a time to a */
274 /* state vector. */
275
276 /* $ Examples */
277
278 /* Let VINIT contain the initial state of a spacecraft relative to */
279 /* the center of a planet at epoch ET, and let GM be the gravitation */
280 /* parameter of the planet. The call */
281
282 /* CALL OSCLTX ( VINIT, ET, GM, ELTS ) */
283
284 /* produces a set of osculating elements describing the nominal */
285 /* orbit that the spacecraft would follow in the absence of all */
286 /* other bodies in the solar system. */
287
288 /* Now let STATE contain the state of the same spacecraft at some */
289 /* other epoch, LATER. The difference between this state and the */
290 /* state predicted by the nominal orbit at the same epoch can be */
291 /* computed as follows. */
292
293 /* CALL CONICS ( ELTS, LATER, NOMINAL ) */
294 /* CALL VSUBG ( NOMINAL, STATE, 6, DIFF ) */
295
296 /* WRITE (*,*) 'Perturbation in x, dx/dt = ', DIFF(1), DIFF(4) */
297 /* WRITE (*,*) ' y, dy/dt = ', DIFF(2), DIFF(5) */
298 /* WRITE (*,*) ' z, dz/dt = ', DIFF(3), DIFF(6) */
299
300 /* $ Restrictions */
301
302 /* 1) The input state vector must be expressed relative to an */
303 /* inertial reference frame. */
304
305 /* 2) Osculating elements are generally not useful for */
306 /* high-accuracy work. */
307
308 /* 3) Accurate osculating elements may be difficult to derive for */
309 /* near-circular or near-equatorial orbits. Osculating elements */
310 /* for such orbits should be used with caution. */
311
312 /* 4) Extracting osculating elements from a state vector is a */
313 /* mathematically simple but numerically challenging task. The */
314 /* mapping from a state vector to equivalent elements is */
315 /* undefined for certain state vectors, and the mapping is */
316 /* difficult to implement with finite precision arithmetic for */
317 /* states near the subsets of R6 where singularities occur. */
318
319 /* In general, the elements found by this routine can have */
320 /* two kinds of problems: */
321
322 /* - The elements are not accurate but still represent */
323 /* the input state accurately. The can happen in */
324 /* cases where the inclination is near zero or 180 */
325 /* degrees, or for near-circular orbits. */
326
327 /* - The elements are garbage. This can occur when */
328 /* the eccentricity of the orbit is close to but */
329 /* not equal to 1. In general, any inputs that cause */
330 /* great loss of precision in the computation of the */
331 /* specific angular momentum vector or the eccentricity */
332 /* vector will result in invalid outputs. */
333
334 /* For further details, see the Exceptions section. */
335
336 /* Users of this routine should carefully consider whether */
337 /* it is suitable for their applications. One recommended */
338 /* "sanity check" on the outputs is to supply them to the */
339 /* SPICELIB routine CONICS and compare the resulting state */
340 /* vector with the one supplied to this routine. */
341
342 /* $ Literature_References */
343
344 /* [1] Roger Bate, Fundamentals of Astrodynamics, Dover, 1971. */
345
346 /* $ Author_and_Institution */
347
348 /* N.J. Bachman (JPL) */
349 /* K.R. Gehringer (JPL) */
350 /* I.M. Underwood (JPL) */
351 /* E.D. Wright (JPL) */
352
353 /* $ Version */
354
355 /* - SPICELIB Version 1.0.0, 02-FEB-2017 (NJB) */
356
357 /* 12-MAR-2015 (NJB) */
358
359 /* Re-arranged test for small E to avoid overflow. */
360 /* Changed definition of B to make the maximum value */
361 /* of TAU equal to LIMIT. Removed test comparing */
362 /* E/LIMIT to RMAG. */
363
364 /* 11-NOV-2014 (NJB) */
365
366 /* Original version. Based on OSCELT version 1.3.1, */
367 /* 28-FEB-2008 */
368
369 /* -& */
370 /* $ Index_Entries */
371
372 /* extended conic elements from state */
373 /* extended osculating elements from state */
374 /* convert state to extended osculating elements */
375
376 /* -& */
377
378 /* SPICELIB functions */
379
380
381 /* Local parameters */
382
383
384 /* Local variables */
385
386 if (return_()) {
387 return 0;
388 }
389 chkin_("OSCLTX", (ftnlen)6);
390 if (pass1) {
391 limit = dpmax_() / 200.;
392 pass1 = FALSE_;
393 }
394
395 /* Initialize the members of ELTS that won't be set by OSCELT. This */
396 /* is necessary because this routine may return early if A or TAU */
397 /* cannot be computed. */
398
399 cleard_(&c__12, &elts[8]);
400 nu = 0.;
401 tau = 0.;
402 a = 0.;
403
404 /* Compute osculating elements using OSCELT. Take advantage */
405 /* of OSCELT's error handling. */
406
407 oscelt_(state, et, mu, elts);
408 if (failed_()) {
409 chkout_("OSCLTX", (ftnlen)6);
410 return 0;
411 }
412
413 /* Due to checks made by OSCELT, we know that MU and RMAG */
414 /* are strictly positive. */
415
416
417 /* Compute the true anomaly at ET. For the non-degenerate case, we */
418 /* use equation 2.4-5 from [1] to compute the eccentricity vector: */
419
420 /* _ 1 2 mu _ _ _ _ */
421 /* e = ---- * [ ( v - ---- ) r - <r, v> v ] (1) */
422 /* mu r */
423
424
425 /* where */
426
427 /* _ */
428 /* e is the eccentricity vector */
429
430 /* _ _ */
431 /* r, v are, respectively, the object's position and */
432 /* velocity vectors expressed in the inertial */
433 /* reference frame of the input state. */
434 /* _ _ */
435 /* r, v are, respectively, the magnitudes of r and v */
436
437 /* mu is the GM value of the center of motion */
438
439
440
441 unorm_(state, rhat, &rmag);
442 vequ_(&state[3], vel);
443 vmag = vnorm_(vel);
444 ecc = exact_(&elts[1], &c_b5, &c_b6);
445 if (ecc != 0.) {
446 /* _ _ */
447 /* C1 is the coefficient of r/||r|| in (1) above: */
448
449 /* Computing 2nd power */
450 d__1 = vmag;
451 c1 = 1. / *mu * (rmag * (d__1 * d__1) - *mu);
452 /* _ */
453 /* C2 is the coefficient of v in (1) above: */
454
455 c2 = -(1. / *mu) * vdot_(state, vel);
456
457 /* ECCVEC is the eccentricity vector: */
458
459 vlcom_(&c1, rhat, &c2, vel, eccvec);
460
461 /* POLE is the orbital pole unit vector. */
462
463 ucrss_(state, vel, pole);
464
465 /* Compute the transformation from the frame of the input state */
466 /* to the perifocal frame. */
467
468 twovec_(eccvec, &c__1, pole, &c__3, xmat);
469 if (failed_()) {
470 chkout_("OSCLTX", (ftnlen)6);
471 return 0;
472 }
473
474 /* Transform the input position to the perifocal frame and */
475 /* compute its "right ascension." This is the true anomaly NU, */
476 /* confined to the range [0, 2*pi) radians. */
477
478 mxv_(xmat, state, rperi);
479 recrad_(rperi, &range, &nu, &dec);
480 } else {
481
482 /* The orbit is circular or nearly so. The mean anomaly can */
483 /* be used as the true anomaly. */
484
485 nu = elts[5];
486 }
487 elts[8] = nu;
488
489 /* The semi-major axis A satisfies the relationship */
490
491 /* mu */
492 /* E = - ---- (2) */
493 /* 2A */
494
495 /* where E represents the specific energy of the orbit, as long */
496 /* as A is computable. */
497
498 /* If the orbit is not parabolic or too close to parabolic, */
499 /* we can recover A from (2) above: */
500
501 /* mu */
502 /* A = - ---- */
503 /* 2E */
504
505 /* To make sure that A is computable, we require */
506
507
508
509 /* |mu/E| < LIMIT */
510
511 /* so */
512
513 /* |mu/LIMIT| < |E| */
514
515 /* We compute E from the specific kinetic and potential */
516 /* energy of the orbit: */
517
518 /* 2 */
519 /* v mu */
520 /* E = --- - ---- (3) */
521 /* 2 r */
522
523 /* The second term on the right hand side is computable if */
524
525 /* mu/LIMIT < r */
526
527
528 if (elts[1] == 1.) {
529
530 /* This is the parabolic case. */
531
532 chkout_("OSCLTX", (ftnlen)6);
533 return 0;
534 }
535 if (rmag <= *mu / limit) {
536
537 /* We need RMAG > MU/LIMIT to make E computable. */
538
539 /* We can't safely compute E. */
540
541 chkout_("OSCLTX", (ftnlen)6);
542 return 0;
543 }
544
545 /* We assume V and MU are small enough not to cause overflow. */
546
547 /* Computing 2nd power */
548 d__1 = vmag;
549 e = d__1 * d__1 * .5 - *mu / rmag;
550 if (abs(e) < abs(*mu) / limit) {
551
552 /* We can't safely compute A. The case of E equal to 0 */
553 /* gets the code to this point. */
554
555 chkout_("OSCLTX", (ftnlen)6);
556 return 0;
557 }
558
559 /* At this point we can compute A (which is stored in ELTS(10)). */
560
561 a = -(*mu) / (e * 2);
562 elts[9] = a;
563
564 /* If the orbit is elliptical, compute the period. */
565
566 ecc = elts[1];
567 if (ecc < 1.) {
568
569 /* The period is given by the formula */
570
571 /* 3 1/2 */
572 /* tau = 2 * pi * ( a / mu ) */
573
574 /* We need to make sure this computation doesn't */
575 /* overflow. */
576
577 /* If mu >= 1 then we can express the constraint */
578 /* on a as */
579
580 /* 1/3 2/3 */
581 /* ( a / mu ) < ( LIMIT / (2*pi) ) */
582
583
584 /* Otherwise, we require that */
585
586 /* 2/3 1/3 */
587 /* a < ( LIMIT / (2*pi) ) * mu */
588
589
590
591 /* Note that the case of non-positive MU has already */
592 /* been ruled out. */
593
594 d__1 = limit / twopi_();
595 b = pow_dd(&d__1, &c_b13);
596 mucubr = pow_dd(mu, &c_b14);
597 if (*mu >= 1.) {
598 cmptau = a / mucubr < b;
599 } else {
600 cmptau = a < b * mucubr;
601 }
602 if (cmptau) {
603 d__1 = a / mucubr;
604 tau = twopi_() * pow_dd(&d__1, &c_b15);
605 }
606 }
607
608 /* Assign the remaining members of ELTS. */
609
610 elts[10] = tau;
611 chkout_("OSCLTX", (ftnlen)6);
612 return 0;
613 } /* oscltx_ */
614
615