docs: 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;