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