def sgp4(satrec, tsince, whichconst=None):

     mrt = 0.0
     if whichconst is None:
          whichconst = satrec.whichconst

     """
     /* ------------------ set mathematical constants --------------- */
     // sgp4fix divisor for divide by zero check on inclination
     // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
     // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency
     """
     temp4 =   1.5e-12;
     twopi = 2.0 * pi;
     x2o3  = 2.0 / 3.0;
     #  sgp4fix identify constants and allow alternate values
     tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 = whichconst
     vkmpersec     = radiusearthkm * xke/60.0;

     #  --------------------- clear sgp4 error flag -----------------
     satrec.t     = tsince;
     satrec.error = 0;
     satrec.error_message = None

     #  ------- update for secular gravity and atmospheric drag -----
     xmdf    = satrec.mo + satrec.mdot * satrec.t;
     argpdf  = satrec.argpo + satrec.argpdot * satrec.t;
     nodedf  = satrec.nodeo + satrec.nodedot * satrec.t;
     argpm   = argpdf;
     mm      = xmdf;
     t2      = satrec.t * satrec.t;
     nodem   = nodedf + satrec.nodecf * t2;
     tempa   = 1.0 - satrec.cc1 * satrec.t;
     tempe   = satrec.bstar * satrec.cc4 * satrec.t;
     templ   = satrec.t2cof * t2;

     if satrec.isimp != 1:

         delomg = satrec.omgcof * satrec.t;
         #  sgp4fix use mutliply for speed instead of pow
         delmtemp =  1.0 + satrec.eta * cos(xmdf);
         delm   = satrec.xmcof * \
                  (delmtemp * delmtemp * delmtemp -
                  satrec.delmo);
         temp   = delomg + delm;
         mm     = xmdf + temp;
         argpm  = argpdf - temp;
         t3     = t2 * satrec.t;
         t4     = t3 * satrec.t;
         tempa  = tempa - satrec.d2 * t2 - satrec.d3 * t3 - \
                          satrec.d4 * t4;
         tempe  = tempe + satrec.bstar * satrec.cc5 * (sin(mm) -
                          satrec.sinmao);
         templ  = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof +
                          satrec.t * satrec.t5cof);

     nm    = satrec.no;
     em    = satrec.ecco;
     inclm = satrec.inclo;
     if satrec.method == 'd':

         tc = satrec.t;
         (
             atime, em,    argpm,  inclm, xli,
             mm,    xni,   nodem,  dndt,  nm,
         ) = _dspace(
               satrec.irez,
               satrec.d2201, satrec.d2211, satrec.d3210,
               satrec.d3222, satrec.d4410, satrec.d4422,
               satrec.d5220, satrec.d5232, satrec.d5421,
               satrec.d5433, satrec.dedt,  satrec.del1,
               satrec.del2,  satrec.del3,  satrec.didt,
               satrec.dmdt,  satrec.dnodt, satrec.domdt,
               satrec.argpo, satrec.argpdot, satrec.t, tc,
               satrec.gsto, satrec.xfact, satrec.xlamo,
               satrec.no, satrec.atime,
               em, argpm, inclm, satrec.xli, mm, satrec.xni,
               nodem, nm
             );

     if nm <= 0.0:

         satrec.error_message = ('mean motion {0:f} is less than zero'
                                 .format(nm))
         satrec.error = 2;
         #  sgp4fix add return
         return false, false;

     am = pow((xke / nm),x2o3) * tempa * tempa;
     nm = xke / pow(am, 1.5);
     em = em - tempe;

     #  fix tolerance for error recognition
     #  sgp4fix am is fixed from the previous nm check
     if em >= 1.0 or em < -0.001:  # || (am < 0.95)

         satrec.error_message = ('mean eccentricity {0:f} not within'
                                 ' range 0.0 <= e < 1.0'.format(em))
         satrec.error = 1;
         #  sgp4fix to return if there is an error in eccentricity
         return false, false;

     #  sgp4fix fix tolerance to avoid a divide by zero
     if em < 1.0e-6:
         em  = 1.0e-6;
     mm     = mm + satrec.no * templ;
     xlm    = mm + argpm + nodem;
     emsq   = em * em;
     temp   = 1.0 - emsq;

     nodem  = nodem % twopi if nodem >= 0.0 else -(-nodem % twopi)
     argpm  = argpm % twopi
     xlm    = xlm % twopi
     mm     = (xlm - argpm - nodem) % twopi

     #  ----------------- compute extra mean quantities -------------
     sinim = sin(inclm);
     cosim = cos(inclm);

     #  -------------------- add lunar-solar periodics --------------
     ep     = em;
     xincp  = inclm;
     argpp  = argpm;
     nodep  = nodem;
     mp     = mm;
     sinip  = sinim;
     cosip  = cosim;
     if satrec.method == 'd':

         ep, xincp, nodep, argpp, mp = _dpper(
               satrec, satrec.inclo,
               'n', ep, xincp, nodep, argpp, mp, satrec.afspc_mode
             );
         if xincp < 0.0:

             xincp  = -xincp;
             nodep = nodep + pi;
             argpp  = argpp - pi;

         if ep < 0.0 or ep > 1.0:

             satrec.error_message = ('perturbed eccentricity {0:f} not within'
                                     ' range 0.0 <= e <= 1.0'.format(ep))
             satrec.error = 3;
             #  sgp4fix add return
             return false, false;

     #  -------------------- long period periodics ------------------
     if satrec.method == 'd':

         sinip =  sin(xincp);
         cosip =  cos(xincp);
         satrec.aycof = -0.5*j3oj2*sinip;
         #  sgp4fix for divide by zero for xincp = 180 deg
         if fabs(cosip+1.0) > 1.5e-12:
             satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip);
         else:
             satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4;

     axnl = ep * cos(argpp);
     temp = 1.0 / (am * (1.0 - ep * ep));
     aynl = ep* sin(argpp) + temp * satrec.aycof;
     xl   = mp + argpp + nodep + temp * satrec.xlcof * axnl;

     #  --------------------- solve kepler's equation ---------------
     u    = (xl - nodep) % twopi
     eo1  = u;
     tem5 = 9999.9;
     ktr = 1;
     #    sgp4fix for kepler iteration
     #    the following iteration needs better limits on corrections
     while fabs(tem5) >= 1.0e-12 and ktr <= 10:

         sineo1 = sin(eo1);
         coseo1 = cos(eo1);
         tem5   = 1.0 - coseo1 * axnl - sineo1 * aynl;
         tem5   = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5;
         if fabs(tem5) >= 0.95:
             tem5 = 0.95 if tem5 > 0.0 else -0.95;
         eo1    = eo1 + tem5;
         ktr = ktr + 1;

     #  ------------- short period preliminary quantities -----------
     ecose = axnl*coseo1 + aynl*sineo1;
     esine = axnl*sineo1 - aynl*coseo1;
     el2   = axnl*axnl + aynl*aynl;
     pl    = am*(1.0-el2);
     if pl < 0.0:

         satrec.error_message = ('semilatus rectum {0:f} is less than zero'
                                 .format(pl))
         satrec.error = 4;
         #  sgp4fix add return
         return false, false;

     else:

         rl     = am * (1.0 - ecose);
         rdotl  = sqrt(am) * esine/rl;
         rvdotl = sqrt(pl) / rl;
         betal  = sqrt(1.0 - el2);
         temp   = esine / (1.0 + betal);
         sinu   = am / rl * (sineo1 - aynl - axnl * temp);
         cosu   = am / rl * (coseo1 - axnl + aynl * temp);
         su     = atan2(sinu, cosu);
         sin2u  = (cosu + cosu) * sinu;
         cos2u  = 1.0 - 2.0 * sinu * sinu;
         temp   = 1.0 / pl;
         temp1  = 0.5 * j2 * temp;
         temp2  = temp1 * temp;

         #  -------------- update for short period periodics ------------
         if satrec.method == 'd':

             cosisq                 = cosip * cosip;
             satrec.con41  = 3.0*cosisq - 1.0;
             satrec.x1mth2 = 1.0 - cosisq;
             satrec.x7thm1 = 7.0*cosisq - 1.0;

         mrt   = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + \
                 0.5 * temp1 * satrec.x1mth2 * cos2u;
         su    = su - 0.25 * temp2 * satrec.x7thm1 * sin2u;
         xnode = nodep + 1.5 * temp2 * cosip * sin2u;
         xinc  = xincp + 1.5 * temp2 * cosip * sinip * cos2u;
         mvt   = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / xke;
         rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u +
                 1.5 * satrec.con41) / xke;

         #  --------------------- orientation vectors -------------------
         sinsu =  sin(su);
         cossu =  cos(su);
         snod  =  sin(xnode);
         cnod  =  cos(xnode);
         sini  =  sin(xinc);
         cosi  =  cos(xinc);
         xmx   = -snod * cosi;
         xmy   =  cnod * cosi;
         ux    =  xmx * sinsu + cnod * cossu;
         uy    =  xmy * sinsu + snod * cossu;
         uz    =  sini * sinsu;
         vx    =  xmx * cossu - cnod * sinsu;
         vy    =  xmy * cossu - snod * sinsu;
         vz    =  sini * cossu;

         #  --------- position and velocity (in km and km/sec) ----------
         _mr = mrt * radiusearthkm
         r = (_mr * ux, _mr * uy, _mr * uz)
         v = ((mvt * ux + rvdot * vx) * vkmpersec,
              (mvt * uy + rvdot * vy) * vkmpersec,
              (mvt * uz + rvdot * vz) * vkmpersec)

     #  sgp4fix for decaying satellites
     if mrt < 1.0:

         satrec.error_message = ('mrt {0:f} is less than 1.0 indicating'
                                 ' the satellite has decayed'.format(mrt))
         satrec.error = 6;
         return false, false;

     return r, v;
/* ------------------ set mathematical constants --------------- */
     // sgp4fix divisor for divide by zero check on inclination
     // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
     // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency