1 /* spke15.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__6 = 6;
11
12 /* $Procedure SPKE15 ( Evaluate a type 15 SPK data record) */
spke15_(doublereal * et,doublereal * recin,doublereal * state)13 /* Subroutine */ int spke15_(doublereal *et, doublereal *recin, doublereal *
14 state)
15 {
16 /* System generated locals */
17 doublereal d__1;
18
19 /* Builtin functions */
20 double sqrt(doublereal), d_mod(doublereal *, doublereal *), d_sign(
21 doublereal *, doublereal *);
22
23 /* Local variables */
24 doublereal near__, dmdt;
25 extern /* Subroutine */ int vscl_(doublereal *, doublereal *, doublereal *
26 );
27 extern doublereal vdot_(doublereal *, doublereal *), vsep_(doublereal *,
28 doublereal *);
29 extern /* Subroutine */ int vequ_(doublereal *, doublereal *);
30 integer j2flg;
31 doublereal p, angle, dnode, z__;
32 extern /* Subroutine */ int chkin_(char *, ftnlen);
33 doublereal epoch, speed, dperi, theta, manom;
34 extern /* Subroutine */ int moved_(doublereal *, integer *, doublereal *),
35 errdp_(char *, doublereal *, ftnlen), vcrss_(doublereal *,
36 doublereal *, doublereal *);
37 extern doublereal twopi_(void);
38 extern logical vzero_(doublereal *);
39 extern /* Subroutine */ int vrotv_(doublereal *, doublereal *, doublereal
40 *, doublereal *);
41 doublereal oneme2, state0[6];
42 extern /* Subroutine */ int prop2b_(doublereal *, doublereal *,
43 doublereal *, doublereal *);
44 doublereal pa[3], gm, ta, dt;
45 extern doublereal pi_(void);
46 doublereal tp[3], pv[3], cosinc;
47 extern /* Subroutine */ int sigerr_(char *, ftnlen), vhatip_(doublereal *)
48 , chkout_(char *, ftnlen), vsclip_(doublereal *, doublereal *),
49 setmsg_(char *, ftnlen);
50 doublereal tmpsta[6], oj2;
51 extern logical return_(void);
52 doublereal ecc;
53 extern doublereal dpr_(void);
54 doublereal dot, rpl, k2pi;
55
56 /* $ Abstract */
57
58 /* Evaluates a single SPK data record from a segment of type 15 */
59 /* (Precessing Conic Propagation). */
60
61 /* $ Disclaimer */
62
63 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
64 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
65 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
66 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
67 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
68 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
69 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
70 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
71 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
72 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
73
74 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
75 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
76 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
77 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
78 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
79 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
80
81 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
82 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
83 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
84 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
85
86 /* $ Required_Reading */
87
88 /* SPK */
89
90 /* $ Keywords */
91
92 /* EPHEMERIS */
93
94 /* $ Declarations */
95 /* $ Brief_I/O */
96
97 /* Variable I/O Description */
98 /* -------- --- -------------------------------------------------- */
99 /* ET I Target epoch. */
100 /* RECIN I Data record. */
101 /* STATE O State (position and velocity). */
102
103 /* $ Detailed_Input */
104
105 /* ET is a target epoch, specified as ephemeris seconds past */
106 /* J2000, at which a state vector is to be computed. */
107
108 /* RECIN is a data record which, when evaluated at epoch ET, */
109 /* will give the state (position and velocity) of some */
110 /* body, relative to some center, in some inertial */
111 /* reference frame. */
112
113 /* The structure of RECIN is: */
114
115 /* RECIN(1) epoch of periapsis */
116 /* in ephemeris seconds past J2000. */
117 /* RECIN(2)-RECIN(4) unit trajectory pole vector */
118 /* RECIN(5)-RECIN(7) unit periapsis vector */
119 /* RECIN(8) semi-latus rectum---p in the */
120 /* equation: */
121
122 /* r = p/(1 + ECC*COS(Nu)) */
123
124 /* RECIN(9) eccentricity */
125 /* RECIN(10) J2 processing flag describing */
126 /* what J2 corrections are to be */
127 /* applied when the orbit is */
128 /* propagated. */
129
130 /* All J2 corrections are applied */
131 /* if this flag has a value that */
132 /* is not 1,2 or 3. */
133
134 /* If the value of the flag is 3 */
135 /* no corrections are done. */
136
137 /* If the value of the flag is 1 */
138 /* no corrections are computed for */
139 /* the precession of the line */
140 /* of apsides. However, regression */
141 /* of the line of nodes is */
142 /* performed. */
143
144 /* If the value of the flag is 2 */
145 /* no corrections are done for */
146 /* the regression of the line of */
147 /* nodes. However, precession of the */
148 /* line of apsides is performed. */
149
150 /* Note that J2 effects are computed */
151 /* only if the orbit is elliptic and */
152 /* does not intersect the central */
153 /* body. */
154
155 /* RECIN(11)-RECIN(13) unit central body pole vector */
156 /* RECIN(14) central body GM */
157 /* RECIN(15) central body J2 */
158 /* RECIN(16) central body radius */
159
160 /* Units are radians, km, seconds */
161
162 /* $ Detailed_Output */
163
164 /* STATE is the state produced by evaluating RECIN at ET. */
165 /* Units are km and km/sec. */
166
167 /* $ Parameters */
168
169 /* None. */
170
171 /* $ Files */
172
173 /* None. */
174
175 /* $ Exceptions */
176
177 /* 1) If the eccentricity is less than zero, the error */
178 /* 'SPICE(BADECCENTRICITY)' will be signalled. */
179
180 /* 2) If the semi-latus rectum is non-positive, the error */
181 /* 'SPICE(BADLATUSRECTUM)' is signalled. */
182
183 /* 3) If the pole vector, trajectory pole vector or periapsis vector */
184 /* has zero length, the error 'SPICE(BADVECTOR)' is signalled. */
185
186 /* 4) If the trajectory pole vector and the periapsis vector are */
187 /* not orthogonal, the error 'SPICE(BADINITSTATE)' is */
188 /* signalled. The test for orthogonality is very crude. The */
189 /* routine simply checks that the absolute value of the dot */
190 /* product of the unit vectors parallel to the trajectory pole */
191 /* and periapse vectors is less than 0.00001. This check is */
192 /* intended to catch blunders, not to enforce orthogonality to */
193 /* double precision tolerance. */
194
195 /* 5) If the mass of the central body is non-positive, the error */
196 /* 'SPICE(NONPOSITIVEMASS)' is signalled. */
197
198 /* 6) If the radius of the central body is negative, the error */
199 /* 'SPICE(BADRADIUS)' is signalled. */
200
201 /* $ Particulars */
202
203 /* This algorithm applies J2 corrections for precessing the */
204 /* node and argument of periapse for an object orbiting an */
205 /* oblate spheroid. */
206
207 /* Note the effects of J2 are incorporated only for elliptic */
208 /* orbits that do not intersect the central body. */
209
210 /* While the derivation of the effect of the various harmonics */
211 /* of gravitational field are beyond the scope of this header */
212 /* the effect of the J2 term of the gravity model are as follows */
213
214
215 /* The line of node precesses. Over one orbit average rate of */
216 /* precession, DNode/dNu, is given by */
217
218 /* 3 J2 */
219 /* dNode/dNu = - ----------------- DCOS( inc ) */
220 /* 2 (P/RPL)**2 */
221
222 /* (Since this is always less than zero for oblate spheroids, this */
223 /* should be called regression of nodes.) */
224
225 /* The line of apsides precesses. The average rate of precession */
226 /* DPeri/dNu is given by */
227 /* 3 J2 */
228 /* dPeri/dNu = ----------------- ( 5*DCOS ( inc ) - 1 ) */
229 /* 2 (P/RPL)**2 */
230
231 /* Details of these formulae are given in the Battin's book (see */
232 /* literature references below). */
233
234
235 /* It is assumed that this routine is used in conjunction with */
236 /* the routine SPKR15 as shown here: */
237
238 /* CALL SPKR15 ( HANDLE, DESCR, ET, RECIN ) */
239 /* CALL SPKE15 ( ET, RECIN, STATE ) */
240
241 /* where it is known in advance that the HANDLE, DESCR pair points */
242 /* to a type 15 data segment. */
243
244 /* $ Examples */
245
246 /* The SPKEnn routines are almost always used in conjunction with */
247 /* the corresponding SPKRnn routines, which read the records from */
248 /* SPK files. */
249
250 /* The data returned by the SPKRnn routine is in its rawest form, */
251 /* taken directly from the segment. As such, it will be meaningless */
252 /* to a user unless he/she understands the structure of the data type */
253 /* completely. Given that understanding, however, the SPKRnn */
254 /* routines might be used to examine raw segment data before */
255 /* evaluating it with the SPKEnn routines. */
256
257
258 /* C */
259 /* C Get a segment applicable to a specified body and epoch. */
260 /* C */
261 /* CALL SPKSFS ( BODY, ET, HANDLE, DESCR, IDENT, FOUND ) */
262
263 /* C */
264 /* C Look at parts of the descriptor. */
265 /* C */
266 /* CALL DAFUS ( DESCR, 2, 6, DCD, ICD ) */
267 /* CENTER = ICD( 2 ) */
268 /* REF = ICD( 3 ) */
269 /* TYPE = ICD( 4 ) */
270
271 /* IF ( TYPE .EQ. 15 ) THEN */
272
273 /* CALL SPKR15 ( HANDLE, DESCR, ET, RECORD ) */
274 /* . */
275 /* . Look at the RECORD data. */
276 /* . */
277 /* CALL SPKE15 ( ET, RECORD, STATE ) */
278 /* . */
279 /* . Check out the evaluated state. */
280 /* . */
281 /* END IF */
282
283 /* $ Restrictions */
284
285 /* None. */
286
287 /* $ Author_and_Institution */
288
289 /* K.R. Gehringer (JPL) */
290 /* S. Schlaifer (JPL) */
291 /* W.L. Taber (JPL) */
292
293 /* $ Literature_References */
294
295 /* [1] `Fundamentals of Celestial Mechanics', Second Edition 1989 */
296 /* by J.M.A. Danby; Willman-Bell, Inc., P.O. Box 35025 */
297 /* Richmond Virginia; pp 345-347. */
298
299 /* [2] `Astronautical Guidance', by Richard H. Battin. 1964 */
300 /* McGraw-Hill Book Company, San Francisco. pp 199 */
301
302 /* $ Version */
303
304 /* - SPICELIB Version 1.2.0, 02-SEP-2005 (NJB) */
305
306 /* Updated to remove non-standard use of duplicate arguments */
307 /* in VHAT, VROTV, and VSCL calls. */
308
309 /* - SPICELIB Version 1.1.0, 29-FEB-1996 (KRG) */
310
311 /* The declaration for the SPICELIB function PI is now */
312 /* preceded by an EXTERNAL statement declaring PI to be an */
313 /* external function. This removes a conflict with any */
314 /* compilers that have a PI intrinsic function. */
315
316 /* - SPICELIB Version 1.0.0, 15-NOV-1994 (WLT) (SS) */
317
318 /* -& */
319 /* $ Index_Entries */
320
321 /* evaluate type_15 spk segment */
322
323 /* -& */
324 /* $ Revisions */
325
326 /* - SPICELIB Version 1.2.0, 02-SEP-2005 (NJB) */
327
328 /* Updated to remove non-standard use of duplicate arguments */
329 /* in VHAT, VROTV, and VSCL calls. */
330
331 /* - SPICELIB Version 1.1.0, 29-FEB-1996 (KRG) */
332
333 /* The declaration for the SPICELIB function PI is now */
334 /* preceded by an EXTERNAL statement declaring PI to be an */
335 /* external function. This removes a conflict with any */
336 /* compilers that have a PI intrinsic function. */
337
338 /* - SPICELIB Version 1.0.0, 15-NOV-1994 (WLT) (SS) */
339
340 /* -& */
341
342 /* SPICELIB Functions */
343
344
345 /* Local Variables */
346
347
348 /* Standard SPICE error handling. */
349
350 if (return_()) {
351 return 0;
352 }
353 chkin_("SPKE15", (ftnlen)6);
354
355 /* Fetch the various entities from the input record, first the epoch. */
356
357 epoch = recin[0];
358
359 /* The trajectory pole vector. */
360
361 vequ_(&recin[1], tp);
362
363 /* The periapsis vector. */
364
365 vequ_(&recin[4], pa);
366
367 /* Semi-latus rectum ( P in the P/(1 + ECC*COS(Nu) ), */
368 /* and eccentricity. */
369
370 p = recin[7];
371 ecc = recin[8];
372
373 /* J2 processing flag. */
374
375 j2flg = (integer) recin[9];
376
377 /* Central body pole vector. */
378
379 vequ_(&recin[10], pv);
380
381 /* The central mass, J2 and radius of the central body. */
382
383 gm = recin[13];
384 oj2 = recin[14];
385 rpl = recin[15];
386
387 /* Check all the inputs here for obvious failures. Yes, perhaps */
388 /* this is overkill. However, there is a lot more computation */
389 /* going on in this routine so that the small amount of overhead */
390 /* here should not be significant. */
391
392 if (p <= 0.) {
393 setmsg_("The semi-latus rectum supplied to the SPK type 15 evaluator"
394 " was non-positive. This value must be positive. The value s"
395 "upplied was #.", (ftnlen)133);
396 errdp_("#", &p, (ftnlen)1);
397 sigerr_("SPICE(BADLATUSRECTUM)", (ftnlen)21);
398 chkout_("SPKE15", (ftnlen)6);
399 return 0;
400 } else if (ecc < 0.) {
401 setmsg_("The eccentricity supplied for a type 15 segment is negative"
402 ". It must be non-negative. The value supplied to the type 1"
403 "5 evaluator was #. ", (ftnlen)138);
404 errdp_("#", &ecc, (ftnlen)1);
405 sigerr_("SPICE(BADECCENTRICITY)", (ftnlen)22);
406 chkout_("SPKE15", (ftnlen)6);
407 return 0;
408 } else if (gm <= 0.) {
409 setmsg_("The mass supplied for the central body of a type 15 segment"
410 " was non-positive. Masses must be positive. The value suppl"
411 "ied was #. ", (ftnlen)130);
412 errdp_("#", &gm, (ftnlen)1);
413 sigerr_("SPICE(NONPOSITIVEMASS)", (ftnlen)22);
414 chkout_("SPKE15", (ftnlen)6);
415 return 0;
416 } else if (vzero_(tp)) {
417 setmsg_("The trajectory pole vector supplied to SPKE15 had length ze"
418 "ro. The most likely cause of this problem is a corrupted SPK"
419 " (ephemeris) file. ", (ftnlen)138);
420 sigerr_("SPICE(BADVECTOR)", (ftnlen)16);
421 chkout_("SPKE15", (ftnlen)6);
422 return 0;
423 } else if (vzero_(pa)) {
424 setmsg_("The periapse vector supplied to SPKE15 had length zero. The"
425 " most likely cause of this problem is a corrupted SPK (ephem"
426 "eris) file. ", (ftnlen)131);
427 sigerr_("SPICE(BADVECTOR)", (ftnlen)16);
428 chkout_("SPKE15", (ftnlen)6);
429 return 0;
430 } else if (vzero_(pv)) {
431 setmsg_("The central pole vector supplied to SPKE15 had length zero."
432 " The most likely cause of this problem is a corrupted SPK (e"
433 "phemeris) file. ", (ftnlen)135);
434 sigerr_("SPICE(BADVECTOR)", (ftnlen)16);
435 chkout_("SPKE15", (ftnlen)6);
436 return 0;
437 } else if (rpl < 0.) {
438 setmsg_("The central body radius was negative. It must be zero or po"
439 "sitive. The value supplied was #. ", (ftnlen)94);
440 errdp_("#", &rpl, (ftnlen)1);
441 sigerr_("SPICE(BADRADIUS)", (ftnlen)16);
442 chkout_("SPKE15", (ftnlen)6);
443 return 0;
444 }
445
446 /* Convert TP, PV and PA to unit vectors. */
447 /* (It won't hurt to polish them up a bit here if they are already */
448 /* unit vectors.) */
449
450 vhatip_(pa);
451 vhatip_(tp);
452 vhatip_(pv);
453
454 /* One final check. Make sure the pole and periapse vectors are */
455 /* orthogonal. (We will use a very crude check but this should */
456 /* rule out any obvious errors.) */
457
458 dot = vdot_(pa, tp);
459 if (abs(dot) > 1e-5) {
460 angle = vsep_(pa, tp) * dpr_();
461 setmsg_("The periapsis and trajectory pole vectors are not orthogona"
462 "l. The anglebetween them is # degrees. ", (ftnlen)98);
463 errdp_("#", &angle, (ftnlen)1);
464 sigerr_("SPICE(BADINITSTATE)", (ftnlen)19);
465 chkout_("SPKE15", (ftnlen)6);
466 return 0;
467 }
468
469 /* Compute the distance and speed at periapse. */
470
471 near__ = p / (ecc + 1.);
472 speed = sqrt(gm / p) * (ecc + 1.);
473
474 /* Next get the position at periapse ... */
475
476 vscl_(&near__, pa, state0);
477
478 /* ... and the velocity at periapsis. */
479
480 vcrss_(tp, pa, &state0[3]);
481 vsclip_(&speed, &state0[3]);
482
483 /* Determine the elapsed time from periapse to the requested */
484 /* epoch and propagate the state at periapsis to the epoch of */
485 /* interest. */
486
487 /* Note that we are making use of the following fact. */
488
489 /* If R is a rotation, then the states obtained by */
490 /* the following blocks of code are mathematically the */
491 /* same. (In reality they may differ slightly due to */
492 /* roundoff.) */
493
494 /* Code block 1. */
495
496 /* CALL MXV ( R, STATE0, STATE0 ) */
497 /* CALL MXV ( R, STATE0(4), STATE0(4) ) */
498 /* CALL PROP2B( GM, STATE0, DT, STATE ) */
499
500 /* Code block 2. */
501
502 /* CALL PROP2B( GM, STATE0, DT, STATE ) */
503 /* CALL MXV ( R, STATE, STATE ) */
504 /* CALL MXV ( R, STATE(4), STATE(4) ) */
505
506
507 /* This allows us to first compute the propagation of our initial */
508 /* state and then if needed perform the precession of the line */
509 /* of nodes and apsides by simply precessing the resulting state. */
510
511 dt = *et - epoch;
512 prop2b_(&gm, state0, &dt, state);
513
514 /* If called for, handle precession needed due to the J2 term. Note */
515 /* that the motion of the lines of nodes and apsides is formulated */
516 /* in terms of the true anomaly. This means we need the accumulated */
517 /* true anomaly in order to properly transform the state. */
518
519 if (j2flg != 3 && oj2 != 0. && ecc < 1. && near__ > rpl) {
520
521 /* First compute the change in mean anomaly since periapsis. */
522
523 /* Computing 2nd power */
524 d__1 = ecc;
525 oneme2 = 1. - d__1 * d__1;
526 dmdt = oneme2 / p * sqrt(gm * oneme2 / p);
527 manom = dmdt * dt;
528
529 /* Next compute the angle THETA such that THETA is between */
530 /* -pi and pi and such than MANOM = THETA + K*2*pi for */
531 /* some integer K. */
532
533 d__1 = twopi_();
534 theta = d_mod(&manom, &d__1);
535 if (abs(theta) > pi_()) {
536 d__1 = twopi_();
537 theta -= d_sign(&d__1, &theta);
538 }
539 k2pi = manom - theta;
540
541 /* We can get the accumulated true anomaly from the propagated */
542 /* state theta and the accumulated mean anomaly prior to this */
543 /* orbit. */
544
545 ta = vsep_(pa, state);
546 ta = d_sign(&ta, &theta);
547 ta += k2pi;
548
549 /* Determine how far the line of nodes and periapsis have moved. */
550
551 cosinc = vdot_(pv, tp);
552 /* Computing 2nd power */
553 d__1 = rpl / p;
554 z__ = ta * 1.5 * oj2 * (d__1 * d__1);
555 dnode = -z__ * cosinc;
556 /* Computing 2nd power */
557 d__1 = cosinc;
558 dperi = z__ * (d__1 * d__1 * 2.5 - .5);
559
560 /* Precess the periapsis by rotating the state vector about the */
561 /* trajectory pole */
562
563 if (j2flg != 1) {
564 vrotv_(state, tp, &dperi, tmpsta);
565 vrotv_(&state[3], tp, &dperi, &tmpsta[3]);
566 moved_(tmpsta, &c__6, state);
567 }
568
569 /* Regress the line of nodes by rotating the state */
570 /* about the pole of the central body. */
571
572 if (j2flg != 2) {
573 vrotv_(state, pv, &dnode, tmpsta);
574 vrotv_(&state[3], pv, &dnode, &tmpsta[3]);
575 moved_(tmpsta, &c__6, state);
576 }
577
578 /* We could perform the rotations above in the other order, */
579 /* but we would also have to rotate the pole before precessing */
580 /* the line of apsides. */
581
582 }
583
584 /* That's all folks. Check out and return. */
585
586 chkout_("SPKE15", (ftnlen)6);
587 return 0;
588 } /* spke15_ */
589
590