001/*
002 * $Id: SGP4unit.java,v 1.2 2011/03/24 16:06:33 davep Exp $
003 *
004 * This file is part of McIDAS-V
005 *
006 * Copyright 2007-2011
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
086package edu.wisc.ssec.mcidasv.data.adde.sgp4;
087
088/**
089 *
090 * @author Shawn E. Gano, shawn@gano.name
091 */
092public 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