1 /* ev2lin.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_b91 = .66666666666666663;
11 static doublereal c_b92 = 3.5;
12 static doublereal c_b153 = 1.5;
13 static integer c__20 = 20;
14
15 /* $Procedure EV2LIN ( Evaluate "two-line" element data) */
ev2lin_(doublereal * et,doublereal * geophs,doublereal * elems,doublereal * state)16 /* Subroutine */ int ev2lin_(doublereal *et, doublereal *geophs, doublereal *
17 elems, doublereal *state)
18 {
19 /* Initialized data */
20
21 static logical doinit = TRUE_;
22
23 /* System generated locals */
24 integer i__1;
25 doublereal d__1;
26
27 /* Builtin functions */
28 integer s_rnge(char *, integer, char *, integer);
29 double cos(doublereal), sin(doublereal), sqrt(doublereal), pow_dd(
30 doublereal *, doublereal *), d_mod(doublereal *, doublereal *),
31 atan2(doublereal, doublereal);
32
33 /* Local variables */
34 static integer head;
35 static doublereal coef, eeta, delm, aodp, delo, capu, xmdf, aynl, elsq,
36 temp;
37 static integer last;
38 static doublereal rdot, cosu, tokm;
39 static integer list[12] /* was [2][6] */;
40 static doublereal sinu, coef1, t2cof, t3cof, t4cof, t5cof, temp1, temp2,
41 temp3, temp4, temp5, cos2u, temp6, mov1m, sin2u, a, e, f;
42 static integer i__, j;
43 static doublereal m;
44 static integer n;
45 static doublereal r__, s, u, betal, omega, betao;
46 extern /* Subroutine */ int chkin_(char *, ftnlen);
47 static doublereal epoch, ecose, aycof, delmo, esine, a3ovk2, tcube, cosik,
48 tempa, bstar, cosio, xincl, etasq, rfdot, sinik, a1, rdotk, c1,
49 c2, c3, c4, c5, cosuk, d2, d3, j2, j3, j4, qomso, d4, lower;
50 extern doublereal twopi_(void);
51 static doublereal q1, q2, psisq, qoms24, s4, sinio, sinmo, sinuk, tempe,
52 betao2, betao3, betao4, templ, tfour, upper, x1m5th, x1mth2,
53 x3thm1, x7thm1, fmod2p, theta2, theta4, xinck, xlcof, xmcof,
54 xmdot, xnode, xnodp;
55 static integer count;
56 static doublereal xndd6o;
57 static integer after;
58 static logical recog, unrec;
59 static doublereal ae, xhdot1;
60 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen);
61 static doublereal xndt2o, ke, ao, fl, eo, qoms2t, er, fu, pl, omgadf, rk,
62 qo, uk, so, xl;
63 static integer before;
64 static doublereal xn, omegao, delomg;
65 extern doublereal brcktd_(doublereal *, doublereal *, doublereal *);
66 static doublereal omgcof, perige, ux, uy, uz, fprime, elemnt[60] /*
67 was [10][6] */, tsince, ae2, ae3, ae4, epsiln, xnodeo, cosnok,
68 lstgeo[8], omgdot, ck2, cosepw, ck4, prelim[174] /* was [29][6]
69 */, rfdotk, sinepw, sinnok, vx, tokmps, vy, pinvsq, vz, xnodcf,
70 xnoddf, xnodek, epwnxt, xnodot;
71 static logical newgeo;
72 extern /* Subroutine */ int setmsg_(char *, ftnlen), errint_(char *,
73 integer *, ftnlen), sigerr_(char *, ftnlen), chkout_(char *,
74 ftnlen);
75 static doublereal eta, axn, ayn, epw, est, tsi, xll, xmo, xno, xmp, tsq,
76 xlt, xmx, xmy, del1, c1sq, pix2;
77
78 /* $ Abstract */
79
80 /* This routine evaluates NORAD two-line element data for */
81 /* near-earth orbiting spacecraft (that is spacecraft with */
82 /* orbital periods less than 225 minutes). */
83
84 /* $ Disclaimer */
85
86 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
87 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
88 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
89 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
90 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
91 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
92 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
93 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
94 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
95 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
96
97 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
98 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
99 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
100 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
101 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
102 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
103
104 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
105 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
106 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
107 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
108
109 /* $ Required_Reading */
110
111 /* None. */
112
113 /* $ Keywords */
114
115 /* EPHEMERIS */
116
117 /* $ Declarations */
118 /* $ Brief_I/O */
119
120 /* VARIABLE I/O DESCRIPTION */
121 /* -------- --- -------------------------------------------------- */
122 /* ET I Epoch in seconds past ephemeris epoch J2000. */
123 /* GEOPHS I Geophysical constants */
124 /* ELEMS I Two-line element data */
125 /* STATE O Evaluated state */
126 /* NMODL P Parameter controlling number of buffered elements. */
127
128 /* $ Detailed_Input */
129
130 /* ET is the poch in seconds past ephemeris epoch J2000 */
131 /* at which a state should be produced from the */
132 /* 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 /* GEOPHS(7) = RE: Equatorial radius of the earth in km. */
155
156
157 /* GEOPHS(8) = AE: Distance units/earth radius */
158 /* (normally 1) */
159
160 /* Below are currently recommended values for these */
161 /* items: */
162
163 /* J2 = 1.082616D-3 */
164 /* J3 = -2.53881D-6 */
165 /* J4 = -1.65597D-6 */
166
167 /* The next item is the square root of GM for the */
168 /* earth given in units of earth-radii**1.5/Minute */
169
170 /* KE = 7.43669161D-2 */
171
172 /* The next two items give the top and */
173 /* bottom of the atmospheric drag model */
174 /* used by the type 10 ephemeris type. */
175 /* Don't adjust these unless you understand */
176 /* the full implications of such changes. */
177
178 /* QO = 120.0D0 */
179 /* SO = 78.0D0 */
180
181 /* The following is the equatorial radius */
182 /* of the earth as used by NORAD in km. */
183
184 /* ER = 6378.135D0 */
185
186 /* The value of AE is the number of */
187 /* distance units per earth radii used by */
188 /* the NORAD state propagation software. */
189 /* The value should be 1 unless you've got */
190 /* a very good understanding of the NORAD */
191 /* routine SGP4 and the affect of changing */
192 /* this value.. */
193
194 /* AE = 1.0D0 */
195
196 /* ELEMS is an array containg two-line element data */
197 /* as prescribed below. The elements XNDD6O and BSTAR */
198 /* must already be scaled by the proper exponent stored */
199 /* in the two line elements set. Moreover, the */
200 /* various items must be converted to the units shown */
201 /* here. */
202
203 /* ELEMS ( 1 ) = XNDT2O in radians/minute**2 */
204 /* ELEMS ( 2 ) = XNDD6O in radians/minute**3 */
205 /* ELEMS ( 3 ) = BSTAR */
206 /* ELEMS ( 4 ) = XINCL in radians */
207 /* ELEMS ( 5 ) = XNODEO in radians */
208 /* ELEMS ( 6 ) = EO */
209 /* ELEMS ( 7 ) = OMEGAO in radians */
210 /* ELEMS ( 8 ) = XMO in radians */
211 /* ELEMS ( 9 ) = XNO in radians/minute */
212 /* ELEMS ( 10 ) = EPOCH of the elements in seconds */
213 /* past ephemeris epoch J2000. */
214
215 /* $ Detailed_Output */
216
217 /* STATE is the state produced by evaluating the input elements */
218 /* at the input epoch ET. Units are km and km/sec. */
219
220 /* $ Parameters */
221
222 /* NMODL is a parameter that controls how many element sets */
223 /* can be buffered internally. Since there are a lot */
224 /* of computations that are independent of time these */
225 /* are buffered and only computed if an unbuffered */
226 /* model is supplied. This value should always */
227 /* be at least 2. Increasing it a great deal is not */
228 /* advised since the time needed to search the */
229 /* buffered elements for a match increases linearly */
230 /* with the NMODL. Imperically, 6 seems to be a good */
231 /* break even value for NMODL. */
232
233 /* $ Exceptions */
234
235 /* 1) No checks are made on the reasonableness of the inputs. */
236
237 /* 2) SPICE(ITERATIONEXCEEDED) signals if the EST calculation loop */
238 /* exceds the MXLOOP value. This error should signal only for */
239 /* bad (nonphysical) TLEs. */
240
241 /* $ Files */
242
243 /* None. */
244
245 /* $ Particulars */
246
247 /* This routine evaluates NORAD two-line element sets for */
248 /* near-earth orbitting satellites. Near earth is defined to */
249 /* be a satellite with an orbital period of less than 225 */
250 /* minutes. This code is an adaptation of the NORAD routine */
251 /* SGP4 to elliminate common blocks, allow buffering of models */
252 /* and intermediate parameters and double precision calculations. */
253
254 /* $ Examples */
255
256 /* None. */
257
258 /* $ Restrictions */
259
260 /* None. */
261
262 /* $ Literature_References */
263
264 /* None. */
265
266 /* $ Author_and_Institution */
267
268 /* W.L. Taber (JPL) */
269
270 /* $ Version */
271
272 /* - SPICELIB Version 1.1.0, 15-SEP-2014 (EDW) */
273
274 /* Added error check to prevent infinite loop in */
275 /* calculation of EST. */
276
277 /* - SPICELIB Version 1.0.3, 02-JAN-2008 (EDW) */
278
279 /* Corrected error in the calculation of the C4 term */
280 /* identified by Curtis Haase. */
281
282 /* Minor edit to the COEF1 declaration strictly */
283 /* identifying the constant as a double. */
284
285 /* From: */
286
287 /* COEF1 = COEF / PSISQ**3.5 */
288
289 /* To: */
290
291 /* COEF1 = COEF / PSISQ**3.5D0 */
292
293 /* - SPICELIB Version 1.0.2, 08-JUL-2004 (EDW) */
294
295 /* Corrected error in the calculation of the C2 term. */
296 /* Reordered C1, C2 calculations to avoid division */
297 /* by BSTAR. */
298
299 /* - SPICELIB Version 1.0.1, 10-MAR-1998 (EDW) */
300
301 /* Corrected error in header describing the GEOPHS array. */
302
303 /* - SPICELIB Version 1.0.0, 14-JAN-1994 (WLT) */
304
305 /* -& */
306 /* $ Index_Entries */
307
308 /* Evaluate NORAD two-line element data. */
309
310 /* -& */
311
312 /* Spicelib functions */
313
314
315 /* Local Parameters */
316
317
318 /* The following parameters give the location of the various */
319 /* geophysical parameters needed for the two line element */
320 /* sets. */
321
322 /* KJ2 --- location of J2 */
323 /* KJ3 --- location of J3 */
324 /* KJ4 --- location if J4 */
325 /* KKE --- location of KE = sqrt(GM) in eart-radii**1.5/MIN */
326 /* KQO --- upper bound of atmospheric model in KM */
327 /* KSO --- lower bound of atmospheric model in KM */
328 /* KER --- earth equatorial radius in KM. */
329 /* KAE --- distance units/earth radius */
330
331
332
333 /* An enumeration of the various components of the */
334 /* elements array---ELEMS */
335
336 /* KNDT20 */
337 /* KNDD60 */
338 /* KBSTAR */
339 /* KINCL */
340 /* KNODE0 */
341 /* KECC */
342 /* KOMEGA */
343 /* KMO */
344 /* KNO */
345
346
347 /* The parameters NEXT and PREV are used in our linked list */
348 /* LIST(NEXT,I) points to the list item the occurs after */
349 /* list item I. LIST ( PREV, I ) points to the list item */
350 /* that preceeds list item I. */
351 /* NEXT */
352 /* PREV */
353
354
355 /* There are a number of preliminary quantities that are needed */
356 /* to compute the state. Those that are not time dependent and */
357 /* depend only upon the elements are stored in a buffer so that */
358 /* if an element set matches a saved set, these preliminary */
359 /* quantities will not be recomputed. PRSIZE is the parameter */
360 /* used to declare the needed room. */
361
362
363 /* When we perform bisection in the solution of Kepler's equation */
364 /* we don't want to bisect too many times. */
365
366
367 /* Numerical Constants */
368
369
370 /* Local variables */
371
372 /* Geophysical Quantities */
373
374
375 /* Elements */
376
377
378 /* Intermediate quantities. The time independent quantities */
379 /* are calculated only as new elements come into the routine. */
380
381 chkin_("EV2LIN", (ftnlen)6);
382
383 /* Rather than always making function calls we store the */
384 /* values of the PI dependent constants the first time */
385 /* through the routine. */
386
387 if (doinit) {
388 doinit = FALSE_;
389 pix2 = twopi_();
390 for (i__ = 1; i__ <= 8; ++i__) {
391 lstgeo[(i__1 = i__ - 1) < 8 && 0 <= i__1 ? i__1 : s_rnge("lstgeo",
392 i__1, "ev2lin_", (ftnlen)605)] = 0.;
393 }
394 for (i__ = 1; i__ <= 6; ++i__) {
395 for (j = 1; j <= 10; ++j) {
396 elemnt[(i__1 = j + i__ * 10 - 11) < 60 && 0 <= i__1 ? i__1 :
397 s_rnge("elemnt", i__1, "ev2lin_", (ftnlen)610)] = 0.;
398 }
399 }
400
401 /* Set up our doubly linked list of most recently used */
402 /* models. Here's how things are supposed to be arranged: */
403
404 /* LIST(NEXT,I) points to the ephemeris model that was used */
405 /* most recently after ephemeris model I. */
406 /* LIST(PREV,I) points to the latest ephemeris model used */
407 /* that was used more recently than I. */
408
409 /* HEAD points to the most recently used ephemris */
410 /* model. */
411
412
413 head = 1;
414 list[(i__1 = (head << 1) - 1) < 12 && 0 <= i__1 ? i__1 : s_rnge("list"
415 , i__1, "ev2lin_", (ftnlen)629)] = 0;
416 list[0] = 2;
417 for (i__ = 2; i__ <= 5; ++i__) {
418 list[(i__1 = (i__ << 1) - 2) < 12 && 0 <= i__1 ? i__1 : s_rnge(
419 "list", i__1, "ev2lin_", (ftnlen)634)] = i__ + 1;
420 list[(i__1 = (i__ << 1) - 1) < 12 && 0 <= i__1 ? i__1 : s_rnge(
421 "list", i__1, "ev2lin_", (ftnlen)635)] = i__ - 1;
422 }
423 list[10] = 0;
424 list[11] = 5;
425 }
426
427 /* We update the geophysical parameters only if there */
428 /* has been a change from the last time they were */
429 /* supplied. */
430
431 if (lstgeo[7] != geophs[7] || lstgeo[6] != geophs[6] || lstgeo[0] !=
432 geophs[0] || lstgeo[1] != geophs[1] || lstgeo[2] != geophs[2] ||
433 lstgeo[3] != geophs[3] || lstgeo[4] != geophs[4] || lstgeo[5] !=
434 geophs[5]) {
435 for (i__ = 1; i__ <= 8; ++i__) {
436 lstgeo[(i__1 = i__ - 1) < 8 && 0 <= i__1 ? i__1 : s_rnge("lstgeo",
437 i__1, "ev2lin_", (ftnlen)657)] = geophs[i__ - 1];
438 }
439 j2 = geophs[0];
440 j3 = geophs[1];
441 j4 = geophs[2];
442 ke = geophs[3];
443 qo = geophs[4];
444 so = geophs[5];
445 er = geophs[6];
446 ae = geophs[7];
447 ae2 = ae * ae;
448 ae3 = ae * ae2;
449 ae4 = ae * ae3;
450 ck2 = j2 * .5 * ae2;
451 a3ovk2 = j3 * -2. * ae / j2;
452 ck4 = j4 * -.375 * ae4;
453 qomso = qo - so;
454 q1 = qomso * ae / er;
455 q2 = q1 * q1;
456 qoms2t = q2 * q2;
457 s = ae * (so / er + 1.);
458
459 /* When we've finished up we will need to convert everything */
460 /* back to KM and KM/SEC the two variables below give the */
461 /* factors we shall need to do this. */
462
463 tokm = er / ae;
464 tokmps = tokm / 60.;
465 newgeo = TRUE_;
466 } else {
467 newgeo = FALSE_;
468 }
469
470 /* Fetch all of the pieces of this model. */
471
472 epoch = elems[9];
473 xndt2o = elems[0];
474 xndd6o = elems[1];
475 bstar = elems[2];
476 xincl = elems[3];
477 xnodeo = elems[4];
478 eo = elems[5];
479 omegao = elems[6];
480 xmo = elems[7];
481 xno = elems[8];
482
483 /* See if this model is already buffered, start at the first */
484 /* model in the list (the most recently used model). */
485
486 unrec = TRUE_;
487 n = head;
488 while(n != 0 && unrec) {
489
490 /* The actual order of the elements is such that we can */
491 /* usually tell that a stored model is different from */
492 /* the one under consideration by looking at the */
493 /* end of the list first. Hence we start with I = NELEMS */
494 /* and decrement I until we have looked at everything */
495 /* or found a mismatch. */
496
497 recog = TRUE_;
498 i__ = 10;
499 while(recog && i__ > 0) {
500 recog = recog && elemnt[(i__1 = i__ + n * 10 - 11) < 60 && 0 <=
501 i__1 ? i__1 : s_rnge("elemnt", i__1, "ev2lin_", (ftnlen)
502 733)] == elems[i__ - 1];
503 --i__;
504 }
505 unrec = ! recog;
506 if (unrec) {
507 last = n;
508 n = list[(i__1 = (n << 1) - 2) < 12 && 0 <= i__1 ? i__1 : s_rnge(
509 "list", i__1, "ev2lin_", (ftnlen)741)];
510 }
511 }
512 if (n == 0) {
513 n = last;
514 }
515
516 /* Either N points to a recognized item or it points to the */
517 /* tail of the list where the least recently used items is */
518 /* located. In either case N must be made the head of the */
519 /* list. (If it is already the head of the list we don't */
520 /* have to bother with anything.) */
521
522 if (n != head) {
523
524 /* Find the items that come before and after N and */
525 /* link them together. */
526
527 before = list[(i__1 = (n << 1) - 1) < 12 && 0 <= i__1 ? i__1 : s_rnge(
528 "list", i__1, "ev2lin_", (ftnlen)762)];
529 after = list[(i__1 = (n << 1) - 2) < 12 && 0 <= i__1 ? i__1 : s_rnge(
530 "list", i__1, "ev2lin_", (ftnlen)763)];
531 list[(i__1 = (before << 1) - 2) < 12 && 0 <= i__1 ? i__1 : s_rnge(
532 "list", i__1, "ev2lin_", (ftnlen)765)] = after;
533 if (after != 0) {
534 list[(i__1 = (after << 1) - 1) < 12 && 0 <= i__1 ? i__1 : s_rnge(
535 "list", i__1, "ev2lin_", (ftnlen)768)] = before;
536 }
537
538 /* Now the guy that will come after N is the current */
539 /* head of the list. N will have no predecessor. */
540
541 list[(i__1 = (n << 1) - 2) < 12 && 0 <= i__1 ? i__1 : s_rnge("list",
542 i__1, "ev2lin_", (ftnlen)774)] = head;
543 list[(i__1 = (n << 1) - 1) < 12 && 0 <= i__1 ? i__1 : s_rnge("list",
544 i__1, "ev2lin_", (ftnlen)775)] = 0;
545
546 /* The predecessor the current head of the list becomes N */
547
548 list[(i__1 = (head << 1) - 1) < 12 && 0 <= i__1 ? i__1 : s_rnge("list"
549 , i__1, "ev2lin_", (ftnlen)779)] = n;
550
551 /* and finally, N becomes the head of the list. */
552
553 head = n;
554 }
555 if (recog && ! newgeo) {
556
557 /* We can just look up the intermediate values from */
558 /* computations performed on a previous call to this */
559 /* routine. */
560
561 aodp = prelim[(i__1 = n * 29 - 29) < 174 && 0 <= i__1 ? i__1 : s_rnge(
562 "prelim", i__1, "ev2lin_", (ftnlen)794)];
563 aycof = prelim[(i__1 = n * 29 - 28) < 174 && 0 <= i__1 ? i__1 :
564 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)795)];
565 c1 = prelim[(i__1 = n * 29 - 27) < 174 && 0 <= i__1 ? i__1 : s_rnge(
566 "prelim", i__1, "ev2lin_", (ftnlen)796)];
567 c4 = prelim[(i__1 = n * 29 - 26) < 174 && 0 <= i__1 ? i__1 : s_rnge(
568 "prelim", i__1, "ev2lin_", (ftnlen)797)];
569 c5 = prelim[(i__1 = n * 29 - 25) < 174 && 0 <= i__1 ? i__1 : s_rnge(
570 "prelim", i__1, "ev2lin_", (ftnlen)798)];
571 cosio = prelim[(i__1 = n * 29 - 24) < 174 && 0 <= i__1 ? i__1 :
572 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)799)];
573 d2 = prelim[(i__1 = n * 29 - 23) < 174 && 0 <= i__1 ? i__1 : s_rnge(
574 "prelim", i__1, "ev2lin_", (ftnlen)800)];
575 d3 = prelim[(i__1 = n * 29 - 22) < 174 && 0 <= i__1 ? i__1 : s_rnge(
576 "prelim", i__1, "ev2lin_", (ftnlen)801)];
577 d4 = prelim[(i__1 = n * 29 - 21) < 174 && 0 <= i__1 ? i__1 : s_rnge(
578 "prelim", i__1, "ev2lin_", (ftnlen)802)];
579 delmo = prelim[(i__1 = n * 29 - 20) < 174 && 0 <= i__1 ? i__1 :
580 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)803)];
581 eta = prelim[(i__1 = n * 29 - 19) < 174 && 0 <= i__1 ? i__1 : s_rnge(
582 "prelim", i__1, "ev2lin_", (ftnlen)804)];
583 omgcof = prelim[(i__1 = n * 29 - 18) < 174 && 0 <= i__1 ? i__1 :
584 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)805)];
585 omgdot = prelim[(i__1 = n * 29 - 17) < 174 && 0 <= i__1 ? i__1 :
586 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)806)];
587 perige = prelim[(i__1 = n * 29 - 16) < 174 && 0 <= i__1 ? i__1 :
588 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)807)];
589 sinio = prelim[(i__1 = n * 29 - 15) < 174 && 0 <= i__1 ? i__1 :
590 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)808)];
591 sinmo = prelim[(i__1 = n * 29 - 14) < 174 && 0 <= i__1 ? i__1 :
592 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)809)];
593 t2cof = prelim[(i__1 = n * 29 - 13) < 174 && 0 <= i__1 ? i__1 :
594 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)810)];
595 t3cof = prelim[(i__1 = n * 29 - 12) < 174 && 0 <= i__1 ? i__1 :
596 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)811)];
597 t4cof = prelim[(i__1 = n * 29 - 11) < 174 && 0 <= i__1 ? i__1 :
598 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)812)];
599 t5cof = prelim[(i__1 = n * 29 - 10) < 174 && 0 <= i__1 ? i__1 :
600 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)813)];
601 x1mth2 = prelim[(i__1 = n * 29 - 9) < 174 && 0 <= i__1 ? i__1 :
602 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)814)];
603 x3thm1 = prelim[(i__1 = n * 29 - 8) < 174 && 0 <= i__1 ? i__1 :
604 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)815)];
605 x7thm1 = prelim[(i__1 = n * 29 - 7) < 174 && 0 <= i__1 ? i__1 :
606 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)816)];
607 xlcof = prelim[(i__1 = n * 29 - 6) < 174 && 0 <= i__1 ? i__1 : s_rnge(
608 "prelim", i__1, "ev2lin_", (ftnlen)817)];
609 xmcof = prelim[(i__1 = n * 29 - 5) < 174 && 0 <= i__1 ? i__1 : s_rnge(
610 "prelim", i__1, "ev2lin_", (ftnlen)818)];
611 xmdot = prelim[(i__1 = n * 29 - 4) < 174 && 0 <= i__1 ? i__1 : s_rnge(
612 "prelim", i__1, "ev2lin_", (ftnlen)819)];
613 xnodcf = prelim[(i__1 = n * 29 - 3) < 174 && 0 <= i__1 ? i__1 :
614 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)820)];
615 xnodot = prelim[(i__1 = n * 29 - 2) < 174 && 0 <= i__1 ? i__1 :
616 s_rnge("prelim", i__1, "ev2lin_", (ftnlen)821)];
617 xnodp = prelim[(i__1 = n * 29 - 1) < 174 && 0 <= i__1 ? i__1 : s_rnge(
618 "prelim", i__1, "ev2lin_", (ftnlen)822)];
619 } else {
620
621 /* Compute all of the intermediate items needed. */
622 /* First, the inclination dependent constants. */
623
624 cosio = cos(xincl);
625 sinio = sin(xincl);
626 theta2 = cosio * cosio;
627 theta4 = theta2 * theta2;
628 x3thm1 = theta2 * 3. - 1.;
629 x7thm1 = theta2 * 7. - 1.;
630 x1mth2 = 1. - theta2;
631 x1m5th = 1. - theta2 * 5.;
632
633 /* Eccentricity dependent constants */
634
635 betao = sqrt(1. - eo * eo);
636 betao2 = 1. - eo * eo;
637 betao3 = betao * betao2;
638 betao4 = betao2 * betao2;
639
640 /* Semi-major axis and ascending node related constants. */
641
642 d__1 = ke / xno;
643 a1 = pow_dd(&d__1, &c_b91);
644 del1 = ck2 * 1.5 * x3thm1 / (a1 * a1 * betao3);
645 ao = a1 * (1. - del1 * (del1 * (del1 * 134. / 81. + 1.) +
646 .33333333333333331));
647 delo = ck2 * 1.5 * x3thm1 / (ao * ao * betao3);
648 xnodp = xno / (delo + 1.);
649 aodp = ao / (1. - delo);
650 s4 = s;
651 qoms24 = qoms2t;
652 perige = er * (aodp * (1. - eo) - ae);
653
654 /* For perigee below 156 km, the values of S and QOMS2T are */
655 /* altered. */
656
657 if (perige < 156.) {
658 s4 = perige - 78.;
659 if (perige <= 98.) {
660 s4 = 20.;
661 }
662 /* Computing 4th power */
663 d__1 = (120. - s4) * ae / er, d__1 *= d__1;
664 qoms24 = d__1 * d__1;
665 s4 = ae + s4 / er;
666 }
667
668 /* The next block is simply a pretty print of the code in */
669 /* sgp4 from label number 10 through the label 90. */
670
671 pinvsq = 1. / (aodp * aodp * betao4);
672 tsi = 1. / (aodp - s4);
673 eta = aodp * eo * tsi;
674 etasq = eta * eta;
675 eeta = eo * eta;
676 /* Computing 4th power */
677 d__1 = tsi, d__1 *= d__1;
678 coef = qoms24 * (d__1 * d__1);
679 psisq = (d__1 = 1. - etasq, abs(d__1));
680 coef1 = coef / pow_dd(&psisq, &c_b92);
681 c2 = coef1 * xnodp * (aodp * (etasq * 1.5 + 1. + eeta * (etasq + 4.))
682 + ck2 * .75 * (tsi / psisq) * x3thm1 * (etasq * (etasq * 3. +
683 24.) + 8.));
684 c1 = c2 * bstar;
685 c3 = coef * tsi * a3ovk2 * xnodp * ae * sinio / eo;
686 c4 = xnodp * 2. * coef1 * aodp * betao2 * (eta * (etasq * .5 + 2.) +
687 eo * (etasq * 2. + .5) - ck2 * tsi / (aodp * psisq) * 2. * (
688 x3thm1 * -3. * (1. - eeta * 2. + etasq * (1.5 - eeta * .5)) +
689 cos(omegao * 2.) * .75 * x1mth2 * (etasq * 2. - eeta * (etasq
690 + 1.))));
691 c5 = coef1 * 2. * aodp * betao2 * ((etasq + eeta) * 2.75 + 1. + eeta *
692 etasq);
693 temp1 = ck2 * 3. * pinvsq * xnodp;
694 temp2 = temp1 * ck2 * pinvsq;
695 temp3 = ck4 * 1.25 * pinvsq * pinvsq * xnodp;
696 xmdot = xnodp + temp1 * .5 * betao * x3thm1 + temp2 * .0625 * betao *
697 (13. - theta2 * 78. + theta4 * 137.);
698 omgdot = temp1 * -.5 * x1m5th + temp2 * .0625 * (7. - theta2 * 114. +
699 theta4 * 395.) + temp3 * (3. - theta2 * 36. + theta4 * 49.);
700 xhdot1 = -temp1 * cosio;
701 xnodot = xhdot1 + cosio * (temp2 * .5 * (4. - theta2 * 19.) + temp3 *
702 2. * (3. - theta2 * 7.));
703 omgcof = bstar * c3 * cos(omegao);
704 xmcof = -bstar * .66666666666666663 * coef * ae / eeta;
705 xnodcf = betao2 * 3.5 * xhdot1 * c1;
706 t2cof = c1 * 1.5;
707 aycof = a3ovk2 * .25 * sinio;
708 xlcof = aycof * .5 * (cosio * 5. + 3.) / (cosio + 1.);
709 /* Computing 3rd power */
710 d__1 = eta * cos(xmo) + 1.;
711 delmo = d__1 * (d__1 * d__1);
712 sinmo = sin(xmo);
713
714 /* For perigee less than 220 kilometers, the ISIMP flag is set */
715 /* and the equations are truncated to linear variation in SQRT */
716 /* A and quadratic variation in mean anomaly. Also, the C3 */
717 /* term, the Delta OMEGA term, and the Delta M term are */
718 /* dropped. (Note: Normally we would just use */
719
720 if (perige >= 220.) {
721 c1sq = c1 * c1;
722 d2 = tsi * 4. * c1sq * aodp;
723 temp = d2 * tsi * c1 * .33333333333333331;
724 d3 = temp * (s4 + aodp * 17.);
725 d4 = temp * tsi * c1 * aodp * .5 * (aodp * 221. + s4 * 31.);
726 t3cof = d2 + c1sq * 2.;
727 t4cof = (d3 * 3. + c1 * (d2 * 12. + c1sq * 10.)) * .25;
728 t5cof = (d4 * 3. + c1 * 12. * d3 + d2 * 6. * d2 + c1sq * 15. * (
729 d2 * 2. + c1sq)) * .2;
730 }
731
732 /* Now store the intermediate computations so that if we */
733 /* should hit this model again we can just look up the needed */
734 /* results from the above computations. */
735
736 prelim[(i__1 = n * 29 - 29) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
737 "im", i__1, "ev2lin_", (ftnlen)992)] = aodp;
738 prelim[(i__1 = n * 29 - 28) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
739 "im", i__1, "ev2lin_", (ftnlen)993)] = aycof;
740 prelim[(i__1 = n * 29 - 27) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
741 "im", i__1, "ev2lin_", (ftnlen)994)] = c1;
742 prelim[(i__1 = n * 29 - 26) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
743 "im", i__1, "ev2lin_", (ftnlen)995)] = c4;
744 prelim[(i__1 = n * 29 - 25) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
745 "im", i__1, "ev2lin_", (ftnlen)996)] = c5;
746 prelim[(i__1 = n * 29 - 24) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
747 "im", i__1, "ev2lin_", (ftnlen)997)] = cosio;
748 prelim[(i__1 = n * 29 - 23) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
749 "im", i__1, "ev2lin_", (ftnlen)998)] = d2;
750 prelim[(i__1 = n * 29 - 22) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
751 "im", i__1, "ev2lin_", (ftnlen)999)] = d3;
752 prelim[(i__1 = n * 29 - 21) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
753 "im", i__1, "ev2lin_", (ftnlen)1000)] = d4;
754 prelim[(i__1 = n * 29 - 20) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
755 "im", i__1, "ev2lin_", (ftnlen)1001)] = delmo;
756 prelim[(i__1 = n * 29 - 19) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
757 "im", i__1, "ev2lin_", (ftnlen)1002)] = eta;
758 prelim[(i__1 = n * 29 - 18) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
759 "im", i__1, "ev2lin_", (ftnlen)1003)] = omgcof;
760 prelim[(i__1 = n * 29 - 17) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
761 "im", i__1, "ev2lin_", (ftnlen)1004)] = omgdot;
762 prelim[(i__1 = n * 29 - 16) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
763 "im", i__1, "ev2lin_", (ftnlen)1005)] = perige;
764 prelim[(i__1 = n * 29 - 15) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
765 "im", i__1, "ev2lin_", (ftnlen)1006)] = sinio;
766 prelim[(i__1 = n * 29 - 14) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
767 "im", i__1, "ev2lin_", (ftnlen)1007)] = sinmo;
768 prelim[(i__1 = n * 29 - 13) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
769 "im", i__1, "ev2lin_", (ftnlen)1008)] = t2cof;
770 prelim[(i__1 = n * 29 - 12) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
771 "im", i__1, "ev2lin_", (ftnlen)1009)] = t3cof;
772 prelim[(i__1 = n * 29 - 11) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
773 "im", i__1, "ev2lin_", (ftnlen)1010)] = t4cof;
774 prelim[(i__1 = n * 29 - 10) < 174 && 0 <= i__1 ? i__1 : s_rnge("prel"
775 "im", i__1, "ev2lin_", (ftnlen)1011)] = t5cof;
776 prelim[(i__1 = n * 29 - 9) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
777 , i__1, "ev2lin_", (ftnlen)1012)] = x1mth2;
778 prelim[(i__1 = n * 29 - 8) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
779 , i__1, "ev2lin_", (ftnlen)1013)] = x3thm1;
780 prelim[(i__1 = n * 29 - 7) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
781 , i__1, "ev2lin_", (ftnlen)1014)] = x7thm1;
782 prelim[(i__1 = n * 29 - 6) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
783 , i__1, "ev2lin_", (ftnlen)1015)] = xlcof;
784 prelim[(i__1 = n * 29 - 5) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
785 , i__1, "ev2lin_", (ftnlen)1016)] = xmcof;
786 prelim[(i__1 = n * 29 - 4) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
787 , i__1, "ev2lin_", (ftnlen)1017)] = xmdot;
788 prelim[(i__1 = n * 29 - 3) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
789 , i__1, "ev2lin_", (ftnlen)1018)] = xnodcf;
790 prelim[(i__1 = n * 29 - 2) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
791 , i__1, "ev2lin_", (ftnlen)1019)] = xnodot;
792 prelim[(i__1 = n * 29 - 1) < 174 && 0 <= i__1 ? i__1 : s_rnge("prelim"
793 , i__1, "ev2lin_", (ftnlen)1020)] = xnodp;
794
795 /* Finally, move these elements in the storage area */
796 /* for checking the next time through. */
797
798 for (i__ = 1; i__ <= 10; ++i__) {
799 elemnt[(i__1 = i__ + n * 10 - 11) < 60 && 0 <= i__1 ? i__1 :
800 s_rnge("elemnt", i__1, "ev2lin_", (ftnlen)1026)] = elems[
801 i__ - 1];
802 }
803 }
804
805 /* Now that all of the introductions are out of the way */
806 /* we can get down to business. */
807
808 /* Compute the time since the epoch for this model. */
809
810 tsince = *et - epoch;
811
812 /* and convert it to minutes */
813
814 tsince /= 60.;
815 xmdf = xmo + xmdot * tsince;
816 omgadf = omegao + omgdot * tsince;
817 xnoddf = xnodeo + xnodot * tsince;
818 omega = omgadf;
819 xmp = xmdf;
820 tsq = tsince * tsince;
821 xnode = xnoddf + xnodcf * tsq;
822 tempa = 1. - c1 * tsince;
823 tempe = bstar * c4 * tsince;
824 templ = t2cof * tsq;
825 if (perige > 220.) {
826 tcube = tsq * tsince;
827 tfour = tcube * tsince;
828 delomg = omgcof * tsince;
829 /* Computing 3rd power */
830 d__1 = eta * cos(xmdf) + 1.;
831 delm = xmcof * (d__1 * (d__1 * d__1) - delmo);
832 temp = delomg + delm;
833 xmp = xmdf + temp;
834 omega = omgadf - temp;
835 tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour;
836 tempe += bstar * c5 * (sin(xmp) - sinmo);
837 templ = templ + tcube * t3cof + tfour * (t4cof + tsince * t5cof);
838 }
839 /* Computing 2nd power */
840 d__1 = tempa;
841 a = aodp * (d__1 * d__1);
842 xl = xmp + omega + xnode + xnodp * templ;
843 e = eo - tempe;
844
845 /* The parameter BETA used to be needed, but it's only use */
846 /* was in the computation of TEMP below where it got squared */
847 /* so we'll remove it from the list of things to compute. */
848
849 /* BETA = DSQRT( 1.0D0 - E*E ) */
850
851 xn = ke / pow_dd(&a, &c_b153);
852
853 /* Long period periodics */
854
855 temp = 1. / (a * (1. - e * e));
856 aynl = temp * aycof;
857 ayn = e * sin(omega) + aynl;
858 axn = e * cos(omega);
859 xll = temp * xlcof * axn;
860 xlt = xl + xll;
861
862 /* Solve keplers equation. */
863
864 /* We are going to solve for the roots of this equation by */
865 /* using a mixture of Newton's method and the prescription for */
866 /* root finding outlined in the SPICE routine UNITIM. */
867
868 /* We are going to solve the equation */
869
870 /* U = EPW - AXN * SIN(EPW) + AYN * COS(EPW) */
871
872 /* Where */
873
874 /* AYN = E * SIN(OMEGA) + AYNL */
875 /* AXN = E * COS(OMEGA) */
876
877 /* And */
878
879 /* AYNL = -0.50D0 * SINIO * AE * J3 / (J2*A*(1.0D0 - E*E)) */
880
881 /* Since this is a low earth orbiter (period less than 225 minutes) */
882 /* The maximum value E can take (without having the orbiter */
883 /* plowing fields) is approximately 0.47 and AYNL will not be */
884 /* more than about .01. ( Typically E will be much smaller */
885 /* on the order of about .1 ) Thus we can initially */
886 /* view the problem of solving the equation for EPW as a */
887 /* function of the form */
888
889 /* U = EPW + F ( EPW ) (1) */
890
891 /* Where F( EPW ) = -AXN*SIN(EPW) + AYN*COS(EPW) */
892
893 /* Note that | F'(EPW) | < M = DSQRT( AXN**2 + AYN**2 ) < 0.48 */
894
895 /* From the above discussion it is evident that F is a contraction */
896 /* mapping. So that we can employ the same techniques as were */
897 /* used in the routine UNITIM to get our first approximations of */
898 /* the root. Once we have some good first approximations, we */
899 /* will speed up the root finding by using Newton's method for */
900 /* finding a zero of a function. The function we will work on */
901 /* is */
902
903 /* f (x) = x - U - AXN*SIN(x) + AYN*COS(x) (2) */
904
905 /* By applying Newton's method we will go from linear to */
906 /* quadratic convergence. */
907
908 /* We will keep track of our error bounds along the way so */
909 /* that we will know how many iterations to perform in each */
910 /* phase of the root extraction. */
911
912 /* few steps using bisection. */
913
914
915 /* For the benefit of those interested */
916 /* here's the basics of what we'll do. */
917
918 /* Whichever EPW satisfies equation (1) will be */
919 /* unique. The uniqueness of the solution is ensured because the */
920 /* expression on the right-hand side of the equation is */
921 /* monotone increasing in EPW. */
922
923 /* Let's suppose that EPW is the solution, then the following */
924 /* is true. */
925
926 /* EPW = U - F(EPW) */
927
928 /* but we can also replace the EPW on the right hand side of the */
929 /* equation by U - F(EPW). Thus */
930
931 /* EPW = U - F( U - F(EPW)) */
932
933 /* = U - F( U - F( U - F(EPW))) */
934
935 /* = U - F( U - F( U - F( U - F(EPW)))) */
936
937 /* = U - F( U - F( U - F( U - F( U - F(EPW))))) */
938 /* . */
939 /* . */
940 /* . */
941 /* = U - F( U - F( U - F( U - F( U - F(U - ... ))) */
942
943 /* and so on, for as long as we have patience to perform the */
944 /* substitutions. */
945
946 /* The point of doing this recursive substitution is that we */
947 /* hope to move EPW to an insignificant part of the computation. */
948 /* This would seem to have a reasonable chance of success since */
949 /* F is a bounded and has a small derivative. */
950
951 /* Following this idea, we will attempt to solve for EPW using */
952 /* the recursive method outlined below. */
953
954 /* We will make our first guess at EPW, call it EPW_0. */
955
956 /* EPW_0 = U */
957
958 /* Our next guess, EPW_1, is given by: */
959
960 /* EPW_1 = U - F(EPW_0) */
961
962 /* And so on: */
963
964 /* EPW_2 = U - F(EPW_1) [ = U - F(U - F(U)) ] */
965 /* EPW_3 = U - F(EPW_2) [ = U - F(U - F(U - F(U))) ] */
966 /* . */
967 /* . */
968 /* . */
969 /* EPW_n = U - F(EPW_(n-1)) [ = U - F(U - F(U - F(U...)))] */
970
971 /* The questions to ask at this point are: */
972
973 /* 1) Do the EPW_i's converge? */
974 /* 2) If they converge, do they converge to EPW? */
975 /* 3) If they converge to EPW, how fast do they get there? */
976
977 /* 1) The sequence of approximations converges. */
978
979 /* | EPW_n - EPW_(n-1) | = [ U - F( EPW_(n-1) ) ] */
980 /* - [ U - F( EPW_(n-2) ) ] */
981
982 /* = [ F( EPW_(n-2) ) - F( EPW_(n-1)) ] */
983
984 /* The function F has an important property. The absolute */
985 /* value of its derivative is always less than M. */
986 /* This means that for any pair of real numbers s,t */
987
988 /* | F(t) - F(s) | < M*| t - s |. */
989
990 /* From this observation, we can see that */
991
992 /* | EPW_n - EPW_(n-1) | < M*| EPW_(n-1) - EPW_(n-2) | */
993
994 /* With this fact available, we could (with a bit more work) */
995 /* conclude that the sequence of EPW_i's converges and that */
996 /* it converges at a rate that is at least as fast as the */
997 /* sequence M, M**2, M**3. In fact the difference */
998 /* |EPW - EPW_N| < M/(1-M) * | EPW_N - EPW_(N-1) | */
999
1000 /* < M/(1-M) * M**N | EPW_1 - EPW_0 | */
1001
1002 /* 2) If we let EPW be the limit of the EPW_i's then it follows */
1003 /* that */
1004
1005 /* EPW = U - F(EPW). */
1006
1007
1008 /* or that */
1009
1010 /* U = EPW + F(EPW). */
1011
1012 /* We will use this technique to get an approximation that */
1013 /* is within a tolerance of EPW and then switch to */
1014 /* a Newton's method. (We'll compute the tolerance using */
1015 /* the value of M given above). */
1016
1017
1018 /* For the Newton's method portion of the problem, recall */
1019 /* from Taylor's formula that: */
1020
1021 /* f(x) = f(x_0) + f'(x_0)(x-x_0) + f''(c)/2 (x-x_0)**2 */
1022
1023 /* for some c between x and x_0 */
1024
1025 /* If x happens to be a zero of f then we can rearrange the */
1026 /* terms above to get */
1027
1028 /* f(x_0) f''(c) */
1029 /* x = x_0 - ------- + -------- ( x - x_0)**2 */
1030 /* f'(x_0) f'(x_0) */
1031
1032 /* Thus the error in the Newton approximation */
1033
1034
1035 /* f(x_0) */
1036 /* x = x_0 - ------- */
1037 /* f'(x_0) */
1038
1039 /* is */
1040
1041 /* f''(c) */
1042 /* -------- ( x - x_0)**2 */
1043 /* f'(x_0) */
1044
1045 /* Thus if we can bound f'' and pick a good first */
1046 /* choice for x_0 (using the first method outlined */
1047 /* above we can get quadratic convergence.) */
1048
1049 /* In our case we have */
1050
1051 /* f (x) = x - U - AXN*SIN(x) + AYN*COS(x) */
1052 /* f' (x) = 1 - AXN*COS(x) - AYN*SIN(x) */
1053 /* f''(x) = AXN*SIN(x) - AYN*COS(x) */
1054
1055 /* So that: */
1056
1057 /* f' (x) > 1 - M */
1058
1059 /* f''(x) < M */
1060
1061 /* Thus the error in the Newton's approximation is */
1062 /* at most */
1063
1064 /* M/(1-M) * ( x - x_0 )**2 */
1065
1066 /* Thus as long as our original estimate (determined using */
1067 /* the contraction method) gets within a reasonable tolerance */
1068 /* of x, we can use Newton's method to acheive faster */
1069 /* convergence. */
1070
1071 m = sqrt(axn * axn + ayn * ayn);
1072 mov1m = (d__1 = m / (1. - m), abs(d__1));
1073 d__1 = xlt - xnode;
1074 fmod2p = d_mod(&d__1, &pix2);
1075 if (fmod2p < 0.) {
1076 fmod2p += pix2;
1077 }
1078 capu = fmod2p;
1079 epw = capu;
1080 est = 1.;
1081 count = 0;
1082 while(est > .125) {
1083 ++count;
1084 if (count > 20) {
1085 setmsg_("EST iteration count of #1 exceeded at time ET #2. This "
1086 "error may indicate a bad TLE set.", (ftnlen)88);
1087 errint_("#1", &c__20, (ftnlen)2);
1088 errdp_("#2", et, (ftnlen)2);
1089 sigerr_("SPICE(ITERATIONEXCEEDED)", (ftnlen)24);
1090 chkout_("EV2LIN", (ftnlen)6);
1091 return 0;
1092 }
1093 epwnxt = capu - axn * sin(epw) + ayn * cos(epw);
1094 est = mov1m * (d__1 = epwnxt - epw, abs(d__1));
1095 epw = epwnxt;
1096 }
1097
1098 /* We need to be able to add something to EPW and not */
1099 /* get EPW (but not too much). */
1100
1101 epsiln = est;
1102 if (epsiln + epw != epw) {
1103
1104 /* Now we switch over to Newton's method. Note that */
1105 /* since our error estimate is less than 1/8, six iterations */
1106 /* of Newton's method should get us to within 1/2**96 of */
1107 /* the correct answer (If there were no round off to contend */
1108 /* with). */
1109
1110 for (i__ = 1; i__ <= 5; ++i__) {
1111 sinepw = sin(epw);
1112 cosepw = cos(epw);
1113 f = epw - capu - axn * sinepw + ayn * cosepw;
1114 fprime = 1. - axn * cosepw - ayn * sinepw;
1115 epwnxt = epw - f / fprime;
1116
1117 /* Our new error estimate comes from the discussion */
1118 /* of convergence of Newton's method. */
1119
1120 epw = epwnxt;
1121 if (epw + est != epw) {
1122 epsiln = est;
1123 est = mov1m * est * est;
1124 }
1125 }
1126 }
1127
1128 /* Finally, we use bisection to avoid the problems of */
1129 /* round-off that may be present in Newton's method. Since */
1130 /* we've gotten quite close to the answer (theoretically */
1131 /* anyway) we won't have to perform many bisection passes. */
1132
1133 /* First we must bracket the root. Note that we will */
1134 /* increase EPSILN so that we don't spend much time */
1135 /* determining the bracketing interval. Also if the first */
1136 /* addition of EPSILN to EPW doesn't modify it, were set up */
1137 /* to just quit. This happens only if F is sufficiently */
1138 /* close to zero that it can't alter EPW by adding it to */
1139 /* or subtracting it from EPW. */
1140
1141 sinepw = sin(epw);
1142 cosepw = cos(epw);
1143 f = epw - capu - axn * sinepw + ayn * cosepw;
1144 /* Computing MAX */
1145 d__1 = abs(f);
1146 epsiln = max(d__1,epsiln);
1147 if (f == 0.) {
1148 lower = epw;
1149 upper = epw;
1150 } else if (f > 0.) {
1151 fu = f;
1152 upper = epw;
1153 lower = epw - epsiln;
1154 epw = lower;
1155 while(f > 0. && lower != upper) {
1156 epw -= epsiln;
1157 f = epw - capu - axn * sin(epw) + ayn * cos(epw);
1158 epsiln *= 2.;
1159 }
1160 lower = epw;
1161 fl = f;
1162 if (f == 0.) {
1163 upper = lower;
1164 }
1165 } else if (f < 0.) {
1166 fl = f;
1167 lower = epw;
1168 upper = epw + epsiln;
1169 epw = upper;
1170 while(f < 0. && lower != upper) {
1171 epw += epsiln;
1172 f = epw - capu - axn * sin(epw) + ayn * cos(epw);
1173 epsiln *= 2.;
1174 }
1175 upper = epw;
1176 fu = f;
1177 if (f == 0.) {
1178 lower = epw;
1179 }
1180 }
1181
1182 /* Finally, bisect until we can do no more. */
1183
1184 count = 0;
1185 while(upper > lower && count < 20) {
1186 ++count;
1187 d__1 = (upper + lower) * .5;
1188 epw = brcktd_(&d__1, &lower, &upper);
1189
1190 /* EPW eventually will not be different from one of the */
1191 /* two bracketing values. If this is the time, we need */
1192 /* to decide on a value for EPW. That's done below. */
1193
1194 if (epw == upper || epw == lower) {
1195 if (-fl < fu) {
1196 epw = lower;
1197 upper = lower;
1198 } else {
1199 epw = upper;
1200 lower = upper;
1201 }
1202 } else {
1203 f = epw - capu - axn * sin(epw) + ayn * cos(epw);
1204 if (f > 0.) {
1205 upper = epw;
1206 fu = f;
1207 } else if (f < 0.) {
1208 lower = epw;
1209 fl = f;
1210 } else {
1211 lower = epw;
1212 upper = epw;
1213 }
1214 }
1215 }
1216
1217 /* Short period preliminary quantities */
1218
1219 sinepw = sin(epw);
1220 cosepw = cos(epw);
1221 temp3 = axn * sinepw;
1222 temp4 = ayn * cosepw;
1223 temp5 = axn * cosepw;
1224 temp6 = ayn * sinepw;
1225 ecose = temp5 + temp6;
1226 esine = temp3 - temp4;
1227 elsq = axn * axn + ayn * ayn;
1228 temp = 1. - elsq;
1229 pl = a * temp;
1230 r__ = a * (1. - ecose);
1231 temp1 = 1. / r__;
1232 rdot = ke * temp1 * sqrt(a) * esine;
1233 rfdot = ke * temp1 * sqrt(pl);
1234 temp2 = a * temp1;
1235 betal = sqrt(temp);
1236 temp3 = 1. / (betal + 1.);
1237 cosu = temp2 * (cosepw - axn + ayn * esine * temp3);
1238 sinu = temp2 * (sinepw - ayn - axn * esine * temp3);
1239
1240 /* Compute the angle from the x-axis of the point ( COSU, SINU ) */
1241
1242 if (sinu != 0. || cosu != 0.) {
1243 u = atan2(sinu, cosu);
1244 if (u < 0.) {
1245 u += pix2;
1246 }
1247 } else {
1248 u = 0.;
1249 }
1250 sin2u = sinu * 2. * cosu;
1251 cos2u = cosu * 2. * cosu - 1.;
1252 temp = 1. / pl;
1253 temp1 = ck2 * temp;
1254 temp2 = temp1 * temp;
1255
1256 /* Update for short periodics */
1257
1258 rk = r__ * (1. - temp2 * 1.5 * betal * x3thm1) + temp1 * .5 * x1mth2 *
1259 cos2u;
1260 uk = u - temp2 * .25 * x7thm1 * sin2u;
1261 xnodek = xnode + temp2 * 1.5 * cosio * sin2u;
1262 xinck = xincl + temp2 * 1.5 * cosio * cos2u * sinio;
1263 rdotk = rdot - xn * temp1 * x1mth2 * sin2u;
1264 rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + x3thm1 * 1.5);
1265
1266 /* Orientation vectors */
1267
1268 sinuk = sin(uk);
1269 cosuk = cos(uk);
1270 sinik = sin(xinck);
1271 cosik = cos(xinck);
1272 sinnok = sin(xnodek);
1273 cosnok = cos(xnodek);
1274 xmx = -sinnok * cosik;
1275 xmy = cosnok * cosik;
1276 ux = xmx * sinuk + cosnok * cosuk;
1277 uy = xmy * sinuk + sinnok * cosuk;
1278 uz = sinik * sinuk;
1279 vx = xmx * cosuk - cosnok * sinuk;
1280 vy = xmy * cosuk - sinnok * sinuk;
1281 vz = sinik * cosuk;
1282
1283 /* Position and velocity */
1284
1285 state[0] = tokm * rk * ux;
1286 state[1] = tokm * rk * uy;
1287 state[2] = tokm * rk * uz;
1288 state[3] = tokmps * (rdotk * ux + rfdotk * vx);
1289 state[4] = tokmps * (rdotk * uy + rfdotk * vy);
1290 state[5] = tokmps * (rdotk * uz + rfdotk * vz);
1291 chkout_("EV2LIN", (ftnlen)6);
1292 return 0;
1293 } /* ev2lin_ */
1294
1295