1 /* dpspce.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 doublereal c_b19 = .66666666666666663;
11 static doublereal c_b20 = 3.5;
12 static doublereal c_b22 = 1.5;
13 static doublereal c_b23 = 1.;
14 static doublereal c_b25 = 0.;
15
16 /* $Procedure DPSPCE ( Propagate a two line element set for deep space ) */
dpspce_(doublereal * time,doublereal * geophs,doublereal * elems,doublereal * state)17 /* Subroutine */ int dpspce_(doublereal *time, doublereal *geophs, doublereal
18 *elems, doublereal *state)
19 {
20 /* Initialized data */
21
22 static logical doinit = TRUE_;
23 static logical first = TRUE_;
24
25 /* System generated locals */
26 integer i__1, i__2;
27 doublereal d__1, d__2;
28
29 /* Builtin functions */
30 integer s_rnge(char *, integer, char *, integer);
31 double pow_dd(doublereal *, doublereal *), cos(doublereal), sqrt(
32 doublereal), sin(doublereal), d_mod(doublereal *, doublereal *),
33 atan2(doublereal, doublereal);
34
35 /* Local variables */
36 static doublereal coef, eeta, aodp, delo, capu, uang, xmdf, xinc, xmam,
37 aynl, elsq, temp;
38 static logical cont;
39 static doublereal rdot, cosu, sinu, coef1, t2cof, temp1, temp2, temp3,
40 temp4, temp5, cos2u, temp6;
41 extern /* Subroutine */ int zzdpinit_(doublereal *, doublereal *,
42 doublereal *, doublereal *, doublereal *, doublereal *);
43 static doublereal sin2u, a, e;
44 static integer i__;
45 static doublereal m[3], n[3], s, u[3], v[3], betal, scale, betao;
46 extern /* Subroutine */ int chkin_(char *, ftnlen);
47 static doublereal epoch, ecose, aycof, esine, a3ovk2, tempa, tempe, bstar,
48 cosio, xincl, etasq, rfdot, sinio, a1, rdotk, c1, c2, cosuk, c4,
49 qoms24, sinuk, templ, x1m5th, x1mth2, x3thm1, x7thm1, psisq,
50 xinck, xlcof, xmdot, xnode, xnodp;
51 extern doublereal twopi_(void);
52 static doublereal s4;
53 extern /* Subroutine */ int vlcom_(doublereal *, doublereal *, doublereal
54 *, doublereal *, doublereal *);
55 static doublereal betao2, theta2, ae, xhdot1, ao, em, eo, qoms2t, pl,
56 omgadf, rk, qo, uk, so;
57 extern doublereal halfpi_(void);
58 static doublereal xl, xn, omegao;
59 extern /* Subroutine */ int latrec_(doublereal *, doublereal *,
60 doublereal *, doublereal *);
61 static doublereal perige, xnodcf, xnoddf, tsince, xnodek, omgdot, rfdotk,
62 xnodeo;
63 extern /* Subroutine */ int chkout_(char *, ftnlen);
64 static doublereal ck2, lstelm[10], ck4, cosepw, sinepw, xkmper, xnodot,
65 lstphs[8];
66 extern logical return_(void);
67 static doublereal pinvsq, xj2, xj3, xj4, eta, axn, xke, ayn, epw, tsi,
68 xll, xmo, xno, tsq, xlt, del1;
69 extern /* Subroutine */ int zzdpsec_(doublereal *, doublereal *,
70 doublereal *, doublereal *, doublereal *, doublereal *,
71 doublereal *, doublereal *, doublereal *);
72 static doublereal pio2;
73 extern /* Subroutine */ int zzdpper_(doublereal *, doublereal *,
74 doublereal *, doublereal *, doublereal *, doublereal *);
75 static doublereal pix2;
76
77 /* $ Abstract */
78
79 /* This routine propagates NORAD two-line element data for */
80 /* earth orbiting deep space vehicles (a vehicle with an */
81 /* orbital period more than 225 minutes). */
82
83 /* $ Disclaimer */
84
85 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
86 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
87 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
88 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
89 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
90 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
91 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
92 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
93 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
94 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
95
96 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
97 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
98 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
99 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
100 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
101 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
102
103 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
104 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
105 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
106 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
107
108 /* $ Required_Reading */
109
110 /* None. */
111
112 /* $ Keywords */
113
114 /* EPHEMERIS */
115 /* TWO LINE ELEMENTS */
116 /* DEEP SPACE PROPAGATOR */
117
118 /* $ Declarations */
119 /* $ Brief_I/O */
120
121 /* VARIABLE I/O DESCRIPTION */
122 /* -------- --- -------------------------------------------------- */
123 /* TIME I Time for state evaluation in seconds past ephemeris */
124 /* epoch J2000. */
125 /* GEOPHS I The array of geophysical constants */
126 /* ELEMS I Array of orbit elements */
127 /* STATE O State vector at TIME */
128
129 /* $ Detailed_Input */
130
131 /* TIME is the epoch in seconds past ephemeris epoch J2000 */
132 /* to produced a state from the input elements. */
133
134 /* GEOPHS is a collection of 8 geophysical constants needed */
135 /* for computing a state. The order of these */
136 /* constants must be: */
137
138 /* GEOPHS(1) = J2 gravitational harmonic for earth */
139 /* GEOPHS(2) = J3 gravitational harmonic for earth */
140 /* GEOPHS(3) = J4 gravitational harmonic for earth */
141
142 /* These first three constants are dimensionless. */
143
144 /* GEOPHS(4) = KE: Square root of the GM for earth where */
145 /* GM is expressed in earth radii cubed per */
146 /* minutes squared. */
147
148 /* GEOPHS(5) = QO: Low altitude bound for atmospheric */
149 /* model in km. */
150
151 /* GEOPHS(6) = SO: High altitude bound for atmospheric */
152 /* model in km. */
153
154
155 /* GEOPHS(7) = RE: Equatorial radius of the earth in km. */
156
157
158 /* GEOPHS(8) = AE: Distance units/earth radius */
159 /* (normally 1) */
160
161 /* Below are currently recommended values for these */
162 /* items: */
163
164 /* J2 = 1.082616D-3 */
165 /* J3 = -2.53881D-6 */
166 /* J4 = -1.65597D-6 */
167
168 /* The next item is the square root of GM for the */
169 /* earth given in units of earth-radii**1.5/Minute */
170
171 /* KE = 7.43669161D-2 */
172
173 /* The next two items define the top and */
174 /* bottom of the atmospheric drag model */
175 /* used by the type 10 ephemeris type. */
176 /* Don't adjust these unless you understand */
177 /* the full implications of such changes. */
178
179 /* QO = 120.0D0 */
180 /* SO = 78.0D0 */
181
182 /* The ER value is the equatorial radius in km */
183 /* of the earth as used by NORAD. */
184
185 /* ER = 6378.135D0 */
186
187 /* The value of AE is the number of */
188 /* distance units per earth radii used by */
189 /* the NORAD state propagation software. */
190 /* The value is 1 unless you've got */
191 /* a very good understanding of the NORAD */
192 /* routine SGP4 and the affect of changing */
193 /* this value.. */
194
195 /* AE = 1.0D0 */
196
197 /* ELEMS is an array containing two-line element data */
198 /* as prescribed below. The elements XNDD6O and BSTAR */
199 /* must have been scaled by the proper exponent stored */
200 /* in the two line elements set. Moreover, the */
201 /* various items must be converted to the units shown */
202 /* here. */
203
204 /* ELEMS ( 1 ) = XNDT2O in radians/minute**2 */
205 /* ELEMS ( 2 ) = XNDD6O in radians/minute**3 */
206 /* ELEMS ( 3 ) = BSTAR */
207 /* ELEMS ( 4 ) = XINCL in radians */
208 /* ELEMS ( 5 ) = XNODEO in radians */
209 /* ELEMS ( 6 ) = EO */
210 /* ELEMS ( 7 ) = OMEGAO in radians */
211 /* ELEMS ( 8 ) = XMO in radians */
212 /* ELEMS ( 9 ) = XNO in radians/minute */
213 /* ELEMS ( 10 ) = EPOCH of the elements in seconds */
214 /* past ephemeris epoch J2000. */
215
216 /* $ Detailed_Output */
217
218 /* STATE A 6 vector containing the X, Y, Z, Vx, Vy, Vz */
219 /* coordinates in the inertial frame (double */
220 /* precision). */
221
222 /* $ Parameters */
223
224 /* None. */
225
226 /* $ Exceptions */
227
228 /* Error free. */
229
230 /* $ Files */
231
232 /* None. */
233
234 /* $ Particulars */
235
236 /* This subroutine is an extensive rewrite of the SDP4 */
237 /* routine as described in the Spacetrack 3 report. All common */
238 /* blocks were removed and all variables are explicitly defined. */
239
240 /* The removal of common blocks causes the set of routines to */
241 /* execute slower than the original version of SDP4. However the */
242 /* stability improves especially as concerns memory and */
243 /* expanded internal documentation. */
244
245 /* Trivial or redundant variables have been eliminated. */
246
247 /* R removed, occurrence replaced with RK */
248 /* E6A renamed TOL */
249 /* THETA4 removed, relevant equation recast in Horner's form */
250 /* i.e. something like x^4 + x^2 -> x^2 ( x^2 + 1 ) */
251 /* U renamed UANG, U is now a euclidean 3 vector. */
252 /* Ux,Uy,Uz removed, replaced with 3-vector U */
253 /* Vx,Vy,Vz removed, replaced with 3-vector V */
254 /* OMEGAQ removed, usage replaced with OMEGAO */
255 /* OMGDT removed, same variable as OMGDOT, so all occurrences */
256 /* replaced with OMGDOT */
257 /* SSL,SSG replaced with the 5-vector SSX */
258 /* SSH,SSE */
259 /* SSI */
260
261 /* Three functions present in the original Spacetrack report, ACTAN, */
262 /* FMOD2P and THETAG, have been either replaced with an intrinsic */
263 /* FORTRAN function (ACTAN -> DATAN2, FMOD2P -> DMOD) or recoded */
264 /* using SPICELIB calls (THETAG). */
265
266 /* The code at the end of this subroutine which calculates */
267 /* orientation vectors, was replaced with a set of calls to */
268 /* SPICELIB vector routines. */
269
270 /* A direct comparison of output from the original Spacetrack 3 code */
271 /* and these NAIF routines for the same elements and time parameters */
272 /* will produce unacceptably different results. */
273
274 /* $ Examples */
275
276
277 /* C--- Load the geophysical constants kernel and the leapsecond */
278 /* kernel */
279 /* CALL FURNSH( '/Users/ewright/lib/geophysical.ker' ) */
280 /* CALL FURNSH( '/kernels/gen/lsk/naif0008.tls' ) */
281
282
283 /* C--- Define a vehicle element array, TDRS 4 Geosynch */
284 /* LINES( 1 ) = '1 19883U 89021B 97133.05943164 -.00000277 ' */
285 /* .// '00000-0 10000-3 0 3315' */
286 /* LINES( 2 ) = '2 19883 0.5548 86.7278 0001786 312.2904 ' */
287 /* .// '172.2391 1.00269108202415' */
288
289
290 /* C--- Identify the earliest first year for the elements */
291 /* FRSTYR = 1988 */
292
293
294 /* C--- Parse the elements to something SPICE can use */
295 /* CALL GETELM ( FRSTYR, LINES, EPOCH, ELEMS ) */
296
297
298 /* C--- Final time past epoch, 1400 mins (in seconds) */
299 /* TF = 1440.D0 * 60.D0 */
300
301 /* C--- Step size for elements output 360 mins (in seconds) */
302 /* DELT = 360.D0 * 60.D0 */
303
304 /* C--- Start time keyed off epoch */
305 /* TIME = EPOCH - 2.D0 * DELT */
306
307 /* DO WHILE ( DABS(TIME - EPOCH) .LE. DABS(TF) ) */
308
309 /* CALL DPSPCE ( TIME, GEOPHS, ELEMS, STATE ) */
310
311 /* WRITE(*, FMT ='(7F17.8)' ) (TIME-EPOCH)/60.D0, */
312 /* . (STATE(I),I=1,6) */
313
314 /* TIME = TIME + DELT */
315
316 /* END DO */
317
318
319 /* $ Restrictions */
320
321 /* None. */
322
323 /* $ Literature_References */
324
325 /* Hoots, Felix R., Ronald L. Roehrich (31 December 1988). "Models */
326 /* for Propagation of NORAD Element Sets". United States Department */
327 /* of Defense Spacetrack Report (3). */
328
329 /* $ Author_and_Institution */
330
331 /* E.D. Wright (JPL) */
332
333 /* $ Version */
334
335 /* - SPICELIB Version 2.0.0, 23-JAN-2013 (EDW) */
336
337 /* Corrected initialization block error. The ZZDPINIT call */
338 /* causes a side-effect required for each DPSPCE call. */
339 /* The ZZDPINIT call now occurs outside the initialization */
340 /* block. Note from designer, side-effects are bad. */
341
342 /* Added proper citation for Hoots paper. */
343
344 /* - SPICELIB Version 1.2.2, 22-AUG-2006 (EDW) */
345
346 /* Replaced references to LDPOOL with references */
347 /* to FURNSH. */
348
349 /* - SPICELIB Version 1.2.1, DEC-27-2000 (EDW) */
350
351 /* Corrected error in header documentation. Horner's Rule */
352 /* not Butcher's. */
353
354 /* - SPICELIB Version 1.2.0, MAR-24-1999 (EDW) */
355
356 /* Documentation expanded to include modifications made */
357 /* to private routines. Some english errors corrected. */
358
359 /* Alphabetized variable declaration lists. */
360
361 /* Temporary variable TEMP removed. OMGDOT argument added to */
362 /* ZZDPSEC call. */
363
364 /* - SPICELIB Version 1.1.0, OCT-05-1998 (WLT) */
365
366 /* Forced initialization section until we can figure out */
367 /* why it doesn't work on SUNs. */
368
369 /* - SPICELIB Version 1.0.1, MAR-11-1998 (EDW) */
370
371 /* Corrected error in header describing GEOPHS array. */
372
373 /* - SPICELIB Version 1.0.0, NOV-11-1998 (EDW) */
374
375 /* -& */
376 /* $ Index_Entries */
377
378 /* NORAD two line elements deep space evaluator */
379
380 /* -& */
381
382 /* Local variables */
383
384
385 /* Define parameters for convergence tolerance and the value for 2/3, */
386 /* 0 and 1. */
387
388
389 /* The geophysical Quantities */
390
391
392 /* Elements */
393
394
395 /* Other quantities */
396
397
398 /* SPICELIB routines */
399
400
401 /* Save everything. */
402
403
404 /* Set initialization flags */
405
406
407 /* Standard SPICE error handling. */
408
409 if (return_()) {
410 return 0;
411 } else {
412 chkin_("DPSPCE", (ftnlen)6);
413 }
414
415 /* If this is the very first time into this routine, set these */
416 /* values. */
417
418 if (first) {
419 pix2 = twopi_();
420 pio2 = halfpi_();
421 first = FALSE_;
422 }
423
424 /* If initialization flag is FALSE, then this is not the first */
425 /* call to this routine. Check the stuff. */
426
427 if (! doinit) {
428
429 /* Check whether the current and last constants and elements */
430 /* match. If not, we need to reinitialize everything */
431 /* since the propagation is dependent on the value of these */
432 /* arrays. */
433
434 for (i__ = 1; i__ <= 8; ++i__) {
435 if (lstphs[(i__1 = i__ - 1) < 8 && 0 <= i__1 ? i__1 : s_rnge(
436 "lstphs", i__1, "dpspce_", (ftnlen)547)] != geophs[(i__2 =
437 i__ - 1) < 8 && 0 <= i__2 ? i__2 : s_rnge("geophs", i__2,
438 "dpspce_", (ftnlen)547)]) {
439 doinit = TRUE_;
440 }
441 }
442 for (i__ = 1; i__ <= 10; ++i__) {
443 if (lstelm[(i__1 = i__ - 1) < 10 && 0 <= i__1 ? i__1 : s_rnge(
444 "lstelm", i__1, "dpspce_", (ftnlen)556)] != elems[(i__2 =
445 i__ - 1) < 10 && 0 <= i__2 ? i__2 : s_rnge("elems", i__2,
446 "dpspce_", (ftnlen)556)]) {
447 doinit = TRUE_;
448 }
449 }
450 }
451
452 /* Initialization block. Always called on the initial entry and */
453 /* anytime the geophysical or elements array changes. */
454
455 if (doinit) {
456 doinit = FALSE_;
457
458 /* Retrieve the geophysical constants from the GEOPHS array */
459
460 xj2 = geophs[0];
461 xj3 = geophs[1];
462 xj4 = geophs[2];
463 xke = geophs[3];
464 qo = geophs[4];
465 so = geophs[5];
466 xkmper = geophs[6];
467 ae = geophs[7];
468
469 /* Save the geophysical constants for later comparison */
470
471 for (i__ = 1; i__ <= 8; ++i__) {
472 lstphs[(i__1 = i__ - 1) < 8 && 0 <= i__1 ? i__1 : s_rnge("lstphs",
473 i__1, "dpspce_", (ftnlen)590)] = geophs[(i__2 = i__ - 1)
474 < 8 && 0 <= i__2 ? i__2 : s_rnge("geophs", i__2, "dpspce_"
475 , (ftnlen)590)];
476 }
477
478 /* Unpack the elements array. */
479
480 bstar = elems[2];
481 xincl = elems[3];
482 xnodeo = elems[4];
483 eo = elems[5];
484 omegao = elems[6];
485 xmo = elems[7];
486 xno = elems[8];
487 epoch = elems[9];
488
489 /* Save the elements for later comparison */
490
491 for (i__ = 1; i__ <= 10; ++i__) {
492 lstelm[(i__1 = i__ - 1) < 10 && 0 <= i__1 ? i__1 : s_rnge("lstelm"
493 , i__1, "dpspce_", (ftnlen)610)] = elems[(i__2 = i__ - 1)
494 < 10 && 0 <= i__2 ? i__2 : s_rnge("elems", i__2, "dpspce_"
495 , (ftnlen)610)];
496 }
497
498 /* Set common variables, the init flag and calculate the */
499 /* WGS-72 physical and geopotential constants */
500
501 /* CK2 = 0.5 * J2 * AE^2 */
502 /* CK4 = -0.375 * J4 * AE^4 */
503
504 /* These are values calculated only once and then saved for */
505 /* future access. */
506
507 /* Computing 2nd power */
508 d__1 = ae;
509 ck2 = xj2 * .5 * (d__1 * d__1);
510 /* Computing 4th power */
511 d__1 = ae, d__1 *= d__1;
512 ck4 = xj4 * -.375 * (d__1 * d__1);
513 /* Computing 4th power */
514 d__1 = (qo - so) * ae / xkmper, d__1 *= d__1;
515 qoms2t = d__1 * d__1;
516 s = ae * (so / xkmper + 1.);
517
518 /* Recover original mean motion (XNODP) and semimajor axis (AODP) */
519 /* from input elements */
520
521 d__1 = xke / xno;
522 a1 = pow_dd(&d__1, &c_b19);
523 cosio = cos(xincl);
524 /* Computing 2nd power */
525 d__1 = cosio;
526 theta2 = d__1 * d__1;
527 x3thm1 = theta2 * 3. - 1.;
528 /* Computing 2nd power */
529 d__1 = eo;
530 betao2 = 1. - d__1 * d__1;
531 betao = sqrt(betao2);
532 /* Computing 2nd power */
533 d__1 = a1;
534 del1 = ck2 * 1.5 * x3thm1 / (d__1 * d__1 * betao * betao2);
535 ao = a1 * (1. - del1 * (del1 * (del1 * 1.654320987654321 + 1.) +
536 .33333333333333331));
537 /* Computing 2nd power */
538 d__1 = ao;
539 delo = ck2 * 1.5 * x3thm1 / (d__1 * d__1 * betao * betao2);
540 xnodp = xno / (delo + 1.);
541 aodp = ao / (1. - delo);
542
543 /* For perigee below 156 km, the values of S and QOMS2T are */
544 /* altered */
545
546 s4 = s;
547 qoms24 = qoms2t;
548 perige = (aodp * (1. - eo) - ae) * xkmper;
549 if (perige < 156.) {
550 s4 = perige - 78.;
551 if (perige > 98.) {
552 /* Computing 4th power */
553 d__1 = (120. - s4) * ae / xkmper, d__1 *= d__1;
554 qoms24 = d__1 * d__1;
555 s4 = s4 / xkmper + ae;
556 } else {
557 s4 = 20.;
558 }
559 }
560 /* Computing 2nd power */
561 d__1 = aodp;
562 /* Computing 2nd power */
563 d__2 = betao2;
564 pinvsq = 1. / (d__1 * d__1 * (d__2 * d__2));
565 tsi = 1. / (aodp - s4);
566 eta = aodp * eo * tsi;
567 /* Computing 2nd power */
568 d__1 = eta;
569 etasq = d__1 * d__1;
570 eeta = eo * eta;
571 psisq = (d__1 = 1. - etasq, abs(d__1));
572 /* Computing 4th power */
573 d__1 = tsi, d__1 *= d__1;
574 coef = qoms24 * (d__1 * d__1);
575 coef1 = coef / pow_dd(&psisq, &c_b20);
576 c2 = coef1 * xnodp * (aodp * (etasq * 1.5 + 1. + eeta * (etasq + 4.))
577 + ck2 * .75 * tsi / psisq * x3thm1 * (etasq * 3. * (etasq +
578 8.) + 8.));
579 c1 = bstar * c2;
580 sinio = sin(xincl);
581 /* Computing 3rd power */
582 d__1 = ae;
583 a3ovk2 = -xj3 / ck2 * (d__1 * (d__1 * d__1));
584 x1mth2 = 1. - theta2;
585 c4 = xnodp * 2. * coef1 * aodp * betao2 * (eta * (etasq * .5 + 2.) +
586 eo * (etasq * 2. + .5) - ck2 * 2. * tsi / (aodp * psisq) * (
587 x3thm1 * -3. * (1. - eeta * 2. + etasq * (1.5 - eeta * .5)) +
588 x1mth2 * .75 * (etasq * 2. - eeta * (etasq + 1.)) * cos(
589 omegao * 2.)));
590 temp1 = ck2 * 3. * pinvsq * xnodp;
591 temp2 = temp1 * ck2 * pinvsq;
592 temp3 = ck4 * 1.25 * pinvsq * pinvsq * xnodp;
593 xmdot = xnodp + temp1 * .5 * betao * x3thm1 + temp2 * .0625 * betao *
594 (theta2 * (theta2 * 137. - 78.) + 13.);
595 x1m5th = 1. - theta2 * 5.;
596 omgdot = temp1 * -.5 * x1m5th + temp2 * .0625 * (theta2 * (theta2 *
597 395. - 114.) + 7.) + temp3 * (theta2 * (theta2 * 49. - 36.) +
598 3.);
599 xhdot1 = -temp1 * cosio;
600 xnodot = xhdot1 + (temp2 * .5 * (4. - theta2 * 19.) + temp3 * 2. * (
601 3. - theta2 * 7.)) * cosio;
602 xnodcf = betao2 * 3.5 * xhdot1 * c1;
603 t2cof = c1 * 1.5;
604 xlcof = a3ovk2 * .125 * sinio * (cosio * 5. + 3.) / (cosio + 1.);
605 aycof = a3ovk2 * .25 * sinio;
606 x7thm1 = theta2 * 7. - 1.;
607 }
608 zzdpinit_(&aodp, &xmdot, &omgdot, &xnodot, &xnodp, elems);
609
610 /* Get the time since the EPOCH in minutes. */
611
612 tsince = (*time - epoch) / 60.;
613
614 /* Update for secular gravity and atmospheric drag */
615
616 xmdf = xmo + xmdot * tsince;
617 omgadf = omegao + omgdot * tsince;
618 xnoddf = xnodeo + xnodot * tsince;
619 tsq = tsince * tsince;
620 xnode = xnoddf + xnodcf * tsq;
621 tempa = 1. - c1 * tsince;
622 tempe = bstar * c4 * tsince;
623 templ = t2cof * tsq;
624 xn = xnodp;
625
626 /* Calculate the secular terms. */
627
628 zzdpsec_(&xmdf, &omgadf, &xnode, &em, &xinc, &xn, &tsince, elems, &omgdot)
629 ;
630 d__1 = xke / xn;
631 /* Computing 2nd power */
632 d__2 = tempa;
633 a = pow_dd(&d__1, &c_b19) * (d__2 * d__2);
634 e = em - tempe;
635 xmam = xmdf + xnodp * templ;
636
637 /* Calculate the periodic terms. */
638
639 zzdpper_(&tsince, &e, &xinc, &omgadf, &xnode, &xmam);
640 xl = xmam + omgadf + xnode;
641 xn = xke / pow_dd(&a, &c_b22);
642
643 /* Long period periodics */
644
645 axn = e * cos(omgadf);
646 /* Computing 2nd power */
647 d__1 = e;
648 temp = 1. / (a * (1. - d__1 * d__1));
649 xll = temp * xlcof * axn;
650 aynl = temp * aycof;
651 xlt = xl + xll;
652 ayn = e * sin(omgadf) + aynl;
653
654 /* Solve Kepler's equation */
655
656 /* U = EPW - AXN * SIN(EPW) + AYN * COS(EPW) */
657
658 /* Where */
659
660 /* AYN = E * SIN(OMEGA) + AYNL */
661 /* AXN = E * COS(OMEGA) */
662
663 /* And */
664
665 /* AYNL = -0.50D0 * SINIO * AE * J3 / (J2 * A * (1.0D0 - E^2)) */
666
667
668 /* Get the mod division of CAPU with 2 Pi */
669
670 d__1 = xlt - xnode;
671 capu = d_mod(&d__1, &pix2);
672 if (capu < 0.) {
673 capu += pix2;
674 }
675
676 /* Set initial states for the Kepler solution */
677
678 epw = capu;
679 cont = TRUE_;
680 while(cont) {
681 temp2 = epw;
682 sinepw = sin(temp2);
683 cosepw = cos(temp2);
684 temp3 = axn * sinepw;
685 temp4 = ayn * cosepw;
686 temp5 = axn * cosepw;
687 temp6 = ayn * sinepw;
688 epw = (capu - temp4 + temp3 - temp2) / (1. - temp5 - temp6) + temp2;
689
690 /* Test for convergence against the defined tolerance */
691
692 if ((d__1 = epw - temp2, abs(d__1)) <= 1e-6) {
693 cont = FALSE_;
694 }
695 }
696
697 /* Short period preliminary quantities */
698
699 ecose = temp5 + temp6;
700 esine = temp3 - temp4;
701 elsq = axn * axn + ayn * ayn;
702 temp = 1. - elsq;
703 pl = a * temp;
704 rk = a * (1. - ecose);
705 temp1 = 1. / rk;
706 rdot = xke * sqrt(a) * esine * temp1;
707 rfdot = xke * sqrt(pl) * temp1;
708 temp2 = a * temp1;
709 betal = sqrt(temp);
710 temp3 = 1. / (betal + 1.);
711 cosu = temp2 * (cosepw - axn + ayn * esine * temp3);
712 sinu = temp2 * (sinepw - ayn - axn * esine * temp3);
713
714 /* Compute the angle from the x-axis of the point ( COSU, SINU ) */
715
716 if (sinu != 0. || cosu != 0.) {
717 uang = atan2(sinu, cosu);
718 if (uang < 0.) {
719 uang += pix2;
720 }
721 } else {
722 uang = 0.;
723 }
724 sin2u = sinu * 2. * cosu;
725 cos2u = cosu * 2. * cosu - 1.;
726 temp1 = ck2 * (1. / pl);
727 temp2 = temp1 * (1. / pl);
728
729 /* Update for short periodics */
730
731 rk = rk * (1. - temp2 * 1.5 * betal * x3thm1) + temp1 * .5 * x1mth2 *
732 cos2u;
733 uk = uang - temp2 * .25 * x7thm1 * sin2u;
734 xnodek = xnode + temp2 * 1.5 * cosio * sin2u;
735 xinck = xinc + temp2 * 1.5 * cosio * sinio * cos2u;
736 rdotk = rdot - xn * temp1 * x1mth2 * sin2u;
737 rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + x3thm1 * 1.5);
738
739 /* Orientation vectors are calculated by */
740
741 /* U = M sin(uk) + N cos(uk) */
742 /* V = M cos(uk) - N sin(uk) */
743
744 /* Where M and N are euclidean 3 vectors */
745
746 /* M = (-sin(xnodek)cos(xinck), cos(xnodek)cos(xinck), sin(xinck) ) */
747 /* N = ( cos(xnodek), sin(xnodek) , 0 ) */
748
749 sinuk = sin(uk);
750 cosuk = cos(uk);
751
752 /* Use LATREC to generate M and N. M is a latitude to rectangle */
753 /* conversion of a unit vector where PI/2 + XNODEK is the longitude */
754
755 d__1 = pio2 + xnodek;
756 latrec_(&c_b23, &d__1, &xinck, m);
757 latrec_(&c_b23, &xnodek, &c_b25, n);
758
759 /* Sum the components to obtain U and V */
760
761 vlcom_(&sinuk, m, &cosuk, n, u);
762 d__1 = -sinuk;
763 vlcom_(&cosuk, m, &d__1, n, v);
764
765 /* Determine the position and velocity then pack the STATE vector */
766 /* with value scaled to KM and KPS. */
767
768 /* R = RK U + 0 V */
769 /* V = RKDOT U + RK RFDOT V */
770
771 scale = xkmper / ae;
772 d__1 = rk * scale;
773 vlcom_(&d__1, u, &c_b25, v, state);
774
775 /* Now scale to KPS for the velocity component */
776
777 scale /= 60.;
778 d__1 = rdotk * scale;
779 d__2 = rfdotk * scale;
780 vlcom_(&d__1, u, &d__2, v, &state[3]);
781
782 /* All done now.... */
783
784 chkout_("DPSPCE", (ftnlen)6);
785 return 0;
786 } /* dpspce_ */
787
788