001 /*
002 * $Id: SGP4unit.java,v 1.3 2012/02/19 17:35:39 davep Exp $
003 *
004 * This file is part of McIDAS-V
005 *
006 * Copyright 2007-2012
007 * Space Science and Engineering Center (SSEC)
008 * University of Wisconsin - Madison
009 * 1225 W. Dayton Street, Madison, WI 53706, USA
010 * https://www.ssec.wisc.edu/mcidas
011 *
012 * All Rights Reserved
013 *
014 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
015 * some McIDAS-V source code is based on IDV and VisAD source code.
016 *
017 * McIDAS-V is free software; you can redistribute it and/or modify
018 * it under the terms of the GNU Lesser Public License as published by
019 * the Free Software Foundation; either version 3 of the License, or
020 * (at your option) any later version.
021 *
022 * McIDAS-V is distributed in the hope that it will be useful,
023 * but WITHOUT ANY WARRANTY; without even the implied warranty of
024 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
025 * GNU Lesser Public License for more details.
026 *
027 * You should have received a copy of the GNU Lesser Public License
028 * along with this program. If not, see http://www.gnu.org/licenses.
029 */
030 /**
031 * Vallado/CSSI SGP4 propogator
032 * Converted to Java from C++ by: Shawn E. Gano, 19 June 2009
033 * The goal of this conversion was to stick as closely to the orginal code
034 * as possible and not trying to convert it to OO design.
035 * This code has been compared using the verification TLEs and "SGP4verification.java"
036 * and the results were virtually identical to the c++ version.
037 * Only rare differences were found in the results on the order of 0.001 cm.
038 * ----------------------------------------------------------------
039 *
040 * sgp4unit.cpp
041 *
042 * this file contains the sgp4 procedures for analytical propagation
043 * of a satellite. the code was originally released in the 1980 and 1986
044 * spacetrack papers. a detailed discussion of the theory and history
045 * may be found in the 2006 aiaa paper by vallado, crawford, hujsak,
046 * and kelso.
047 *
048 * companion code for
049 * fundamentals of astrodynamics and applications
050 * 2007
051 * by david vallado
052 *
053 * (w) 719-573-2600, email dvallado@agi.com
054 *
055 * current :
056 * 3 Nov 08 david vallado
057 * put returns in for error codes
058 * changes :
059 * 29 sep 08 david vallado
060 * fix atime for faster operation in dspace
061 * add operationmode for afspc (a) or improved (i)
062 * performance mode
063 * 16 jun 08 david vallado
064 * update small eccentricity check
065 * 16 nov 07 david vallado
066 * misc fixes for better compliance
067 * 20 apr 07 david vallado
068 * misc fixes for constants
069 * 11 aug 06 david vallado
070 * chg lyddane choice back to strn3, constants, misc doc
071 * 15 dec 05 david vallado
072 * misc fixes
073 * 26 jul 05 david vallado
074 * fixes for paper
075 * note that each fix is preceded by a
076 * comment with "sgp4fix" and an explanation of
077 * what was changed
078 * 10 aug 04 david vallado
079 * 2nd printing baseline working
080 * 14 may 01 david vallado
081 * 2nd edition baseline
082 * 80 norad
083 * original baseline
084 * ---------------------------------------------------------------- */
085
086 package edu.wisc.ssec.mcidasv.data.adde.sgp4;
087
088 /**
089 *
090 * @author Shawn E. Gano, shawn@gano.name
091 */
092 public class SGP4unit
093 {
094 public static double pi = 3.14159265358979323846;
095
096 public static enum Gravconsttype
097 {
098 wgs72old,
099 wgs72,
100 wgs84
101 }
102
103 /** -----------------------------------------------------------------------------
104 * Java version outputs double array: [ep,inclp,nodep,argpp,mp]
105 * -------------------------------------------------------------------------------
106 *
107 * procedure dpper
108 *
109 * this procedure provides deep space long period periodic contributions
110 * to the mean elements. by design, these periodics are zero at epoch.
111 * this used to be dscom which included initialization, but it's really a
112 * recurring function.
113 *
114 * author : david vallado 719-573-2600 28 jun 2005
115 *
116 * inputs :
117 * e3 -
118 * ee2 -
119 * peo -
120 * pgho -
121 * pho -
122 * pinco -
123 * plo -
124 * se2 , se3 , sgh2, sgh3, sgh4, sh2, sh3, si2, si3, sl2, sl3, sl4 -
125 * t -
126 * xh2, xh3, xi2, xi3, xl2, xl3, xl4 -
127 * zmol -
128 * zmos -
129 * ep - eccentricity 0.0 - 1.0
130 * inclo - inclination - needed for lyddane modification
131 * nodep - right ascension of ascending node
132 * argpp - argument of perigee
133 * mp - mean anomaly
134 *
135 * outputs :
136 * ep - eccentricity 0.0 - 1.0
137 * inclp - inclination
138 * nodep - right ascension of ascending node
139 * argpp - argument of perigee
140 * mp - mean anomaly
141 *
142 * locals :
143 * alfdp -
144 * betdp -
145 * cosip , sinip , cosop , sinop ,
146 * dalf -
147 * dbet -
148 * dls -
149 * f2, f3 -
150 * pe -
151 * pgh -
152 * ph -
153 * pinc -
154 * pl -
155 * sel , ses , sghl , sghs , shl , shs , sil , sinzf , sis ,
156 * sll , sls
157 * xls -
158 * xnoh -
159 * zf -
160 * zm -
161 *
162 * coupling :
163 * none.
164 *
165 * references :
166 * hoots, roehrich, norad spacetrack report #3 1980
167 * hoots, norad spacetrack report #6 1986
168 * hoots, schumacher and glover 2004
169 * vallado, crawford, hujsak, kelso 2006
170 ----------------------------------------------------------------------------*/
171 // outputs an array with values for, all outputs are also inputs
172 // [ep,inclp,nodep,argpp,mp]
173 private static double[] dpper(
174 double e3, double ee2, double peo, double pgho, double pho,
175 double pinco, double plo, double se2, double se3, double sgh2,
176 double sgh3, double sgh4, double sh2, double sh3, double si2,
177 double si3, double sl2, double sl3, double sl4, double t,
178 double xgh2, double xgh3, double xgh4, double xh2, double xh3,
179 double xi2, double xi3, double xl2, double xl3, double xl4,
180 double zmol, double zmos, double inclo,
181 char init,
182 double ep, double inclp, double nodep, double argpp, double mp,
183 char opsmode)
184 {
185 // return variables -- all also inputs
186 //double inclp,nodep,argpp,mp; // ep -- input and output
187
188 /* --------------------- local variables ------------------------ */
189 final double twopi = 2.0 * pi;
190 double alfdp, betdp, cosip, cosop, dalf, dbet, dls,
191 f2, f3, pe, pgh, ph, pinc, pl,
192 sel, ses, sghl, sghs, shll, shs, sil,
193 sinip, sinop, sinzf, sis, sll, sls, xls,
194 xnoh, zf, zm, zel, zes, znl, zns;
195
196 /* ---------------------- constants ----------------------------- */
197 zns = 1.19459e-5;
198 zes = 0.01675;
199 znl = 1.5835218e-4;
200 zel = 0.05490;
201
202 /* --------------- calculate time varying periodics ----------- */
203 zm = zmos + zns * t;
204 // be sure that the initial call has time set to zero
205 if(init == 'y')
206 {
207 zm = zmos;
208 }
209 zf = zm + 2.0 * zes * Math.sin(zm);
210 sinzf = Math.sin(zf);
211 f2 = 0.5 * sinzf * sinzf - 0.25;
212 f3 = -0.5 * sinzf * Math.cos(zf);
213 ses = se2 * f2 + se3 * f3;
214 sis = si2 * f2 + si3 * f3;
215 sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
216 sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
217 shs = sh2 * f2 + sh3 * f3;
218 zm = zmol + znl * t;
219 if(init == 'y')
220 {
221 zm = zmol;
222 }
223 zf = zm + 2.0 * zel * Math.sin(zm);
224 sinzf = Math.sin(zf);
225 f2 = 0.5 * sinzf * sinzf - 0.25;
226 f3 = -0.5 * sinzf * Math.cos(zf);
227 sel = ee2 * f2 + e3 * f3;
228 sil = xi2 * f2 + xi3 * f3;
229 sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf;
230 sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
231 shll = xh2 * f2 + xh3 * f3;
232 pe = ses + sel;
233 pinc = sis + sil;
234 pl = sls + sll;
235 pgh = sghs + sghl;
236 ph = shs + shll;
237
238 if(init == 'n')
239 {
240 pe = pe - peo;
241 pinc = pinc - pinco;
242 pl = pl - plo;
243 pgh = pgh - pgho;
244 ph = ph - pho;
245 inclp = inclp + pinc;
246 ep = ep + pe;
247 sinip = Math.sin(inclp);
248 cosip = Math.cos(inclp);
249
250 /* ----------------- apply periodics directly ------------ */
251 // sgp4fix for lyddane choice
252 // strn3 used original inclination - this is technically feasible
253 // gsfc used perturbed inclination - also technically feasible
254 // probably best to readjust the 0.2 limit value and limit discontinuity
255 // 0.2 rad = 11.45916 deg
256 // use next line for original strn3 approach and original inclination
257 // if (inclo >= 0.2)
258 // use next line for gsfc version and perturbed inclination
259 if(inclp >= 0.2)
260 {
261 ph = ph / sinip;
262 pgh = pgh - cosip * ph;
263 argpp = argpp + pgh;
264 nodep = nodep + ph;
265 mp = mp + pl;
266 }
267 else
268 {
269 /* ---- apply periodics with lyddane modification ---- */
270 sinop = Math.sin(nodep);
271 cosop = Math.cos(nodep);
272 alfdp = sinip * sinop;
273 betdp = sinip * cosop;
274 dalf = ph * cosop + pinc * cosip * sinop;
275 dbet = -ph * sinop + pinc * cosip * cosop;
276 alfdp = alfdp + dalf;
277 betdp = betdp + dbet;
278 nodep = (nodep % twopi);
279 // sgp4fix for afspc written intrinsic functions
280 // nodep used without a trigonometric function ahead
281 if((nodep < 0.0) && (opsmode == 'a'))
282 {
283 nodep = nodep + twopi;
284 }
285 xls = mp + argpp + cosip * nodep;
286 dls = pl + pgh - pinc * nodep * sinip;
287 xls = xls + dls;
288 xnoh = nodep;
289 nodep = Math.atan2(alfdp, betdp);
290 // sgp4fix for afspc written intrinsic functions
291 // nodep used without a trigonometric function ahead
292 if((nodep < 0.0) && (opsmode == 'a'))
293 {
294 nodep = nodep + twopi;
295 }
296 if(Math.abs(xnoh - nodep) > pi)
297 {
298 if(nodep < xnoh)
299 {
300 nodep = nodep + twopi;
301 }
302 else
303 {
304 nodep = nodep - twopi;
305 }
306 }
307 mp = mp + pl;
308 argpp = xls - mp - cosip * nodep;
309 }
310 } // if init == 'n'
311
312 return new double[]
313 {
314 ep, inclp, nodep, argpp, mp
315 };
316 //#include "debug1.cpp"
317 } // end dpper
318
319 /*-----------------------------------------------------------------------------
320 *
321 * procedure dscom
322 *
323 * this procedure provides deep space common items used by both the secular
324 * and periodics subroutines. input is provided as shown. this routine
325 * used to be called dpper, but the functions inside weren't well organized.
326 *
327 * author : david vallado 719-573-2600 28 jun 2005
328 *
329 * inputs :
330 * epoch -
331 * ep - eccentricity
332 * argpp - argument of perigee
333 * tc -
334 * inclp - inclination
335 * nodep - right ascension of ascending node
336 * np - mean motion
337 *
338 * outputs :
339 * sinim , cosim , sinomm , cosomm , snodm , cnodm
340 * day -
341 * e3 -
342 * ee2 -
343 * em - eccentricity
344 * emsq - eccentricity squared
345 * gam -
346 * peo -
347 * pgho -
348 * pho -
349 * pinco -
350 * plo -
351 * rtemsq -
352 * se2, se3 -
353 * sgh2, sgh3, sgh4 -
354 * sh2, sh3, si2, si3, sl2, sl3, sl4 -
355 * s1, s2, s3, s4, s5, s6, s7 -
356 * ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3 -
357 * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 -
358 * xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4 -
359 * nm - mean motion
360 * z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 -
361 * zmol -
362 * zmos -
363 *
364 * locals :
365 * a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 -
366 * betasq -
367 * cc -
368 * ctem, stem -
369 * x1, x2, x3, x4, x5, x6, x7, x8 -
370 * xnodce -
371 * xnoi -
372 * zcosg , zsing , zcosgl , zsingl , zcosh , zsinh , zcoshl , zsinhl ,
373 * zcosi , zsini , zcosil , zsinil ,
374 * zx -
375 * zy -
376 *
377 * coupling :
378 * none.
379 *
380 * references :
381 * hoots, roehrich, norad spacetrack report #3 1980
382 * hoots, norad spacetrack report #6 1986
383 * hoots, schumacher and glover 2004
384 * vallado, crawford, hujsak, kelso 2006
385 * ----------------------------------------------------------------------------
386 * constructor modified to return an array with all the values that are not contained in the SGP4SatData object
387 * Array returned with these values:
388 * [snodm, cnodm, sinim, cosim, sinomm, cosomm, day, em, emsq, gam, rtemsq, s1, s2, s3,
389 * s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3, sz11, sz12, sz13,
390 * sz21, sz22, sz23, sz31, sz32, sz33, nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 ]
391 * --------------------------------------------------------------*/
392 private static double[] dscom(
393 // inputs
394 double epoch, double ep, double argpp, double tc, double inclp,
395 double nodep, double np,
396 SGP4SatData satrec // // not
397 // double& snodm, double& cnodm, double& sinim, double& cosim, double& sinomm,
398 // double& cosomm,double& day,
399 // //SGP4SatData
400 // double& e3, double& ee2,
401 // // not
402 // double& em, double& emsq, double& gam,
403 // // SGP4SatData
404 // double& peo, double& pgho, double& pho,
405 // double& pinco, double& plo,
406 // // not
407 // double& rtemsq,
408 // // SGP4SatData
409 // double& se2, double& se3,
410 // double& sgh2, double& sgh3, double& sgh4, double& sh2, double& sh3,
411 // double& si2, double& si3, double& sl2, double& sl3, double& sl4,
412 // // not
413 // double& s1, double& s2, double& s3, double& s4, double& s5,
414 // double& s6, double& s7, double& ss1, double& ss2, double& ss3,
415 // double& ss4, double& ss5, double& ss6, double& ss7, double& sz1,
416 // double& sz2, double& sz3, double& sz11, double& sz12, double& sz13,
417 // double& sz21, double& sz22, double& sz23, double& sz31, double& sz32,
418 // double& sz33,
419 // // SGP4SatData
420 // double& xgh2, double& xgh3, double& xgh4, double& xh2,
421 // double& xh3, double& xi2, double& xi3, double& xl2, double& xl3,
422 // double& xl4,
423 // // not
424 // double& nm, double& z1, double& z2, double& z3,
425 // double& z11, double& z12, double& z13, double& z21, double& z22,
426 // double& z23, double& z31, double& z32, double& z33,
427 // // SGP4SatData
428 // double& zmol, double& zmos
429 )
430 {
431 // Return variables not in SGP4SatData ----------------------------
432 double snodm, cnodm, sinim, cosim, sinomm, cosomm, day, em, emsq, gam, rtemsq,
433 s1, s2, s3, s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7,
434 sz1, sz2, sz3, sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33,
435 nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33;
436
437 // SEG : it is okay to initalize these here as they are not assigned before this function is called
438 s1 = 0;
439 s2 = 0;
440 s3 = 0;
441 s4 = 0;
442 s5 = 0;
443 s6 = 0;
444 s7 = 0;
445 ss1 = 0;
446 ss2 = 0;
447 ss3 = 0;
448 ss4 = 0;
449 ss5 = 0;
450 ss6 = 0;
451 ss7 = 0;
452 sz1 = 0;
453 sz2 = 0;
454 sz3 = 0;
455 sz11 = 0;
456 sz12 = 0;
457 sz13 = 0;
458 sz21 = 0;
459 sz22 = 0;
460 sz23 = 0;
461 sz31 = 0;
462 sz32 = 0;
463 sz33 = 0;
464 z1 = 0;
465 z2 = 0;
466 z3 = 0;
467 z11 = 0;
468 z12 = 0;
469 z13 = 0;
470 z21 = 0;
471 z22 = 0;
472 z23 = 0;
473 z31 = 0;
474 z32 = 0;
475 z33 = 0;
476
477 /* -------------------------- constants ------------------------- */
478 final double zes = 0.01675;
479 final double zel = 0.05490;
480 final double c1ss = 2.9864797e-6;
481 final double c1l = 4.7968065e-7;
482 final double zsinis = 0.39785416;
483 final double zcosis = 0.91744867;
484 final double zcosgs = 0.1945905;
485 final double zsings = -0.98088458;
486 final double twopi = 2.0 * pi;
487
488 /* --------------------- local variables ------------------------ */
489 int lsflg;
490 double a1, a2, a3, a4, a5, a6, a7,
491 a8, a9, a10, betasq, cc, ctem, stem,
492 x1, x2, x3, x4, x5, x6, x7,
493 x8, xnodce, xnoi, zcosg, zcosgl, zcosh, zcoshl,
494 zcosi, zcosil, zsing, zsingl, zsinh, zsinhl, zsini,
495 zsinil, zx, zy;
496
497 nm = np;
498 em = ep;
499 snodm = Math.sin(nodep);
500 cnodm = Math.cos(nodep);
501 sinomm = Math.sin(argpp);
502 cosomm = Math.cos(argpp);
503 sinim = Math.sin(inclp);
504 cosim = Math.cos(inclp);
505 emsq = em * em;
506 betasq = 1.0 - emsq;
507 rtemsq = Math.sqrt(betasq);
508
509 /* ----------------- initialize lunar solar terms --------------- */
510 satrec.peo = 0.0;
511 satrec.pinco = 0.0;
512 satrec.plo = 0.0;
513 satrec.pgho = 0.0;
514 satrec.pho = 0.0;
515 day = epoch + 18261.5 + tc / 1440.0;
516 xnodce = ((4.5236020 - 9.2422029e-4 * day) % twopi);
517 stem = Math.sin(xnodce);
518 ctem = Math.cos(xnodce);
519 zcosil = 0.91375164 - 0.03568096 * ctem;
520 zsinil = Math.sqrt(1.0 - zcosil * zcosil);
521 zsinhl = 0.089683511 * stem / zsinil;
522 zcoshl = Math.sqrt(1.0 - zsinhl * zsinhl);
523 gam = 5.8351514 + 0.0019443680 * day;
524 zx = 0.39785416 * stem / zsinil;
525 zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
526 zx = Math.atan2(zx, zy);
527 zx = gam + zx - xnodce;
528 zcosgl = Math.cos(zx);
529 zsingl = Math.sin(zx);
530
531 /* ------------------------- do solar terms --------------------- */
532 zcosg = zcosgs;
533 zsing = zsings;
534 zcosi = zcosis;
535 zsini = zsinis;
536 zcosh = cnodm;
537 zsinh = snodm;
538 cc = c1ss;
539 xnoi = 1.0 / nm;
540
541 for(lsflg = 1; lsflg <= 2; lsflg++)
542 {
543 a1 = zcosg * zcosh + zsing * zcosi * zsinh;
544 a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
545 a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
546 a8 = zsing * zsini;
547 a9 = zsing * zsinh + zcosg * zcosi * zcosh;
548 a10 = zcosg * zsini;
549 a2 = cosim * a7 + sinim * a8;
550 a4 = cosim * a9 + sinim * a10;
551 a5 = -sinim * a7 + cosim * a8;
552 a6 = -sinim * a9 + cosim * a10;
553
554 x1 = a1 * cosomm + a2 * sinomm;
555 x2 = a3 * cosomm + a4 * sinomm;
556 x3 = -a1 * sinomm + a2 * cosomm;
557 x4 = -a3 * sinomm + a4 * cosomm;
558 x5 = a5 * sinomm;
559 x6 = a6 * sinomm;
560 x7 = a5 * cosomm;
561 x8 = a6 * cosomm;
562
563 z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
564 z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
565 z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
566 z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq;
567 z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq;
568 z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq;
569 z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5);
570 z12 = -6.0 * (a1 * a6 + a3 * a5) + emsq *
571 (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5));
572 z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6);
573 z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7);
574 z22 = 6.0 * (a4 * a5 + a2 * a6) + emsq *
575 (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8));
576 z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8);
577 z1 = z1 + z1 + betasq * z31;
578 z2 = z2 + z2 + betasq * z32;
579 z3 = z3 + z3 + betasq * z33;
580 s3 = cc * xnoi;
581 s2 = -0.5 * s3 / rtemsq;
582 s4 = s3 * rtemsq;
583 s1 = -15.0 * em * s4;
584 s5 = x1 * x3 + x2 * x4;
585 s6 = x2 * x3 + x1 * x4;
586 s7 = x2 * x4 - x1 * x3;
587
588 /* ----------------------- do lunar terms ------------------- */
589 if(lsflg == 1)
590 {
591 ss1 = s1;
592 ss2 = s2;
593 ss3 = s3;
594 ss4 = s4;
595 ss5 = s5;
596 ss6 = s6;
597 ss7 = s7;
598 sz1 = z1;
599 sz2 = z2;
600 sz3 = z3;
601 sz11 = z11;
602 sz12 = z12;
603 sz13 = z13;
604 sz21 = z21;
605 sz22 = z22;
606 sz23 = z23;
607 sz31 = z31;
608 sz32 = z32;
609 sz33 = z33;
610 zcosg = zcosgl;
611 zsing = zsingl;
612 zcosi = zcosil;
613 zsini = zsinil;
614 zcosh = zcoshl * cnodm + zsinhl * snodm;
615 zsinh = snodm * zcoshl - cnodm * zsinhl;
616 cc = c1l;
617 }
618 }
619
620 satrec.zmol = ((4.7199672 + 0.22997150 * day - gam) % twopi);
621 satrec.zmos = ((6.2565837 + 0.017201977 * day) % twopi);
622
623 /* ------------------------ do solar terms ---------------------- */
624 satrec.se2 = 2.0 * ss1 * ss6;
625 satrec.se3 = 2.0 * ss1 * ss7;
626 satrec.si2 = 2.0 * ss2 * sz12;
627 satrec.si3 = 2.0 * ss2 * (sz13 - sz11);
628 satrec.sl2 = -2.0 * ss3 * sz2;
629 satrec.sl3 = -2.0 * ss3 * (sz3 - sz1);
630 satrec.sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes;
631 satrec.sgh2 = 2.0 * ss4 * sz32;
632 satrec.sgh3 = 2.0 * ss4 * (sz33 - sz31);
633 satrec.sgh4 = -18.0 * ss4 * zes;
634 satrec.sh2 = -2.0 * ss2 * sz22;
635 satrec.sh3 = -2.0 * ss2 * (sz23 - sz21);
636
637 /* ------------------------ do lunar terms ---------------------- */
638 satrec.ee2 = 2.0 * s1 * s6;
639 satrec.e3 = 2.0 * s1 * s7;
640 satrec.xi2 = 2.0 * s2 * z12;
641 satrec.xi3 = 2.0 * s2 * (z13 - z11);
642 satrec.xl2 = -2.0 * s3 * z2;
643 satrec.xl3 = -2.0 * s3 * (z3 - z1);
644 satrec.xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel;
645 satrec.xgh2 = 2.0 * s4 * z32;
646 satrec.xgh3 = 2.0 * s4 * (z33 - z31);
647 satrec.xgh4 = -18.0 * s4 * zel;
648 satrec.xh2 = -2.0 * s2 * z22;
649 satrec.xh3 = -2.0 * s2 * (z23 - z21);
650
651 return new double[]
652 {
653 snodm, cnodm, sinim, cosim, sinomm, cosomm, day, em, emsq, gam, rtemsq, s1, s2, s3,
654 s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3, sz11, sz12, sz13,
655 sz21, sz22, sz23, sz31, sz32, sz33, nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33
656 };
657
658 //#include "debug2.cpp"
659 } // end dscom
660
661 /*-----------------------------------------------------------------------------
662 * Java version returns a double array with these values: [ em, argpm, inclm, mm, nm, nodem, dndt]
663 * -----------------------------------------------------------------------------
664 * procedure dsinit
665 *
666 * this procedure provides deep space contributions to mean motion dot due
667 * to geopotential resonance with half day and one day orbits.
668 *
669 * author : david vallado 719-573-2600 28 jun 2005
670 *
671 * inputs :
672 * cosim, sinim-
673 * emsq - eccentricity squared
674 * argpo - argument of perigee
675 * s1, s2, s3, s4, s5 -
676 * ss1, ss2, ss3, ss4, ss5 -
677 * sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33 -
678 * t - time
679 * tc -
680 * gsto - greenwich sidereal time rad
681 * mo - mean anomaly
682 * mdot - mean anomaly dot (rate)
683 * no - mean motion
684 * nodeo - right ascension of ascending node
685 * nodedot - right ascension of ascending node dot (rate)
686 * xpidot -
687 * z1, z3, z11, z13, z21, z23, z31, z33 -
688 * eccm - eccentricity
689 * argpm - argument of perigee
690 * inclm - inclination
691 * mm - mean anomaly
692 * xn - mean motion
693 * nodem - right ascension of ascending node
694 *
695 * outputs :
696 * em - eccentricity
697 * argpm - argument of perigee
698 * inclm - inclination
699 * mm - mean anomaly
700 * nm - mean motion
701 * nodem - right ascension of ascending node
702 * irez - flag for resonance 0-none, 1-one day, 2-half day
703 * atime -
704 * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 -
705 * dedt -
706 * didt -
707 * dmdt -
708 * dndt -
709 * dnodt -
710 * domdt -
711 * del1, del2, del3 -
712 * ses , sghl , sghs , sgs , shl , shs , sis , sls
713 * theta -
714 * xfact -
715 * xlamo -
716 * xli -
717 * xni
718 *
719 * locals :
720 * ainv2 -
721 * aonv -
722 * cosisq -
723 * eoc -
724 * f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543 -
725 * g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533 -
726 * sini2 -
727 * temp -
728 * temp1 -
729 * theta -
730 * xno2 -
731 *
732 * coupling :
733 * getgravconst
734 *
735 * references :
736 * hoots, roehrich, norad spacetrack report #3 1980
737 * hoots, norad spacetrack report #6 1986
738 * hoots, schumacher and glover 2004
739 * vallado, crawford, hujsak, kelso 2006
740 ----------------------------------------------------------------------------*/
741 // returns array containing:
742 // [ em, argpm, inclm, mm, nm, nodem, dndt]
743 private static double[] dsinit(
744 Gravconsttype whichconst,
745 double cosim, double emsq, double argpo, double s1, double s2,
746 double s3, double s4, double s5, double sinim, double ss1,
747 double ss2, double ss3, double ss4, double ss5, double sz1,
748 double sz3, double sz11, double sz13, double sz21, double sz23,
749 double sz31, double sz33, double t, double tc, double gsto,
750 double mo, double mdot, double no, double nodeo, double nodedot,
751 double xpidot, double z1, double z3, double z11, double z13,
752 double z21, double z23, double z31, double z33, double ecco,
753 double eccsq,
754 // return
755 SGP4SatData satrec,
756 double em, // variable is also an INPUT and output - SEG!!
757 double argpm,
758 double inclm,
759 double mm,
760 double nm, // variable is also an INPUT and output - SEG!!
761 double nodem // // not in SGP4SatData
762 // double& em, double& argpm, double& inclm, double& mm,
763 // double& nm, double& nodem,
764 // // in
765 // int& irez,
766 // double& atime, double& d2201, double& d2211, double& d3210, double& d3222,
767 // double& d4410, double& d4422, double& d5220, double& d5232, double& d5421,
768 // double& d5433, double& dedt, double& didt, double& dmdt,
769 // // not
770 // double& dndt,
771 // // in
772 // double& dnodt, double& domdt, double& del1, double& del2, double& del3,
773 // double& xfact, double& xlamo, double& xli, double& xni
774 )
775 {
776 // return values ------------------
777 double dndt; // em,nm, inclm, argpm, nodem, mm are also inputs
778
779 /* --------------------- local variables ------------------------ */
780 final double twopi = 2.0 * pi;
781
782 double ainv2, aonv = 0.0, cosisq, eoc, f220, f221, f311,
783 f321, f322, f330, f441, f442, f522, f523,
784 f542, f543, g200, g201, g211, g300, g310,
785 g322, g410, g422, g520, g521, g532, g533,
786 ses, sgs, sghl, sghs, shs, shll, sis,
787 sini2, sls, temp, temp1, theta, xno2, q22,
788 q31, q33, root22, root44, root54, rptim, root32,
789 root52, x2o3, xke, znl, emo, zns, emsqo,
790 tumin, mu, radiusearthkm, j2, j3, j4, j3oj2;
791
792 q22 = 1.7891679e-6;
793 q31 = 2.1460748e-6;
794 q33 = 2.2123015e-7;
795 root22 = 1.7891679e-6;
796 root44 = 7.3636953e-9;
797 root54 = 2.1765803e-9;
798 rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec
799 root32 = 3.7393792e-7;
800 root52 = 1.1428639e-7;
801 x2o3 = 2.0 / 3.0;
802 znl = 1.5835218e-4;
803 zns = 1.19459e-5;
804
805 // sgp4fix identify constants and allow alternate values
806 double[] temp2 = getgravconst(whichconst);
807 tumin = temp2[0];
808 mu = temp2[1];
809 radiusearthkm = temp2[2];
810 xke = temp2[3];
811 j2 = temp2[4];
812 j3 = temp2[5];
813 j4 = temp2[6];
814 j3oj2 = temp2[7];
815
816 /* -------------------- deep space initialization ------------ */
817 satrec.irez = 0;
818 if((nm < 0.0052359877) && (nm > 0.0034906585))
819 {
820 satrec.irez = 1;
821 }
822 if((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5))
823 {
824 satrec.irez = 2;
825 }
826
827 /* ------------------------ do solar terms ------------------- */
828 ses = ss1 * zns * ss5;
829 sis = ss2 * zns * (sz11 + sz13);
830 sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
831 sghs = ss4 * zns * (sz31 + sz33 - 6.0);
832 shs = -zns * ss2 * (sz21 + sz23);
833 // sgp4fix for 180 deg incl
834 if((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2))
835 {
836 shs = 0.0;
837 }
838 if(sinim != 0.0)
839 {
840 shs = shs / sinim;
841 }
842 sgs = sghs - cosim * shs;
843
844 /* ------------------------- do lunar terms ------------------ */
845 satrec.dedt = ses + s1 * znl * s5;
846 satrec.didt = sis + s2 * znl * (z11 + z13);
847 satrec.dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
848 sghl = s4 * znl * (z31 + z33 - 6.0);
849 shll = -znl * s2 * (z21 + z23);
850 // sgp4fix for 180 deg incl
851 if((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2))
852 {
853 shll = 0.0;
854 }
855 satrec.domdt = sgs + sghl;
856 satrec.dnodt = shs;
857 if(sinim != 0.0)
858 {
859 satrec.domdt = satrec.domdt - cosim / sinim * shll;
860 satrec.dnodt = satrec.dnodt + shll / sinim;
861 }
862
863 /* ----------- calculate deep space resonance effects -------- */
864 dndt = 0.0;
865 theta = ((gsto + tc * rptim) % twopi);
866 em = em + satrec.dedt * t;
867 inclm = inclm + satrec.didt * t;
868 argpm = argpm + satrec.domdt * t;
869 nodem = nodem + satrec.dnodt * t;
870 mm = mm + satrec.dmdt * t;
871 // sgp4fix for negative inclinations
872 // the following if statement should be commented out
873 //if (inclm < 0.0)
874 // {
875 // inclm = -inclm;
876 // argpm = argpm - pi;
877 // nodem = nodem + pi;
878 // }
879
880 /* -------------- initialize the resonance terms ------------- */
881 if(satrec.irez != 0)
882 {
883 aonv = Math.pow(nm / xke, x2o3);
884
885 /* ---------- geopotential resonance for 12 hour orbits ------ */
886 if(satrec.irez == 2)
887 {
888 cosisq = cosim * cosim;
889 emo = em;
890 em = ecco;
891 emsqo = emsq;
892 emsq = eccsq;
893 eoc = em * emsq;
894 g201 = -0.306 - (em - 0.64) * 0.440;
895
896 if(em <= 0.65)
897 {
898 g211 = 3.616 - 13.2470 * em + 16.2900 * emsq;
899 g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc;
900 g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc;
901 g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc;
902 g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
903 g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;
904 }
905 else
906 {
907 g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc;
908 g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc;
909 g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc;
910 g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc;
911 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
912 if(em > 0.715)
913 {
914 g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
915 }
916 else
917 {
918 g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq;
919 }
920 }
921 if(em < 0.7)
922 {
923 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc;
924 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
925 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc;
926 }
927 else
928 {
929 g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
930 g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
931 g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;
932 }
933
934 sini2 = sinim * sinim;
935 f220 = 0.75 * (1.0 + 2.0 * cosim + cosisq);
936 f221 = 1.5 * sini2;
937 f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq);
938 f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq);
939 f441 = 35.0 * sini2 * f220;
940 f442 = 39.3750 * sini2 * sini2;
941 f522 = 9.84375 * sinim * (sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) +
942 0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq));
943 f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim +
944 10.0 * cosisq) + 6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq));
945 f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim + cosisq *
946 (-12.0 + 8.0 * cosim + 10.0 * cosisq));
947 f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim + cosisq *
948 (12.0 + 8.0 * cosim - 10.0 * cosisq));
949 xno2 = nm * nm;
950 ainv2 = aonv * aonv;
951 temp1 = 3.0 * xno2 * ainv2;
952 temp = temp1 * root22;
953 satrec.d2201 = temp * f220 * g201;
954 satrec.d2211 = temp * f221 * g211;
955 temp1 = temp1 * aonv;
956 temp = temp1 * root32;
957 satrec.d3210 = temp * f321 * g310;
958 satrec.d3222 = temp * f322 * g322;
959 temp1 = temp1 * aonv;
960 temp = 2.0 * temp1 * root44;
961 satrec.d4410 = temp * f441 * g410;
962 satrec.d4422 = temp * f442 * g422;
963 temp1 = temp1 * aonv;
964 temp = temp1 * root52;
965 satrec.d5220 = temp * f522 * g520;
966 satrec.d5232 = temp * f523 * g532;
967 temp = 2.0 * temp1 * root54;
968 satrec.d5421 = temp * f542 * g521;
969 satrec.d5433 = temp * f543 * g533;
970 satrec.xlamo = ((mo + nodeo + nodeo - theta - theta) % twopi);
971 satrec.xfact = mdot + satrec.dmdt + 2.0 * (nodedot + satrec.dnodt - rptim) - no;
972 em = emo;
973 emsq = emsqo;
974 }
975
976 /* ---------------- synchronous resonance terms -------------- */
977 if(satrec.irez == 1)
978 {
979 g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq);
980 g310 = 1.0 + 2.0 * emsq;
981 g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq);
982 f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim);
983 f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim);
984 f330 = 1.0 + cosim;
985 f330 = 1.875 * f330 * f330 * f330;
986 satrec.del1 = 3.0 * nm * nm * aonv * aonv;
987 satrec.del2 = 2.0 * satrec.del1 * f220 * g200 * q22;
988 satrec.del3 = 3.0 * satrec.del1 * f330 * g300 * q33 * aonv;
989 satrec.del1 = satrec.del1 * f311 * g310 * q31 * aonv;
990 satrec.xlamo = ((mo + nodeo + argpo - theta) % twopi);
991 satrec.xfact = mdot + xpidot - rptim + satrec.dmdt + satrec.domdt + satrec.dnodt - no;
992 }
993
994 /* ------------ for sgp4, initialize the integrator ---------- */
995 satrec.xli = satrec.xlamo;
996 satrec.xni = no;
997 satrec.atime = 0.0;
998 nm = no + dndt;
999 }
1000
1001 return new double[]
1002 {
1003 em, argpm, inclm, mm, nm, nodem, dndt
1004 };
1005
1006 //#include "debug3.cpp"
1007 } // end dsinit
1008
1009 /*-----------------------------------------------------------------------------
1010 *
1011 * procedure dspace
1012 *
1013 * this procedure provides deep space contributions to mean elements for
1014 * perturbing third body. these effects have been averaged over one
1015 * revolution of the sun and moon. for earth resonance effects, the
1016 * effects have been averaged over no revolutions of the satellite.
1017 * (mean motion)
1018 *
1019 * author : david vallado 719-573-2600 28 jun 2005
1020 *
1021 * inputs :
1022 * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 -
1023 * dedt -
1024 * del1, del2, del3 -
1025 * didt -
1026 * dmdt -
1027 * dnodt -
1028 * domdt -
1029 * irez - flag for resonance 0-none, 1-one day, 2-half day
1030 * argpo - argument of perigee
1031 * argpdot - argument of perigee dot (rate)
1032 * t - time
1033 * tc -
1034 * gsto - gst
1035 * xfact -
1036 * xlamo -
1037 * no - mean motion
1038 * atime -
1039 * em - eccentricity
1040 * ft -
1041 * argpm - argument of perigee
1042 * inclm - inclination
1043 * xli -
1044 * mm - mean anomaly
1045 * xni - mean motion
1046 * nodem - right ascension of ascending node
1047 *
1048 * outputs :
1049 * atime -
1050 * em - eccentricity
1051 * argpm - argument of perigee
1052 * inclm - inclination
1053 * xli -
1054 * mm - mean anomaly
1055 * xni -
1056 * nodem - right ascension of ascending node
1057 * dndt -
1058 * nm - mean motion
1059 *
1060 * locals :
1061 * delt -
1062 * ft -
1063 * theta -
1064 * x2li -
1065 * x2omi -
1066 * xl -
1067 * xldot -
1068 * xnddt -
1069 * xndt -
1070 * xomi -
1071 *
1072 * coupling :
1073 * none -
1074 *
1075 * references :
1076 * hoots, roehrich, norad spacetrack report #3 1980
1077 * hoots, norad spacetrack report #6 1986
1078 * hoots, schumacher and glover 2004
1079 * vallado, crawford, hujsak, kelso 2006
1080 ----------------------------------------------------------------------------*/
1081 // nm - also added as an input since it may not be changed and it needs to retain its old value
1082 // returns array with these values:
1083 // [em, argpm, inclm, mm,nodem, dndt, nm]
1084 private static double[] dspace(
1085 int irez,
1086 double d2201, double d2211, double d3210, double d3222, double d4410,
1087 double d4422, double d5220, double d5232, double d5421, double d5433,
1088 double dedt, double del1, double del2, double del3, double didt,
1089 double dmdt, double dnodt, double domdt, double argpo, double argpdot,
1090 double t, double tc, double gsto, double xfact, double xlamo,
1091 double no,//
1092 // These variables are both inputs and returned as outputs
1093 double em,
1094 double argpm,
1095 double inclm,
1096 double mm,
1097 double nodem,
1098 // data that is retained with the SGP4SatData object
1099 SGP4SatData satrec,
1100 // double& atime, double& em, double& argpm, double& inclm, double& xli,
1101 // double& mm, double& xni, double& nodem, double& dndt, double& nm
1102 double nm // input and output
1103 )
1104 {
1105 // variables that are both inputs and outputs! included in SGP4SatData
1106 // atime, em, argpm, inclm, xli, mm, xni, nodem
1107
1108 // Return variables (that are not also inputs)---------------------
1109 double dndt;
1110
1111 final double twopi = 2.0 * pi;
1112 int iretn, iret;
1113 double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi, g22, g32,
1114 g44, g52, g54, fasx2, fasx4, fasx6, rptim, step2, stepn, stepp;
1115
1116 // SEG -- it is okay to initalize these values here as they are assigned values
1117 xndt = 0;
1118 xnddt = 0;
1119 xldot = 0;
1120
1121 fasx2 = 0.13130908;
1122 fasx4 = 2.8843198;
1123 fasx6 = 0.37448087;
1124 g22 = 5.7686396;
1125 g32 = 0.95240898;
1126 g44 = 1.8014998;
1127 g52 = 1.0508330;
1128 g54 = 4.4108898;
1129 rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec
1130 stepp = 720.0;
1131 stepn = -720.0;
1132 step2 = 259200.0;
1133
1134 /* ----------- calculate deep space resonance effects ----------- */
1135 dndt = 0.0;
1136 theta = ((gsto + tc * rptim) % twopi);
1137 em = em + dedt * t;
1138
1139 inclm = inclm + didt * t;
1140 argpm = argpm + domdt * t;
1141 nodem = nodem + dnodt * t;
1142 mm = mm + dmdt * t;
1143
1144 // sgp4fix for negative inclinations
1145 // the following if statement should be commented out
1146 // if (inclm < 0.0)
1147 // {
1148 // inclm = -inclm;
1149 // argpm = argpm - pi;
1150 // nodem = nodem + pi;
1151 // }
1152
1153 /* - update resonances : numerical (euler-maclaurin) integration - */
1154 /* ------------------------- epoch restart ---------------------- */
1155 // sgp4fix for propagator problems
1156 // the following integration works for negative time steps and periods
1157 // the specific changes are unknown because the original code was so convoluted
1158
1159 // sgp4fix take out atime = 0.0 and fix for faster operation
1160 ft = 0.0;
1161 if(irez != 0)
1162 {
1163 // sgp4fix streamline check
1164 if((satrec.atime == 0.0) || (t * satrec.atime <= 0.0) || (Math.abs(t) < Math.abs(satrec.atime)))
1165 {
1166 satrec.atime = 0.0;
1167 satrec.xni = no;
1168 satrec.xli = xlamo;
1169 }
1170 // sgp4fix move check outside loop
1171 if(t > 0.0)
1172 {
1173 delt = stepp;
1174 }
1175 else
1176 {
1177 delt = stepn;
1178 }
1179
1180 iretn = 381; // added for do loop
1181 iret = 0; // added for loop
1182 while(iretn == 381)
1183 {
1184 /* ------------------- dot terms calculated ------------- */
1185 /* ----------- near - synchronous resonance terms ------- */
1186 if(irez != 2)
1187 {
1188 xndt = del1 * Math.sin(satrec.xli - fasx2) + del2 * Math.sin(2.0 * (satrec.xli - fasx4)) +
1189 del3 * Math.sin(3.0 * (satrec.xli - fasx6));
1190 xldot = satrec.xni + xfact;
1191 xnddt = del1 * Math.cos(satrec.xli - fasx2) +
1192 2.0 * del2 * Math.cos(2.0 * (satrec.xli - fasx4)) +
1193 3.0 * del3 * Math.cos(3.0 * (satrec.xli - fasx6));
1194 xnddt = xnddt * xldot;
1195 }
1196 else
1197 {
1198 /* --------- near - half-day resonance terms -------- */
1199 xomi = argpo + argpdot * satrec.atime;
1200 x2omi = xomi + xomi;
1201 x2li = satrec.xli + satrec.xli;
1202 xndt = d2201 * Math.sin(x2omi + satrec.xli - g22) + d2211 * Math.sin(satrec.xli - g22) +
1203 d3210 * Math.sin(xomi + satrec.xli - g32) + d3222 * Math.sin(-xomi + satrec.xli - g32) +
1204 d4410 * Math.sin(x2omi + x2li - g44) + d4422 * Math.sin(x2li - g44) +
1205 d5220 * Math.sin(xomi + satrec.xli - g52) + d5232 * Math.sin(-xomi + satrec.xli - g52) +
1206 d5421 * Math.sin(xomi + x2li - g54) + d5433 * Math.sin(-xomi + x2li - g54);
1207 xldot = satrec.xni + xfact;
1208 xnddt = d2201 * Math.cos(x2omi + satrec.xli - g22) + d2211 * Math.cos(satrec.xli - g22) +
1209 d3210 * Math.cos(xomi + satrec.xli - g32) + d3222 * Math.cos(-xomi + satrec.xli - g32) +
1210 d5220 * Math.cos(xomi + satrec.xli - g52) + d5232 * Math.cos(-xomi + satrec.xli - g52) +
1211 2.0 * (d4410 * Math.cos(x2omi + x2li - g44) +
1212 d4422 * Math.cos(x2li - g44) + d5421 * Math.cos(xomi + x2li - g54) +
1213 d5433 * Math.cos(-xomi + x2li - g54));
1214 xnddt = xnddt * xldot;
1215 }
1216
1217 /* ----------------------- integrator ------------------- */
1218 // sgp4fix move end checks to end of routine
1219 if(Math.abs(t - satrec.atime) >= stepp)
1220 {
1221 iret = 0;
1222 iretn = 381;
1223 }
1224 else // exit here
1225 {
1226 ft = t - satrec.atime;
1227 iretn = 0;
1228 }
1229
1230 if(iretn == 381)
1231 {
1232 satrec.xli = satrec.xli + xldot * delt + xndt * step2;
1233 satrec.xni = satrec.xni + xndt * delt + xnddt * step2;
1234 satrec.atime = satrec.atime + delt;
1235 }
1236 } // while iretn = 381
1237
1238 nm = satrec.xni + xndt * ft + xnddt * ft * ft * 0.5;
1239 xl = satrec.xli + xldot * ft + xndt * ft * ft * 0.5;
1240 if(irez != 1)
1241 {
1242 mm = xl - 2.0 * nodem + 2.0 * theta;
1243 dndt = nm - no;
1244 }
1245 else
1246 {
1247 mm = xl - nodem - argpm + theta;
1248 dndt = nm - no;
1249 }
1250 nm = no + dndt;
1251 }
1252
1253 return new double[]
1254 {
1255 em, argpm, inclm, mm, nodem, dndt, nm
1256 };
1257
1258 //#include "debug4.cpp"
1259 } // end dsspace
1260
1261 /** -----------------------------------------------------------------------------
1262 *
1263 * procedure initl
1264 *
1265 * this procedure initializes the spg4 propagator. all the initialization is
1266 * consolidated here instead of having multiple loops inside other routines.
1267 *
1268 * author : david vallado 719-573-2600 28 jun 2005
1269 *
1270 * inputs :
1271 * ecco - eccentricity 0.0 - 1.0
1272 * epoch - epoch time in days from jan 0, 1950. 0 hr
1273 * inclo - inclination of satellite
1274 * no - mean motion of satellite
1275 * satn - satellite number
1276 *
1277 * outputs :
1278 * ainv - 1.0 / a
1279 * ao - semi major axis
1280 * con41 -
1281 * con42 - 1.0 - 5.0 cos(i)
1282 * cosio - cosine of inclination
1283 * cosio2 - cosio squared
1284 * eccsq - eccentricity squared
1285 * method - flag for deep space 'd', 'n'
1286 * omeosq - 1.0 - ecco * ecco
1287 * posq - semi-parameter squared
1288 * rp - radius of perigee
1289 * rteosq - square root of (1.0 - ecco*ecco)
1290 * sinio - sine of inclination
1291 * gsto - gst at time of observation rad
1292 * no - mean motion of satellite
1293 *
1294 * locals :
1295 * ak -
1296 * d1 -
1297 * del -
1298 * adel -
1299 * po -
1300 *
1301 * coupling :
1302 * getgravconst
1303 * gstime - find greenwich sidereal time from the julian date
1304 *
1305 * references :
1306 * hoots, roehrich, norad spacetrack report #3 1980
1307 * hoots, norad spacetrack report #6 1986
1308 * hoots, schumacher and glover 2004
1309 * vallado, crawford, hujsak, kelso 2006
1310 * ----------------------------------------------------------------------------
1311 * @param satn satellite number
1312 * @param whichconst which constants set to use
1313 * @param ecco eccentricity (0,1)
1314 * @param epoch epoch time in days from jan 0, 1950. 0 hr
1315 * @param inclo inclination of satellite
1316 * @param satrec satellite object that stores needed SGP4 data
1317 * @return double array with these values: [ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio]
1318 */
1319 // outputs not stored in SGP4SatData and are returned by this function:
1320 // [ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio]
1321 public static double[] initl(
1322 int satn, Gravconsttype whichconst,
1323 double ecco, double epoch, double inclo,
1324 //
1325 SGP4SatData satrec // //double& no,
1326 // //char& method,
1327 // double& ainv, double& ao,
1328 // //double& con41,
1329 // double& con42, double& cosio,
1330 // double& cosio2,double& eccsq, double& omeosq, double& posq,
1331 // double& rp, double& rteosq,double& sinio ,
1332 // //double& gsto, char opsmode
1333 )
1334 {
1335 // no is an input and an output
1336 // return variables ----------------
1337 double ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio;
1338
1339 /* --------------------- local variables ------------------------ */
1340 double ak, d1, del, adel, po, x2o3, j2, xke,
1341 tumin, mu, radiusearthkm, j3, j4, j3oj2;
1342
1343 // sgp4fix use old way of finding gst
1344 double ds70;
1345 double ts70, tfrac, c1, thgr70, fk5r, c1p2p, thgr, thgro;
1346 final double twopi = 2.0 * pi;
1347
1348 /* ----------------------- earth constants ---------------------- */
1349 // sgp4fix identify constants and allow alternate values
1350 double[] temp2 = getgravconst(whichconst);
1351 tumin = temp2[0];
1352 mu = temp2[1];
1353 radiusearthkm = temp2[2];
1354 xke = temp2[3];
1355 j2 = temp2[4];
1356 j3 = temp2[5];
1357 j4 = temp2[6];
1358 j3oj2 = temp2[7];
1359
1360 x2o3 = 2.0 / 3.0;
1361
1362 /* ------------- calculate auxillary epoch quantities ---------- */
1363 eccsq = ecco * ecco;
1364 omeosq = 1.0 - eccsq;
1365 rteosq = Math.sqrt(omeosq);
1366 cosio = Math.cos(inclo);
1367 cosio2 = cosio * cosio;
1368
1369 /* ------------------ un-kozai the mean motion ----------------- */
1370 ak = Math.pow(xke / satrec.no, x2o3);
1371 d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq);
1372 del = d1 / (ak * ak);
1373 adel = ak * (1.0 - del * del - del *
1374 (1.0 / 3.0 + 134.0 * del * del / 81.0));
1375 del = d1 / (adel * adel);
1376 satrec.no = satrec.no / (1.0 + del);
1377
1378 ao = Math.pow(xke / satrec.no, x2o3);
1379 sinio = Math.sin(inclo);
1380 po = ao * omeosq;
1381 con42 = 1.0 - 5.0 * cosio2;
1382 satrec.con41 = -con42 - cosio2 - cosio2;
1383 ainv = 1.0 / ao;
1384 posq = po * po;
1385 rp = ao * (1.0 - ecco);
1386 satrec.method = 'n';
1387
1388 // sgp4fix modern approach to finding sidereal time
1389 if(satrec.operationmode == 'a')
1390 {
1391 // sgp4fix use old way of finding gst
1392 // count integer number of days from 0 jan 1970
1393 ts70 = epoch - 7305.0;
1394 ds70 = Math.floor(ts70 + 1.0e-8);
1395 tfrac = ts70 - ds70;
1396 // find greenwich location at epoch
1397 c1 = 1.72027916940703639e-2;
1398 thgr70 = 1.7321343856509374;
1399 fk5r = 5.07551419432269442e-15;
1400 c1p2p = c1 + twopi;
1401 satrec.gsto = ((thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r) % twopi);
1402 if(satrec.gsto < 0.0)
1403 {
1404 satrec.gsto = satrec.gsto + twopi;
1405 }
1406 }
1407 else
1408 {
1409 satrec.gsto = gstime(epoch + 2433281.5);
1410 }
1411
1412 return new double[]
1413 {
1414 ainv, ao, con42, cosio, cosio2, eccsq, omeosq, posq, rp, rteosq, sinio
1415 };
1416
1417 //#include "debug5.cpp"
1418 } // end initl
1419
1420 /**-----------------------------------------------------------------------------
1421 * This method is called from within SGP4utils.readTLEandIniSGP4 and therefore not
1422 * typically called elsewhere.
1423 * -----------------------------------------------------------------------------
1424 *
1425 * procedure sgp4init
1426 *
1427 * this procedure initializes variables for sgp4.
1428 *
1429 * author : david vallado 719-573-2600 28 jun 2005
1430 *
1431 * inputs :
1432 * opsmode - mode of operation afspc or improved 'a', 'i'
1433 * whichconst - which set of constants to use 72, 84
1434 * satn - satellite number
1435 * bstar - sgp4 type drag coefficient kg/m2er
1436 * ecco - eccentricity
1437 * epoch - epoch time in days from jan 0, 1950. 0 hr
1438 * argpo - argument of perigee (output if ds)
1439 * inclo - inclination
1440 * mo - mean anomaly (output if ds)
1441 * no - mean motion
1442 * nodeo - right ascension of ascending node
1443 *
1444 * outputs :
1445 * satrec - common values for subsequent calls
1446 * return code - non-zero on error.
1447 * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er
1448 * 2 - mean motion less than 0.0
1449 * 3 - pert elements, ecc < 0.0 or ecc > 1.0
1450 * 4 - semi-latus rectum < 0.0
1451 * 5 - epoch elements are sub-orbital
1452 * 6 - satellite has decayed
1453 *
1454 * locals :
1455 * cnodm , snodm , cosim , sinim , cosomm , sinomm
1456 * cc1sq , cc2 , cc3
1457 * coef , coef1
1458 * cosio4 -
1459 * day -
1460 * dndt -
1461 * em - eccentricity
1462 * emsq - eccentricity squared
1463 * eeta -
1464 * etasq -
1465 * gam -
1466 * argpm - argument of perigee
1467 * nodem -
1468 * inclm - inclination
1469 * mm - mean anomaly
1470 * nm - mean motion
1471 * perige - perigee
1472 * pinvsq -
1473 * psisq -
1474 * qzms24 -
1475 * rtemsq -
1476 * s1, s2, s3, s4, s5, s6, s7 -
1477 * sfour -
1478 * ss1, ss2, ss3, ss4, ss5, ss6, ss7 -
1479 * sz1, sz2, sz3
1480 * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 -
1481 * tc -
1482 * temp -
1483 * temp1, temp2, temp3 -
1484 * tsi -
1485 * xpidot -
1486 * xhdot1 -
1487 * z1, z2, z3 -
1488 * z11, z12, z13, z21, z22, z23, z31, z32, z33 -
1489 *
1490 * coupling :
1491 * getgravconst-
1492 * initl -
1493 * dscom -
1494 * dpper -
1495 * dsinit -
1496 * sgp4 -
1497 *
1498 * references :
1499 * hoots, roehrich, norad spacetrack report #3 1980
1500 * hoots, norad spacetrack report #6 1986
1501 * hoots, schumacher and glover 2004
1502 * vallado, crawford, hujsak, kelso 2006
1503 * ----------------------------------------------------------------------------
1504 * @param whichconst
1505 * @param opsmode
1506 * @param satn
1507 * @param epoch
1508 * @param xbstar
1509 * @param xecco
1510 * @param xargpo
1511 * @param xinclo
1512 * @param xmo
1513 * @param xno
1514 * @param xnodeo
1515 * @param satrec
1516 * @return if initialization was successful
1517 */
1518 public static boolean sgp4init(
1519 Gravconsttype whichconst, char opsmode, final int satn, final double epoch,
1520 final double xbstar, final double xecco, final double xargpo,
1521 final double xinclo, final double xmo, final double xno,
1522 final double xnodeo, SGP4SatData satrec)
1523 {
1524 /* --------------------- local variables ------------------------ */
1525 double ao, ainv, con42, cosio, sinio, cosio2, eccsq,
1526 omeosq, posq, rp, rteosq,
1527 cnodm, snodm, cosim, sinim, cosomm, sinomm, cc1sq,
1528 cc2, cc3, coef, coef1, cosio4, day, dndt,
1529 em, emsq, eeta, etasq, gam, argpm, nodem,
1530 inclm, mm, nm, perige, pinvsq, psisq, qzms24,
1531 rtemsq, s1, s2, s3, s4, s5, s6,
1532 s7, sfour, ss1, ss2, ss3, ss4, ss5,
1533 ss6, ss7, sz1, sz2, sz3, sz11, sz12,
1534 sz13, sz21, sz22, sz23, sz31, sz32, sz33,
1535 tc, temp, temp1, temp2, temp3, tsi, xpidot,
1536 xhdot1, z1, z2, z3, z11, z12, z13,
1537 z21, z22, z23, z31, z32, z33,
1538 qzms2t, ss, j2, j3oj2, j4, x2o3, //r[3], v[3],
1539 tumin, mu, radiusearthkm, xke, j3;
1540 double[] r = new double[3];
1541 double[] v = new double[3];
1542
1543 /* ------------------------ initialization --------------------- */
1544 // sgp4fix divisor for divide by zero check on inclination
1545 // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
1546 // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency
1547 final double temp4 = 1.5e-12;
1548
1549 /* ----------- set all near earth variables to zero ------------ */
1550 satrec.isimp = 0;
1551 satrec.method = 'n';
1552 satrec.aycof = 0.0;
1553 satrec.con41 = 0.0;
1554 satrec.cc1 = 0.0;
1555 satrec.cc4 = 0.0;
1556 satrec.cc5 = 0.0;
1557 satrec.d2 = 0.0;
1558 satrec.d3 = 0.0;
1559 satrec.d4 = 0.0;
1560 satrec.delmo = 0.0;
1561 satrec.eta = 0.0;
1562 satrec.argpdot = 0.0;
1563 satrec.omgcof = 0.0;
1564 satrec.sinmao = 0.0;
1565 satrec.t = 0.0;
1566 satrec.t2cof = 0.0;
1567 satrec.t3cof = 0.0;
1568 satrec.t4cof = 0.0;
1569 satrec.t5cof = 0.0;
1570 satrec.x1mth2 = 0.0;
1571 satrec.x7thm1 = 0.0;
1572 satrec.mdot = 0.0;
1573 satrec.nodedot = 0.0;
1574 satrec.xlcof = 0.0;
1575 satrec.xmcof = 0.0;
1576 satrec.nodecf = 0.0;
1577
1578 /* ----------- set all deep space variables to zero ------------ */
1579 satrec.irez = 0;
1580 satrec.d2201 = 0.0;
1581 satrec.d2211 = 0.0;
1582 satrec.d3210 = 0.0;
1583 satrec.d3222 = 0.0;
1584 satrec.d4410 = 0.0;
1585 satrec.d4422 = 0.0;
1586 satrec.d5220 = 0.0;
1587 satrec.d5232 = 0.0;
1588 satrec.d5421 = 0.0;
1589 satrec.d5433 = 0.0;
1590 satrec.dedt = 0.0;
1591 satrec.del1 = 0.0;
1592 satrec.del2 = 0.0;
1593 satrec.del3 = 0.0;
1594 satrec.didt = 0.0;
1595 satrec.dmdt = 0.0;
1596 satrec.dnodt = 0.0;
1597 satrec.domdt = 0.0;
1598 satrec.e3 = 0.0;
1599 satrec.ee2 = 0.0;
1600 satrec.peo = 0.0;
1601 satrec.pgho = 0.0;
1602 satrec.pho = 0.0;
1603 satrec.pinco = 0.0;
1604 satrec.plo = 0.0;
1605 satrec.se2 = 0.0;
1606 satrec.se3 = 0.0;
1607 satrec.sgh2 = 0.0;
1608 satrec.sgh3 = 0.0;
1609 satrec.sgh4 = 0.0;
1610 satrec.sh2 = 0.0;
1611 satrec.sh3 = 0.0;
1612 satrec.si2 = 0.0;
1613 satrec.si3 = 0.0;
1614 satrec.sl2 = 0.0;
1615 satrec.sl3 = 0.0;
1616 satrec.sl4 = 0.0;
1617 satrec.gsto = 0.0;
1618 satrec.xfact = 0.0;
1619 satrec.xgh2 = 0.0;
1620 satrec.xgh3 = 0.0;
1621 satrec.xgh4 = 0.0;
1622 satrec.xh2 = 0.0;
1623 satrec.xh3 = 0.0;
1624 satrec.xi2 = 0.0;
1625 satrec.xi3 = 0.0;
1626 satrec.xl2 = 0.0;
1627 satrec.xl3 = 0.0;
1628 satrec.xl4 = 0.0;
1629 satrec.xlamo = 0.0;
1630 satrec.zmol = 0.0;
1631 satrec.zmos = 0.0;
1632 satrec.atime = 0.0;
1633 satrec.xli = 0.0;
1634 satrec.xni = 0.0;
1635
1636 // sgp4fix - note the following variables are also passed directly via satrec.
1637 // it is possible to streamline the sgp4init call by deleting the "x"
1638 // variables, but the user would need to set the satrec.* values first. we
1639 // include the additional assignments in case twoline2rv is not used.
1640 satrec.bstar = xbstar;
1641 satrec.ecco = xecco;
1642 satrec.argpo = xargpo;
1643 satrec.inclo = xinclo;
1644 satrec.mo = xmo;
1645 satrec.no = xno;
1646 satrec.nodeo = xnodeo;
1647
1648 // sgp4fix add opsmode
1649 satrec.operationmode = opsmode;
1650
1651 // SEG -- also save gravity constant type
1652 satrec.gravconsttype = whichconst;
1653
1654 /* ------------------------ earth constants ----------------------- */
1655 // sgp4fix identify constants and allow alternate values
1656 double[] temp5 = getgravconst(whichconst);//, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 );
1657 tumin = temp5[0];
1658 mu = temp5[1];
1659 radiusearthkm = temp5[2];
1660 xke = temp5[3];
1661 j2 = temp5[4];
1662 j3 = temp5[5];
1663 j4 = temp5[6];
1664 j3oj2 = temp5[7];
1665
1666 ss = 78.0 / radiusearthkm + 1.0;
1667 qzms2t = Math.pow(((120.0 - 78.0) / radiusearthkm), 4);
1668 x2o3 = 2.0 / 3.0;
1669
1670 satrec.init = 'y';
1671 satrec.t = 0.0;
1672
1673 double[] ttemp = initl(satn, whichconst, satrec.ecco, epoch, satrec.inclo, satrec);
1674 ainv = ttemp[0];
1675 ao = ttemp[1];
1676 con42 = ttemp[2];
1677 cosio = ttemp[3];
1678 cosio2 = ttemp[4];
1679 eccsq = ttemp[5];
1680 omeosq = ttemp[6];
1681 posq = ttemp[7];
1682 rp = ttemp[8];
1683 rteosq = ttemp[9];
1684 sinio = ttemp[10];
1685
1686 satrec.error = 0;
1687
1688 // sgp4fix remove this check as it is unnecessary
1689 // the mrt check in sgp4 handles decaying satellite cases even if the starting
1690 // condition is below the surface of te earth
1691 // if (rp < 1.0)
1692 // {
1693 // printf("# *** satn%d epoch elts sub-orbital ***\n", satn);
1694 // satrec.error = 5;
1695 // }
1696
1697 if((omeosq >= 0.0) || (satrec.no >= 0.0))
1698 {
1699 satrec.isimp = 0;
1700 if(rp < (220.0 / radiusearthkm + 1.0))
1701 {
1702 satrec.isimp = 1;
1703 }
1704 sfour = ss;
1705 qzms24 = qzms2t;
1706 perige = (rp - 1.0) * radiusearthkm;
1707
1708 /* - for perigees below 156 km, s and qoms2t are altered - */
1709 if(perige < 156.0)
1710 {
1711 sfour = perige - 78.0;
1712 if(perige < 98.0)
1713 {
1714 sfour = 20.0;
1715 }
1716 qzms24 = Math.pow(((120.0 - sfour) / radiusearthkm), 4.0);
1717 sfour = sfour / radiusearthkm + 1.0;
1718 }
1719 pinvsq = 1.0 / posq;
1720
1721 tsi = 1.0 / (ao - sfour);
1722 satrec.eta = ao * satrec.ecco * tsi;
1723 etasq = satrec.eta * satrec.eta;
1724 eeta = satrec.ecco * satrec.eta;
1725 psisq = Math.abs(1.0 - etasq);
1726 coef = qzms24 * Math.pow(tsi, 4.0);
1727 coef1 = coef / Math.pow(psisq, 3.5);
1728 cc2 = coef1 * satrec.no * (ao * (1.0 + 1.5 * etasq + eeta *
1729 (4.0 + etasq)) + 0.375 * j2 * tsi / psisq * satrec.con41 *
1730 (8.0 + 3.0 * etasq * (8.0 + etasq)));
1731 satrec.cc1 = satrec.bstar * cc2;
1732 cc3 = 0.0;
1733 if(satrec.ecco > 1.0e-4)
1734 {
1735 cc3 = -2.0 * coef * tsi * j3oj2 * satrec.no * sinio / satrec.ecco;
1736 }
1737 satrec.x1mth2 = 1.0 - cosio2;
1738 satrec.cc4 = 2.0 * satrec.no * coef1 * ao * omeosq *
1739 (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecco *
1740 (0.5 + 2.0 * etasq) - j2 * tsi / (ao * psisq) *
1741 (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq *
1742 (1.5 - 0.5 * eeta)) + 0.75 * satrec.x1mth2 *
1743 (2.0 * etasq - eeta * (1.0 + etasq)) * Math.cos(2.0 * satrec.argpo)));
1744 satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 *
1745 (etasq + eeta) + eeta * etasq);
1746 cosio4 = cosio2 * cosio2;
1747 temp1 = 1.5 * j2 * pinvsq * satrec.no;
1748 temp2 = 0.5 * temp1 * j2 * pinvsq;
1749 temp3 = -0.46875 * j4 * pinvsq * pinvsq * satrec.no;
1750 satrec.mdot = satrec.no + 0.5 * temp1 * rteosq * satrec.con41 + 0.0625 *
1751 temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4);
1752 satrec.argpdot = -0.5 * temp1 * con42 + 0.0625 * temp2 *
1753 (7.0 - 114.0 * cosio2 + 395.0 * cosio4) +
1754 temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4);
1755 xhdot1 = -temp1 * cosio;
1756 satrec.nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) +
1757 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio;
1758 xpidot = satrec.argpdot + satrec.nodedot;
1759 satrec.omgcof = satrec.bstar * cc3 * Math.cos(satrec.argpo);
1760 satrec.xmcof = 0.0;
1761 if(satrec.ecco > 1.0e-4)
1762 {
1763 satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta;
1764 }
1765 satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1;
1766 satrec.t2cof = 1.5 * satrec.cc1;
1767 // sgp4fix for divide by zero with xinco = 180 deg
1768 if(Math.abs(cosio + 1.0) > 1.5e-12)
1769 {
1770 satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio);
1771 }
1772 else
1773 {
1774 satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4;
1775 }
1776 satrec.aycof = -0.5 * j3oj2 * sinio;
1777 satrec.delmo = Math.pow((1.0 + satrec.eta * Math.cos(satrec.mo)), 3);
1778 satrec.sinmao = Math.sin(satrec.mo);
1779 satrec.x7thm1 = 7.0 * cosio2 - 1.0;
1780
1781 /* --------------- deep space initialization ------------- */
1782 if((2 * pi / satrec.no) >= 225.0)
1783 {
1784 satrec.method = 'd';
1785 satrec.isimp = 1;
1786 tc = 0.0;
1787 inclm = satrec.inclo;
1788
1789 double[] ttemp2 = dscom(epoch, satrec.ecco, satrec.argpo, tc, satrec.inclo, satrec.nodeo, satrec.no, satrec);
1790 // save ouput vars
1791 snodm = ttemp2[0];
1792 cnodm = ttemp2[1];
1793 sinim = ttemp2[2];
1794 cosim = ttemp2[3];
1795 sinomm = ttemp2[4];
1796 cosomm = ttemp2[5];
1797 day = ttemp2[6];
1798 em = ttemp2[7];
1799 emsq = ttemp2[8];
1800 gam = ttemp2[9];
1801 rtemsq = ttemp2[10];
1802 s1 = ttemp2[11];
1803 s2 = ttemp2[12];
1804 s3 = ttemp2[13];
1805 s4 = ttemp2[14];
1806 s5 = ttemp2[15];
1807 s6 = ttemp2[16];
1808 s7 = ttemp2[17];
1809 ss1 = ttemp2[18];
1810 ss2 = ttemp2[19];
1811 ss3 = ttemp2[20];
1812 ss4 = ttemp2[21];
1813 ss5 = ttemp2[22];
1814 ss6 = ttemp2[23];
1815 ss7 = ttemp2[24];
1816 sz1 = ttemp2[25];
1817 sz2 = ttemp2[26];
1818 sz3 = ttemp2[27];
1819 sz11 = ttemp2[28];
1820 sz12 = ttemp2[29];
1821 sz13 = ttemp2[30];
1822 sz21 = ttemp2[31];
1823 sz22 = ttemp2[32];
1824 sz23 = ttemp2[33];
1825 sz31 = ttemp2[34];
1826 sz32 = ttemp2[35];
1827 sz33 = ttemp2[36];
1828 nm = ttemp2[37];
1829 z1 = ttemp2[38];
1830 z2 = ttemp2[39];
1831 z3 = ttemp2[40];
1832 z11 = ttemp2[41];
1833 z12 = ttemp2[42];
1834 z13 = ttemp2[43];
1835 z21 = ttemp2[44];
1836 z22 = ttemp2[45];
1837 z23 = ttemp2[46];
1838 z31 = ttemp2[47];
1839 z32 = ttemp2[48];
1840 z33 = ttemp2[49];
1841
1842 //dpper(satrec);
1843 ttemp2 = dpper(
1844 satrec.e3, satrec.ee2, satrec.peo, satrec.pgho,
1845 satrec.pho, satrec.pinco, satrec.plo, satrec.se2,
1846 satrec.se3, satrec.sgh2, satrec.sgh3, satrec.sgh4,
1847 satrec.sh2, satrec.sh3, satrec.si2, satrec.si3,
1848 satrec.sl2, satrec.sl3, satrec.sl4, satrec.t,
1849 satrec.xgh2, satrec.xgh3, satrec.xgh4, satrec.xh2,
1850 satrec.xh3, satrec.xi2, satrec.xi3, satrec.xl2,
1851 satrec.xl3, satrec.xl4, satrec.zmol, satrec.zmos, inclm, satrec.init,
1852 satrec.ecco, satrec.inclo, satrec.nodeo, satrec.argpo, satrec.mo,
1853 satrec.operationmode);
1854 satrec.ecco = ttemp2[0];
1855 satrec.inclo = ttemp2[1];
1856 satrec.nodeo = ttemp2[2];
1857 satrec.argpo = ttemp2[3];
1858 satrec.mo = ttemp2[4];
1859
1860 argpm = 0.0;
1861 nodem = 0.0;
1862 mm = 0.0;
1863
1864 double[] ttemp3 = dsinit(whichconst,
1865 cosim, emsq, satrec.argpo, s1, s2, s3, s4, s5, sinim, ss1, ss2, ss3, ss4,
1866 ss5, sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33, satrec.t, tc,
1867 satrec.gsto, satrec.mo, satrec.mdot, satrec.no, satrec.nodeo,
1868 satrec.nodedot, xpidot, z1, z3, z11, z13, z21, z23, z31, z33,
1869 satrec.ecco, eccsq,
1870 satrec,
1871 em, argpm, inclm, mm, nm, nodem);
1872 em = ttemp3[0];
1873 argpm = ttemp3[1];
1874 inclm = ttemp3[2];
1875 mm = ttemp3[3];
1876 nm = ttemp3[4];
1877 nodem = ttemp3[5];
1878 dndt = ttemp3[6];
1879 }
1880
1881 /* ----------- set variables if not deep space ----------- */
1882 if(satrec.isimp != 1)
1883 {
1884 cc1sq = satrec.cc1 * satrec.cc1;
1885 satrec.d2 = 4.0 * ao * tsi * cc1sq;
1886 temp = satrec.d2 * tsi * satrec.cc1 / 3.0;
1887 satrec.d3 = (17.0 * ao + sfour) * temp;
1888 satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) *
1889 satrec.cc1;
1890 satrec.t3cof = satrec.d2 + 2.0 * cc1sq;
1891 satrec.t4cof = 0.25 * (3.0 * satrec.d3 + satrec.cc1 *
1892 (12.0 * satrec.d2 + 10.0 * cc1sq));
1893 satrec.t5cof = 0.2 * (3.0 * satrec.d4 +
1894 12.0 * satrec.cc1 * satrec.d3 +
1895 6.0 * satrec.d2 * satrec.d2 +
1896 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq));
1897 }
1898 } // if omeosq = 0 ...
1899
1900 /* finally propogate to zero epoch to initialize all others. */
1901 // sgp4fix take out check to let satellites process until they are actually below earth surface
1902 // if(satrec.error == 0)
1903 boolean sgp4Error = sgp4(satrec, 0.0, r, v); // SEG removed gravity constant passing
1904
1905 satrec.init = 'n';
1906
1907 //#include "debug6.cpp"
1908 //sgp4fix return boolean. satrec.error contains any error codes
1909 return sgp4Error;
1910 } // end sgp4init
1911
1912
1913
1914 /**
1915 * Similar to sgp4(..) but time parameter is the Julian Date to the propagated to.
1916 * This method was not orgiinally in the CSSI C++ version.
1917 * @param satrec satellite SGP4 data object
1918 * @param jd Julian Date
1919 * @param r position vector [km] return array (needs to be of size 3)
1920 * @param v velocity [km/sec] return array (needs to be of size 3)
1921 * @return
1922 */
1923 public static boolean sgp4Prop2JD(SGP4SatData satrec, double jd, double[] r, double[] v)
1924 {
1925 double tminSinceEpoch = (jd - satrec.jdsatepoch)*24.0*60.0;
1926 return sgp4(satrec, tminSinceEpoch,r, v);
1927 }
1928
1929 /*-----------------------------------------------------------------------------
1930 *
1931 * procedure sgp4
1932 *
1933 * this procedure is the sgp4 prediction model from space command. this is an
1934 * updated and combined version of sgp4 and sdp4, which were originally
1935 * published separately in spacetrack report #3. this version follows the
1936 * methodology from the aiaa paper (2006) describing the history and
1937 * development of the code.
1938 *
1939 * author : david vallado 719-573-2600 28 jun 2005
1940 *
1941 * inputs :
1942 * satrec - initialised structure from sgp4init() call.
1943 * tsince - time eince epoch (minutes)
1944 *
1945 * outputs :
1946 * r - position vector km
1947 * v - velocity km/sec
1948 * return code - non-zero on error.
1949 * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er
1950 * 2 - mean motion less than 0.0
1951 * 3 - pert elements, ecc < 0.0 or ecc > 1.0
1952 * 4 - semi-latus rectum < 0.0
1953 * 5 - epoch elements are sub-orbital
1954 * 6 - satellite has decayed
1955 *
1956 * locals :
1957 * am -
1958 * axnl, aynl -
1959 * betal -
1960 * cosim , sinim , cosomm , sinomm , cnod , snod , cos2u ,
1961 * sin2u , coseo1 , sineo1 , cosi , sini , cosip , sinip ,
1962 * cosisq , cossu , sinsu , cosu , sinu
1963 * delm -
1964 * delomg -
1965 * dndt -
1966 * eccm -
1967 * emsq -
1968 * ecose -
1969 * el2 -
1970 * eo1 -
1971 * eccp -
1972 * esine -
1973 * argpm -
1974 * argpp -
1975 * omgadf -
1976 * pl -
1977 * r -
1978 * rtemsq -
1979 * rdotl -
1980 * rl -
1981 * rvdot -
1982 * rvdotl -
1983 * su -
1984 * t2 , t3 , t4 , tc
1985 * tem5, temp , temp1 , temp2 , tempa , tempe , templ
1986 * u , ux , uy , uz , vx , vy , vz
1987 * inclm - inclination
1988 * mm - mean anomaly
1989 * nm - mean motion
1990 * nodem - right asc of ascending node
1991 * xinc -
1992 * xincp -
1993 * xl -
1994 * xlm -
1995 * mp -
1996 * xmdf -
1997 * xmx -
1998 * xmy -
1999 * nodedf -
2000 * xnode -
2001 * nodep -
2002 * np -
2003 *
2004 * coupling :
2005 * getgravconst-
2006 * dpper
2007 * dpspace
2008 *
2009 * references :
2010 * hoots, roehrich, norad spacetrack report #3 1980
2011 * hoots, norad spacetrack report #6 1986
2012 * hoots, schumacher and glover 2004
2013 * vallado, crawford, hujsak, kelso 2006
2014 * ----------------------------------------------------------------------------
2015 *
2016 * @param satrec satelite sgp4 data object
2017 * @param tsince time eince epoch (minutes)
2018 * @param r position vector [km] return array (needs to be of size 3)
2019 * @param v velocity [km/sec] return array (needs to be of size 3)
2020 * @return true if there were no errors, see satrec.error for error code
2021 */
2022 public static boolean sgp4(
2023 //Gravconsttype whichconst, // SEG removed as it should already be saved into satrec
2024 SGP4SatData satrec, double tsince,
2025 double[] r, double[] v)
2026 {
2027 double am, axnl, aynl, betal, cosim, cnod,
2028 cos2u, coseo1 = 0, cosi, cosip, cosisq, cossu, cosu,
2029 delm, delomg, em, emsq, ecose, el2, eo1,
2030 ep, esine, argpm, argpp, argpdf, pl, mrt = 0.0,
2031 mvt, rdotl, rl, rvdot, rvdotl, sinim,
2032 sin2u, sineo1 = 0, sini, sinip, sinsu, sinu,
2033 snod, su, t2, t3, t4, tem5, temp,
2034 temp1, temp2, tempa, tempe, templ, u, ux,
2035 uy, uz, vx, vy, vz, inclm, mm,
2036 nm, nodem, xinc, xincp, xl, xlm, mp,
2037 xmdf, xmx, xmy, nodedf, xnode, nodep, tc, dndt,
2038 twopi, x2o3, j2, j3, tumin, j4, xke, j3oj2, radiusearthkm,
2039 mu, vkmpersec;
2040 int ktr;
2041
2042 /* ------------------ set mathematical constants --------------- */
2043 // sgp4fix divisor for divide by zero check on inclination
2044 // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
2045 // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency
2046 final double temp4 = 1.5e-12;
2047 twopi = 2.0 * pi;
2048 x2o3 = 2.0 / 3.0;
2049 // sgp4fix identify constants and allow alternate values
2050 double[] temp5 = getgravconst(satrec.gravconsttype);//, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 );
2051 tumin = temp5[0];
2052 mu = temp5[1];
2053 radiusearthkm = temp5[2];
2054 xke = temp5[3];
2055 j2 = temp5[4];
2056 j3 = temp5[5];
2057 j4 = temp5[6];
2058 j3oj2 = temp5[7];
2059
2060 vkmpersec = radiusearthkm * xke / 60.0;
2061
2062 /* --------------------- clear sgp4 error flag ----------------- */
2063 satrec.t = tsince;
2064 satrec.error = 0;
2065
2066 /* ------- update for secular gravity and atmospheric drag ----- */
2067 xmdf = satrec.mo + satrec.mdot * satrec.t;
2068 argpdf = satrec.argpo + satrec.argpdot * satrec.t;
2069 nodedf = satrec.nodeo + satrec.nodedot * satrec.t;
2070 argpm = argpdf;
2071 mm = xmdf;
2072 t2 = satrec.t * satrec.t;
2073 nodem = nodedf + satrec.nodecf * t2;
2074 tempa = 1.0 - satrec.cc1 * satrec.t;
2075 tempe = satrec.bstar * satrec.cc4 * satrec.t;
2076 templ = satrec.t2cof * t2;
2077
2078 if(satrec.isimp != 1)
2079 {
2080 delomg = satrec.omgcof * satrec.t;
2081 delm = satrec.xmcof *
2082 (Math.pow((1.0 + satrec.eta * Math.cos(xmdf)), 3) -
2083 satrec.delmo);
2084 temp = delomg + delm;
2085 mm = xmdf + temp;
2086 argpm = argpdf - temp;
2087 t3 = t2 * satrec.t;
2088 t4 = t3 * satrec.t;
2089 tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 -
2090 satrec.d4 * t4;
2091 tempe = tempe + satrec.bstar * satrec.cc5 * (Math.sin(mm) -
2092 satrec.sinmao);
2093 templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof +
2094 satrec.t * satrec.t5cof);
2095 }
2096
2097 nm = satrec.no;
2098 em = satrec.ecco;
2099 inclm = satrec.inclo;
2100 if(satrec.method == 'd')
2101 {
2102 tc = satrec.t;
2103 double[] ttemp = dspace(
2104 satrec.irez,
2105 satrec.d2201, satrec.d2211, satrec.d3210,
2106 satrec.d3222, satrec.d4410, satrec.d4422,
2107 satrec.d5220, satrec.d5232, satrec.d5421,
2108 satrec.d5433, satrec.dedt, satrec.del1,
2109 satrec.del2, satrec.del3, satrec.didt,
2110 satrec.dmdt, satrec.dnodt, satrec.domdt,
2111 satrec.argpo, satrec.argpdot, satrec.t, tc,
2112 satrec.gsto, satrec.xfact, satrec.xlamo,
2113 satrec.no,
2114 em, argpm, inclm, mm, nodem,
2115 satrec,
2116 nm);
2117 // copy variables back
2118 em = ttemp[0];
2119 argpm = ttemp[1];
2120 inclm = ttemp[2];
2121 mm = ttemp[3];
2122 nodem = ttemp[4];
2123 dndt = ttemp[5];
2124 nm = ttemp[6];
2125
2126 } // if method = d
2127
2128 if(nm <= 0.0)
2129 {
2130 // printf("# error nm %f\n", nm);
2131 satrec.error = 2;
2132 // sgp4fix add return
2133 return false;
2134 }
2135 am = Math.pow((xke / nm), x2o3) * tempa * tempa;
2136 nm = xke / Math.pow(am, 1.5);
2137 em = em - tempe;
2138
2139 // fix tolerance for error recognition
2140 // sgp4fix am is fixed from the previous nm check
2141 if((em >= 1.0) || (em < -0.001)/* || (am < 0.95)*/)
2142 {
2143 // printf("# error em %f\n", em);
2144 satrec.error = 1;
2145 // sgp4fix to return if there is an error in eccentricity
2146 return false;
2147 }
2148 // sgp4fix fix tolerance to avoid a divide by zero
2149 if(em < 1.0e-6)
2150 {
2151 em = 1.0e-6;
2152 }
2153 mm = mm + satrec.no * templ;
2154 xlm = mm + argpm + nodem;
2155 emsq = em * em;
2156 temp = 1.0 - emsq;
2157
2158 nodem = (nodem % twopi);
2159 argpm = (argpm % twopi);
2160 xlm = (xlm % twopi);
2161 mm = ((xlm - argpm - nodem) % twopi);
2162
2163 /* ----------------- compute extra mean quantities ------------- */
2164 sinim = Math.sin(inclm);
2165 cosim = Math.cos(inclm);
2166
2167 /* -------------------- add lunar-solar periodics -------------- */
2168 ep = em;
2169 xincp = inclm;
2170 argpp = argpm;
2171 nodep = nodem;
2172 mp = mm;
2173 sinip = sinim;
2174 cosip = cosim;
2175 if(satrec.method == 'd')
2176 {
2177 //dpper(satrec);
2178 double[] ttemp = dpper(
2179 satrec.e3, satrec.ee2, satrec.peo,
2180 satrec.pgho, satrec.pho, satrec.pinco,
2181 satrec.plo, satrec.se2, satrec.se3,
2182 satrec.sgh2, satrec.sgh3, satrec.sgh4,
2183 satrec.sh2, satrec.sh3, satrec.si2,
2184 satrec.si3, satrec.sl2, satrec.sl3,
2185 satrec.sl4, satrec.t, satrec.xgh2,
2186 satrec.xgh3, satrec.xgh4, satrec.xh2,
2187 satrec.xh3, satrec.xi2, satrec.xi3,
2188 satrec.xl2, satrec.xl3, satrec.xl4,
2189 satrec.zmol, satrec.zmos, satrec.inclo,
2190 'n', ep, xincp, nodep, argpp, mp, satrec.operationmode);
2191 ep = ttemp[0];
2192 xincp = ttemp[1];
2193 nodep = ttemp[2];
2194 argpp = ttemp[3];
2195 mp = ttemp[4];
2196
2197 if(xincp < 0.0)
2198 {
2199 xincp = -xincp;
2200 nodep = nodep + pi;
2201 argpp = argpp - pi;
2202 }
2203 if((ep < 0.0) || (ep > 1.0))
2204 {
2205 // printf("# error ep %f\n", ep);
2206 satrec.error = 3;
2207 // sgp4fix add return
2208 return false;
2209 }
2210 } // if method = d
2211
2212 /* -------------------- long period periodics ------------------ */
2213 if(satrec.method == 'd')
2214 {
2215 sinip = Math.sin(xincp);
2216 cosip = Math.cos(xincp);
2217 satrec.aycof = -0.5 * j3oj2 * sinip;
2218 // sgp4fix for divide by zero for xincp = 180 deg
2219 if(Math.abs(cosip + 1.0) > 1.5e-12)
2220 {
2221 satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip);
2222 }
2223 else
2224 {
2225 satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4;
2226 }
2227 }
2228 axnl = ep * Math.cos(argpp);
2229 temp = 1.0 / (am * (1.0 - ep * ep));
2230 aynl = ep * Math.sin(argpp) + temp * satrec.aycof;
2231 xl = mp + argpp + nodep + temp * satrec.xlcof * axnl;
2232
2233 /* --------------------- solve kepler's equation --------------- */
2234 u = ((xl - nodep) % twopi);
2235 eo1 = u;
2236 tem5 = 9999.9;
2237 ktr = 1;
2238 // sgp4fix for kepler iteration
2239 // the following iteration needs better limits on corrections
2240 while((Math.abs(tem5) >= 1.0e-12) && (ktr <= 10))
2241 {
2242 sineo1 = Math.sin(eo1);
2243 coseo1 = Math.cos(eo1);
2244 tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl;
2245 tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5;
2246 if(Math.abs(tem5) >= 0.95)
2247 {
2248 tem5 = tem5 > 0.0 ? 0.95 : -0.95;
2249 }
2250 eo1 = eo1 + tem5;
2251 ktr = ktr + 1;
2252 }
2253
2254 /* ------------- short period preliminary quantities ----------- */
2255 ecose = axnl * coseo1 + aynl * sineo1;
2256 esine = axnl * sineo1 - aynl * coseo1;
2257 el2 = axnl * axnl + aynl * aynl;
2258 pl = am * (1.0 - el2);
2259 if(pl < 0.0)
2260 {
2261 // printf("# error pl %f\n", pl);
2262 satrec.error = 4;
2263 // sgp4fix add return
2264 return false;
2265 }
2266 else
2267 {
2268 rl = am * (1.0 - ecose);
2269 rdotl = Math.sqrt(am) * esine / rl;
2270 rvdotl = Math.sqrt(pl) / rl;
2271 betal = Math.sqrt(1.0 - el2);
2272 temp = esine / (1.0 + betal);
2273 sinu = am / rl * (sineo1 - aynl - axnl * temp);
2274 cosu = am / rl * (coseo1 - axnl + aynl * temp);
2275 su = Math.atan2(sinu, cosu);
2276 sin2u = (cosu + cosu) * sinu;
2277 cos2u = 1.0 - 2.0 * sinu * sinu;
2278 temp = 1.0 / pl;
2279 temp1 = 0.5 * j2 * temp;
2280 temp2 = temp1 * temp;
2281
2282 /* -------------- update for short period periodics ------------ */
2283 if(satrec.method == 'd')
2284 {
2285 cosisq = cosip * cosip;
2286 satrec.con41 = 3.0 * cosisq - 1.0;
2287 satrec.x1mth2 = 1.0 - cosisq;
2288 satrec.x7thm1 = 7.0 * cosisq - 1.0;
2289 }
2290 mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) +
2291 0.5 * temp1 * satrec.x1mth2 * cos2u;
2292 su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u;
2293 xnode = nodep + 1.5 * temp2 * cosip * sin2u;
2294 xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u;
2295 mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / xke;
2296 rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u +
2297 1.5 * satrec.con41) / xke;
2298
2299 /* --------------------- orientation vectors ------------------- */
2300 sinsu = Math.sin(su);
2301 cossu = Math.cos(su);
2302 snod = Math.sin(xnode);
2303 cnod = Math.cos(xnode);
2304 sini = Math.sin(xinc);
2305 cosi = Math.cos(xinc);
2306 xmx = -snod * cosi;
2307 xmy = cnod * cosi;
2308 ux = xmx * sinsu + cnod * cossu;
2309 uy = xmy * sinsu + snod * cossu;
2310 uz = sini * sinsu;
2311 vx = xmx * cossu - cnod * sinsu;
2312 vy = xmy * cossu - snod * sinsu;
2313 vz = sini * cossu;
2314
2315 /* --------- position and velocity (in km and km/sec) ---------- */
2316 r[0] = (mrt * ux) * radiusearthkm;
2317 r[1] = (mrt * uy) * radiusearthkm;
2318 r[2] = (mrt * uz) * radiusearthkm;
2319 v[0] = (mvt * ux + rvdot * vx) * vkmpersec;
2320 v[1] = (mvt * uy + rvdot * vy) * vkmpersec;
2321 v[2] = (mvt * uz + rvdot * vz) * vkmpersec;
2322 } // if pl > 0
2323
2324 // sgp4fix for decaying satellites
2325 if(mrt < 1.0)
2326 {
2327 // printf("# decay condition %11.6f \n",mrt);
2328 satrec.error = 6;
2329 return false;
2330 }
2331
2332 //#include "debug7.cpp"
2333 return true;
2334 } // end sgp4
2335
2336 /** -----------------------------------------------------------------------------
2337 *
2338 * function gstime
2339 *
2340 * this function finds the greenwich sidereal time.
2341 *
2342 * author : david vallado 719-573-2600 1 mar 2001
2343 *
2344 * inputs description range / units
2345 * jdut1 - julian date in ut1 days from 4713 bc
2346 *
2347 * outputs :
2348 * gstime - greenwich sidereal time 0 to 2pi rad
2349 *
2350 * locals :
2351 * temp - temporary variable for doubles rad
2352 * tut1 - julian centuries from the
2353 * jan 1, 2000 12 h epoch (ut1)
2354 *
2355 * coupling :
2356 * none
2357 *
2358 * references :
2359 * vallado 2004, 191, eq 3-45
2360 * ---------------------------------------------------------------------------
2361 * @param jdut1
2362 * @return
2363 */
2364 public static double gstime(double jdut1)
2365 {
2366 final double twopi = 2.0 * pi;
2367 final double deg2rad = pi / 180.0;
2368 double temp, tut1;
2369
2370 tut1 = (jdut1 - 2451545.0) / 36525.0;
2371 temp = -6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 +
2372 (876600.0 * 3600 + 8640184.812866) * tut1 + 67310.54841; // sec
2373 temp = ((temp * deg2rad / 240.0) % twopi); //360/86400 = 1/240, to deg, to rad
2374
2375 // ------------------------ check quadrants ---------------------
2376 if(temp < 0.0)
2377 {
2378 temp += twopi;
2379 }
2380
2381 return temp;
2382 } // end gstime
2383
2384 /* -----------------------------------------------------------------------------
2385 *
2386 * function getgravconst
2387 *
2388 * this function gets constants for the propagator. note that mu is identified to
2389 * facilitiate comparisons with newer models. the common useage is wgs72.
2390 *
2391 * author : david vallado 719-573-2600 21 jul 2006
2392 *
2393 * inputs :
2394 * whichconst - which set of constants to use wgs72old, wgs72, wgs84
2395 *
2396 * outputs :
2397 * tumin - minutes in one time unit
2398 * mu - earth gravitational parameter
2399 * radiusearthkm - radius of the earth in km
2400 * xke - reciprocal of tumin
2401 * j2, j3, j4 - un-normalized zonal harmonic values
2402 * j3oj2 - j3 divided by j2
2403 *
2404 * locals :
2405 *
2406 * coupling :
2407 * none
2408 *
2409 * references :
2410 * norad spacetrack report #3
2411 * vallado, crawford, hujsak, kelso 2006
2412 ---------------------------------------------------------------------------
2413 * @param whichconst
2414 * @return [tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2]
2415 */
2416 public static double[] getgravconst(
2417 Gravconsttype whichconst // double& tumin,
2418 // double& mu,
2419 // double& radiusearthkm,
2420 // double& xke,
2421 // double& j2,
2422 // double& j3,
2423 // double& j4,
2424 // double& j3oj2
2425 )
2426 {
2427
2428 // return values
2429 double tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2;
2430
2431 switch(whichconst)
2432 {
2433 // -- wgs-72 low precision str#3 constants --
2434 case wgs72old:
2435 mu = 398600.79964; // in km3 / s2
2436 radiusearthkm = 6378.135; // km
2437 xke = 0.0743669161;
2438 tumin = 1.0 / xke;
2439 j2 = 0.001082616;
2440 j3 = -0.00000253881;
2441 j4 = -0.00000165597;
2442 j3oj2 = j3 / j2;
2443 break;
2444 // ------------ wgs-72 constants ------------
2445 case wgs72:
2446 mu = 398600.8; // in km3 / s2
2447 radiusearthkm = 6378.135; // km
2448 xke = 60.0 / Math.sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu);
2449 tumin = 1.0 / xke;
2450 j2 = 0.001082616;
2451 j3 = -0.00000253881;
2452 j4 = -0.00000165597;
2453 j3oj2 = j3 / j2;
2454 break;
2455 case wgs84:
2456 // ------------ wgs-84 constants ------------
2457 mu = 398600.5; // in km3 / s2
2458 radiusearthkm = 6378.137; // km
2459 xke = 60.0 / Math.sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu);
2460 tumin = 1.0 / xke;
2461 j2 = 0.00108262998905;
2462 j3 = -0.00000253215306;
2463 j4 = -0.00000161098761;
2464 j3oj2 = j3 / j2;
2465 break;
2466 default:
2467 System.out.println("unknown gravity option:" + whichconst + ", using wgs84");
2468 // MODIFIED - SHAWN GANO -- Orginal implementation just returned no values!
2469 // ------------ wgs-84 constants ------------
2470 mu = 398600.5; // in km3 / s2
2471 radiusearthkm = 6378.137; // km
2472 xke = 60.0 / Math.sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu);
2473 tumin = 1.0 / xke;
2474 j2 = 0.00108262998905;
2475 j3 = -0.00000253215306;
2476 j4 = -0.00000161098761;
2477 j3oj2 = j3 / j2;
2478 break;
2479 }
2480
2481 return new double[]
2482 {
2483 tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2
2484 };
2485
2486 } // end getgravconst
2487 } // SGP4unit