Notebook source code: examples/kepler-numba/kepler.ipynb
Run the notebook yourself on binder Binder badge

Kepler (minimum time orbit transfer) - Numba

CNES TAS

Minimum time control of the Kepler equation (CNES / TAS / Inria / CNRS collaboration):

\[t_f \to \min,\]
\[\ddot{q} = -\mu\frac{q}{|q|^3}+\frac{u}{m}\,,\quad t \in [0,t_f],\]
\[\dot{m} = -\beta|u|,\quad |u| \leq T_{\mathrm{max}}.\]

Fixed initial and final Keplerian orbits (free final longitude).

Thumbnail

Initializations

[1]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
from nutopy import nle
from nutopy import tools
from nutopy import ocp
from numba import njit

ctmax = (3600**2) / 1e6                                             # Conversion from Newtons
mass0 = 1500.                                                       # Initial mass of the spacecraft
beta = 1.42e-02                                                     # Engine specific impulsion
mu = 5165.8620912                                                   # Earth gravitation constant
t0 = 0.                                                             # Initial time (final time is free)
x0 = np.array([ 11.625, 0.75, 0., 6.12e-02, 0., 3.14159265358979 ]) # Initial state (fixed initial longitude)
xf_fixed = np.array([ 42.165, 0., 0., 0., 0. ])                     # Final state (free final longitude)

# tmax = 60 Newtons
tmax = ctmax * 60.; tf = 15.2055; p0 = -np.array([ .361266, 22.2412, 7.87736, 0., 0., -5.90802 ]); N = 1000

# tmax = 6 Newtons
#tmax = ctmax * 6.; tf = 1.32e2; p0 = -np.array([ -4.743728539366440e+00, -7.171314869854240e+01, -2.750468309804530e+00, 4.505679923365745e+01, -3.026794475592510e+00, 2.248091067047670e+00 ]); N = 1000

# tmax = 0.7 Newtons
#tmax = ctmax * 0.7; tf = 1.210000000000000e+03; p0 = -np.array([ -2.215319700438820e+01, -4.347109477345140e+01, 9.613188807286992e-01, 3.181800985503019e+02, -2.307236094862410e+00, -5.797863110671591e-01 ]); N = 5000

# tmax = 0.14 Newtons
#tmax = ctmax * 0.14; tf = 6.08e3; p0 = -np.array([ -1.234155379067110e+02, -6.207170881591489e+02, 5.742554220129187e-01, 1.629324243017332e+03, -2.373935935351530e+00, -2.854066853269850e-01 ]) ; N = 5000

p0 = p0 / np.linalg.norm(p0) # Normalization |p0|=1 for free final time
y = np.hstack((p0, tf)) # initial guess, y = (p0, tf)

Hamiltonian (python + Numba Just In Time compilation)

[2]:
@njit
def hfun(t, x, p, tmax, mass0, beta, mu):
    pa = x[0]
    ex = x[1]
    ey = x[2]
    hx = x[3]
    hy = x[4]
    lg = x[5]

    p1 = p[0]
    p2 = p[1]
    p3 = p[2]
    p4 = p[3]
    p5 = p[4]
    p6 = p[5]

    pdm = np.sqrt(pa/mu)
    cl = np.cos(lg)
    sl = np.sin(lg)
    w = 1.0+ex*cl+ey*sl
    pdmw = pdm/w
    zz = hx*sl-hy*cl
    uh = (1+hx**2+hy**2)/2.0

    f06 = w**2/(pa*pdm)
    h0  = p6*f06

    f12 = pdm *   sl
    f13 = pdm * (-cl)
    h1  = p2*f12 + p3*f13

    f21 = pdm * 2.0 * pa / w
    f22 = pdm * (cl+(ex+cl)/w)
    f23 = pdm * (sl+(ey+sl)/w)
    h2  = p1*f21 + p2*f22 + p3*f23

    f32 = pdmw * (-zz*ey)
    f33 = pdmw * zz*ex
    f34 = pdmw * uh * cl
    f35 = pdmw * uh * sl
    f36 = pdmw * zz
    h3  = p2*f32 + p3*f33 + p4*f34 + p5*f35 + p6*f36

    mass = mass0 - beta*tmax*t

    psi = np.sqrt(h1**2 + h2**2 + h3**2)

    h = -1.0 + h0 + (tmax/mass) * psi
[3]:
@njit
def hfun_d(t, x, xd, p, pd, tmax, mass0, beta, mu):
    pad = xd[0]
    pa = x[0]
    exd = xd[1]
    ex = x[1]
    eyd = xd[2]
    ey = x[2]
    hxd = xd[3]
    hx = x[3]
    hyd = xd[4]
    hy = x[4]
    lgd = xd[5]
    lg = x[5]
    p1d = pd[0]
    p1 = p[0]
    p2d = pd[1]
    p2 = p[1]
    p3d = pd[2]
    p3 = p[2]
    p4d = pd[3]
    p4 = p[3]
    p5d = pd[4]
    p5 = p[4]
    p6d = pd[5]
    p6 = p[5]
    temp = np.sqrt(pa/mu)
    if pa/mu == 0.0:
        pdmd = 0.0
    else:
        pdmd = pad/(2.0*temp*mu)
    pdm = temp
    cld = -(np.sin(lg)*lgd)
    cl = np.cos(lg)
    sld = np.cos(lg)*lgd
    sl = np.sin(lg)
    wd = cl*exd + ex*cld + sl*eyd + ey*sld
    w = 1.0 + ex*cl + ey*sl
    pdmwd = (pdmd-pdm*wd/w)/w
    pdmw = pdm/w
    zzd = sl*hxd + hx*sld - cl*hyd - hy*cld
    zz = hx*sl - hy*cl
    uhd = (2*hx*hxd+2*hy*hyd)/2.0
    uh = (1+hx**2+hy**2)/2.0
    temp = w*w/(pa*pdm)
    f06d = (2*w*wd-temp*(pdm*pad+pa*pdmd))/(pa*pdm)
    f06 = temp
    h0d = f06*p6d + p6*f06d
    h0 = p6*f06
    f12d = sl*pdmd + pdm*sld
    f12 = pdm*sl
    f13d = -(cl*pdmd+pdm*cld)
    f13 = pdm*(-cl)
    h1d = f12*p2d + p2*f12d + f13*p3d + p3*f13d
    h1 = p2*f12 + p3*f13
    temp = pdm*pa/w
    f21d = 2.0*(pa*pdmd+pdm*pad-temp*wd)/w
    f21 = 2.0*temp
    temp = (ex+cl)/w
    f22d = (cl+temp)*pdmd + pdm*(cld+(exd+cld-temp*wd)/w)
    f22 = pdm*(cl+temp)
    temp = (ey+sl)/w
    f23d = (sl+temp)*pdmd + pdm*(sld+(eyd+sld-temp*wd)/w)
    f23 = pdm*(sl+temp)
    h2d = f21*p1d + p1*f21d + f22*p2d + p2*f22d + f23*p3d + p3*f23d
    h2 = p1*f21 + p2*f22 + p3*f23
    f32d = -(ey*(zz*pdmwd+pdmw*zzd)+pdmw*zz*eyd)
    f32 = pdmw*(-(zz*ey))
    f33d = ex*(zz*pdmwd+pdmw*zzd) + pdmw*zz*exd
    f33 = pdmw*zz*ex
    f34d = cl*(uh*pdmwd+pdmw*uhd) + pdmw*uh*cld
    f34 = pdmw*uh*cl
    f35d = sl*(uh*pdmwd+pdmw*uhd) + pdmw*uh*sld
    f35 = pdmw*uh*sl
    f36d = zz*pdmwd + pdmw*zzd
    f36 = pdmw*zz
    h3d = f32*p2d + p2*f32d + f33*p3d + p3*f33d + f34*p4d + p4*f34d + f35*p5d + p5*f35d + f36*p6d + p6*f36d
    h3 = p2*f32 + p3*f33 + p4*f34 + p5*f35 + p6*f36
    mass = mass0 - beta*tmax*t
    arg1d = 2*h1*h1d + 2*h2*h2d + 2*h3*h3d
    arg1 = h1**2 + h2**2 + h3**2
    temp = np.sqrt(arg1)
    if arg1 == 0.0:
        psid = 0.0
    else:
        psid = arg1d/(2.0*temp)
    psi = temp
    hd = h0d + tmax*psid/mass
    h = -1.0 + h0 + tmax/mass*psi
    return h, hd
[4]:
@njit
def hfun_d_d(t, x, xd, xd0, p, pd, pd0, tmax, mass0, beta, mu):
    pad = xd[0]
    pad0 = xd0[0]
    pa = x[0]
    exd = xd[1]
    exd0 = xd0[1]
    ex = x[1]
    eyd = xd[2]
    eyd0 = xd0[2]
    ey = x[2]
    hxd = xd[3]
    hxd0 = xd0[3]
    hx = x[3]
    hyd = xd[4]
    hyd0 = xd0[4]
    hy = x[4]
    lgd = xd[5]
    lgd0 = xd0[5]
    lg = x[5]
    p1d = pd[0]
    p1d0 = pd0[0]
    p1 = p[0]
    p2d = pd[1]
    p2d0 = pd0[1]
    p2 = p[1]
    p3d = pd[2]
    p3d0 = pd0[2]
    p3 = p[2]
    p4d = pd[3]
    p4d0 = pd0[3]
    p4 = p[3]
    p5d = pd[4]
    p5d0 = pd0[4]
    p5 = p[4]
    p6d = pd[5]
    p6d0 = pd0[5]
    p6 = p[5]
    temp0 = np.sqrt(pa/mu)
    if (pa/mu == 0.0):
        tempd = 0.0
    else:
        tempd = pad0/(2.0*temp0*mu)
    temp = temp0
    if (pa/mu == 0.0):
        pdmd = 0.0
        pdmdd = 0.0
    else:
        temp0 = pad/(2.0*mu*temp)
        pdmdd = -(temp0*tempd/temp)
        pdmd = temp0
    pdmd0 = tempd
    pdm = temp
    cldd = -(lgd*np.cos(lg)*lgd0)
    cld = -(np.sin(lg)*lgd)
    cld0 = -(np.sin(lg)*lgd0)
    cl = np.cos(lg)
    sldd = -(lgd*np.sin(lg)*lgd0)
    sld = np.cos(lg)*lgd
    sld0 = np.cos(lg)*lgd0
    sl = np.sin(lg)
    wdd = exd*cld0 + cld*exd0 + ex*cldd + eyd*sld0 + sld*eyd0 + ey*sldd
    wd = cl*exd + ex*cld + sl*eyd + ey*sld
    wd0 = cl*exd0 + ex*cld0 + sl*eyd0 + ey*sld0
    w = 1.0 + ex*cl + ey*sl
    temp0 = pdm*wd/w
    temp1 = (pdmd-temp0)/w
    pdmwdd = (pdmdd-(wd*pdmd0+pdm*wdd-temp0*wd0)/w-temp1*wd0)/w
    pdmwd = temp1
    pdmwd0 = (pdmd0-pdm*wd0/w)/w
    pdmw = pdm/w
    zzdd = hxd*sld0 + sld*hxd0 + hx*sldd - hyd*cld0 - cld*hyd0 - hy*cldd
    zzd = sl*hxd + hx*sld - cl*hyd - hy*cld
    zzd0 = sl*hxd0 + hx*sld0 - cl*hyd0 - hy*cld0
    zz = hx*sl - hy*cl
    uhdd = (hxd*2*hxd0+hyd*2*hyd0)/2.0
    uhd = (2*hx*hxd+2*hy*hyd)/2.0
    uhd0 = (2*hx*hxd0+2*hy*hyd0)/2.0
    uh = (1+hx**2+hy**2)/2.0
    temp1 = w*w/(pa*pdm)
    tempd = (2*w*wd0-temp1*(pdm*pad0+pa*pdmd0))/(pa*pdm)
    temp = temp1
    temp1 = pad*pdm + pa*pdmd
    temp0 = (2*w*wd-temp*temp1)/(pa*pdm)
    f06dd = (wd*2*wd0+2*w*wdd-temp1*tempd-temp*(pad*pdmd0+pdmd*pad0+pa*pdmdd)-temp0*(pdm*pad0+pa*pdmd0))/(pa*pdm)
    f06d = temp0
    f06d0 = tempd
    f06 = temp
    h0dd = p6d*f06d0 + f06d*p6d0 + p6*f06dd
    h0d = f06*p6d + p6*f06d
    h0 = p6*f06
    f12dd = pdmd*sld0 + sl*pdmdd + sld*pdmd0 + pdm*sldd
    f12d = sl*pdmd + pdm*sld
    f12d0 = sl*pdmd0 + pdm*sld0
    f12 = pdm*sl
    f13dd = -(pdmd*cld0) - cl*pdmdd - cld*pdmd0 - pdm*cldd
    f13d = -(cl*pdmd+pdm*cld)
    f13d0 = -(cl*pdmd0+pdm*cld0)
    f13 = pdm*(-cl)
    h1dd = p2d*f12d0 + f12d*p2d0 + p2*f12dd + p3d*f13d0 + f13d*p3d0 + p3*f13dd
    h1d = f12*p2d + p2*f12d + f13*p3d + p3*f13d
    h1d0 = f12*p2d0 + p2*f12d0 + f13*p3d0 + p3*f13d0
    h1 = p2*f12 + p3*f13
    temp1 = pdm*pa/w
    tempd = (pa*pdmd0+pdm*pad0-temp1*wd0)/w
    temp = temp1
    temp1 = (pa*pdmd+pad*pdm-temp*wd)/w
    f21dd = 2.0*(pdmd*pad0+pa*pdmdd+pad*pdmd0-wd*tempd-temp*wdd-temp1*wd0)/w
    f21d = 2.0*temp1
    f21d0 = 2.0*tempd
    f21 = 2.0*temp
    temp1 = (ex+cl)/w
    tempd = (exd0+cld0-temp1*wd0)/w
    temp = temp1
    temp1 = (exd+cld-temp*wd)/w
    f22dd = pdmd*(cld0+tempd) + (cl+temp)*pdmdd + (cld+temp1)*pdmd0 + pdm*(cldd+(cldd-wd*tempd-temp*wdd-temp1*wd0)/w)
    f22d = (cl+temp)*pdmd + pdm*(cld+temp1)
    f22d0 = (cl+temp)*pdmd0 + pdm*(cld0+tempd)
    f22 = pdm*(cl+temp)
    temp1 = (ey+sl)/w
    tempd = (eyd0+sld0-temp1*wd0)/w
    temp = temp1
    temp1 = (eyd+sld-temp*wd)/w
    f23dd = pdmd*(sld0+tempd) + (sl+temp)*pdmdd + (sld+temp1)*pdmd0 + pdm*(sldd+(sldd-wd*tempd-temp*wdd-temp1*wd0)/w)
    f23d = (sl+temp)*pdmd + pdm*(sld+temp1)
    f23d0 = (sl+temp)*pdmd0 + pdm*(sld0+tempd)
    f23 = pdm*(sl+temp)
    h2dd = p1d*f21d0 + f21d*p1d0 + p1*f21dd + p2d*f22d0 + f22d*p2d0 + p2*f22dd + p3d*f23d0 + f23d*p3d0 + p3*f23dd
    h2d = f21*p1d + p1*f21d + f22*p2d + p2*f22d + f23*p3d + p3*f23d
    h2d0 = f21*p1d0 + p1*f21d0 + f22*p2d0 + p2*f22d0 + f23*p3d0 + p3*f23d0
    h2 = p1*f21 + p2*f22 + p3*f23
    temp1 = zz*pdmwd + pdmw*zzd
    f32dd = -(temp1*eyd0) - ey*(pdmwd*zzd0+zz*pdmwdd+zzd*pdmwd0+pdmw*zzdd) - eyd*(zz*pdmwd0+pdmw*zzd0)
    f32d = -(ey*temp1) - eyd*(pdmw*zz)
    f32d0 = -(ey*(zz*pdmwd0+pdmw*zzd0)+pdmw*zz*eyd0)
    f32 = pdmw*(-(zz*ey))
    temp1 = zz*pdmwd + pdmw*zzd
    f33dd = temp1*exd0 + ex*(pdmwd*zzd0+zz*pdmwdd+zzd*pdmwd0+pdmw*zzdd) + exd*(zz*pdmwd0+pdmw*zzd0)
    f33d = ex*temp1 + exd*(pdmw*zz)
    f33d0 = ex*(zz*pdmwd0+pdmw*zzd0) + pdmw*zz*exd0
    f33 = pdmw*zz*ex
    temp1 = uh*pdmwd + pdmw*uhd
    f34dd = temp1*cld0 + cl*(pdmwd*uhd0+uh*pdmwdd+uhd*pdmwd0+pdmw*uhdd) + cld*(uh*pdmwd0+pdmw*uhd0) + pdmw*uh*cldd
    f34d = cl*temp1 + pdmw*uh*cld
    f34d0 = cl*(uh*pdmwd0+pdmw*uhd0) + pdmw*uh*cld0
    f34 = pdmw*uh*cl
    temp1 = uh*pdmwd + pdmw*uhd
    f35dd = temp1*sld0 + sl*(pdmwd*uhd0+uh*pdmwdd+uhd*pdmwd0+pdmw*uhdd) + sld*(uh*pdmwd0+pdmw*uhd0) + pdmw*uh*sldd
    f35d = sl*temp1 + pdmw*uh*sld
    f35d0 = sl*(uh*pdmwd0+pdmw*uhd0) + pdmw*uh*sld0
    f35 = pdmw*uh*sl
    f36dd = pdmwd*zzd0 + zz*pdmwdd + zzd*pdmwd0 + pdmw*zzdd
    f36d = zz*pdmwd + pdmw*zzd
    f36d0 = zz*pdmwd0 + pdmw*zzd0
    f36 = pdmw*zz
    h3dd = p2d*f32d0 + f32d*p2d0 + p2*f32dd + p3d*f33d0 + f33d*p3d0 + p3*f33dd + p4d*f34d0 + f34d*p4d0 + p4*f34dd + p5d*f35d0 + f35d*p5d0 + p5*f35dd + p6d*f36d0 + f36d*p6d0 + p6*f36dd
    h3d = f32*p2d + p2*f32d + f33*p3d + p3*f33d + f34*p4d + p4*f34d + f35*p5d + p5*f35d + f36*p6d + p6*f36d
    h3d0 = f32*p2d0 + p2*f32d0 + f33*p3d0 + p3*f33d0 + f34*p4d0 + p4*f34d0 + f35*p5d0 + p5*f35d0 + f36*p6d0 + p6*f36d0
    h3 = p2*f32 + p3*f33 + p4*f34 + p5*f35 + p6*f36
    mass = mass0 - beta*tmax*t
    arg1dd = h1d*2*h1d0 + 2*h1*h1dd + h2d*2*h2d0 + 2*h2*h2dd + h3d*2*h3d0 + 2*h3*h3dd
    arg1d = 2*h1*h1d + 2*h2*h2d + 2*h3*h3d
    arg1d0 = 2*h1*h1d0 + 2*h2*h2d0 + 2*h3*h3d0
    arg1 = h1**2 + h2**2 + h3**2
    temp1 = np.sqrt(arg1)
    if (arg1 == 0.0):
        tempd = 0.0
    else:
        tempd = arg1d0/(2.0*temp1)
    temp = temp1
    if (arg1 == 0.0):
        psid = 0.0
        psidd = 0.0
    else:
        temp1 = arg1d/(2.0*temp)
        psidd = (arg1dd-temp1*2.0*tempd)/(2.0*temp)
        psid = temp1
    psi = temp
    hdd = h0dd + tmax*psidd/mass
    hd = h0d + tmax*psid/mass
    h = -1.0 + h0 + tmax/mass*psi
    return h, hd, hdd
[5]:
hfun = tools.tensorize(hfun_d, hfun_d_d, tvars=(2, 3), full=True)(hfun)
h = ocp.Hamiltonian(hfun)
f = ocp.Flow(h)

Shooting function

[6]:
def dshoot(t0, dt0, x0, dx0, p0, dp0, tf, dtf, next=False):
    (xf, dxf), (pf, dpf) = f((t0, dt0), (x0, dx0), (p0, dp0), (tf, dtf), tmax, mass0, beta, mu)
    s = np.zeros(7) # code duplication and full=True
    s[0:5] = xf[0:5] - xf_fixed
    s[5] = pf[5] # free final longitude
    s[6] = p0[0]**2 + p0[1]**2 + p0[2]**2 + p0[3]**2 + p0[4]**2 + p0[5]**2 - 1.
    ds = np.zeros(7)
    ds[0:5] = dxf[0:5]
    ds[5] = dpf[5]
    ds[6] = 2*p0[0]*dp0[0] + 2*p0[1]*dp0[1] + 2*p0[2]*dp0[2] + 2*p0[3]*dp0[3] + 2*p0[4]*dp0[4] + 2*p0[5]*dp0[5]
    if not next: return s, ds
    else: return s, ds, ((tf, dtf), (xf, dxf), (pf, dpf), None)

@tools.vectorize(vvars=(3,))
@tools.vectorize(vvars=(4,), next=True)
@tools.tensorize(dshoot, full=True)
def shoot(t0, x0, p0, tf, next=False):
    """s = shoot(t0, x0, p0, tf)

    Shooting function associated with h
    """
    xf, pf = f(t0, x0, p0, tf, tmax, mass0, beta, mu)
    s = np.zeros(7)
    s[0:5] = xf[0:5] - xf_fixed
    s[5] = pf[5] # free final longitude
    s[6] = p0[0]**2 + p0[1]**2 + p0[2]**2 + p0[3]**2 + p0[4]**2 + p0[5]**2 - 1.
    if not next: return s
    else: return s, (tf, xf, pf, None)

Solve

[7]:
dfoo = lambda y, dy: shoot(t0, x0, (y[:-1], dy[:-1]), (y[-1], dy[-1]))
foo = lambda y: shoot(t0, x0, y[:-1], y[-1])
foo = tools.tensorize(dfoo, full=True)(foo)

nleopt = nle.Options(SolverMethod='hybrj', Display='on', TolX=1e-8)
et = time.time(); sol = nle.solve(foo, y, df=foo, options=nleopt); y_sol = sol.x; et = time.time() - et
print('Elapsed time:', et)
print('y_sol =', y_sol)
print('foo =', foo(y_sol))

     Calls  |f(x)|                 |x|

         1  2.292126134511967e+00  1.523834735954001e+01
         2  4.556611736187258e-01  1.475855901478838e+01
         3  2.151564089092056e-02  1.485241942477452e+01
         4  8.256206372631771e-03  1.483036174390080e+01
         5  5.169276888587317e-03  1.483414821919238e+01
         6  5.596292416643664e-04  1.483439823807463e+01
         7  4.977547273842370e-05  1.483424598086740e+01
         8  3.628885599208353e-05  1.483411574622270e+01
         9  1.988047917363645e-05  1.483410845271092e+01
        10  4.925723971613692e-06  1.483410894452168e+01
        11  1.045416659513755e-07  1.483410881453910e+01
        12  7.952332083287182e-09  1.483410884353815e+01

 Results of the nle solver method:

 xsol    =  [-0.01618971 -0.90848895 -0.32525986 -0.09447217  0.03432296  0.24184433
 14.80036436]
 f(xsol) =  [ 7.79323983e-09 -1.26483179e-09  2.13029667e-10  8.43600263e-10
 -1.94929158e-10  1.56334936e-11 -3.31533134e-10]
 nfev    =  12
 njev    =  1
 status  =  1
 success =  True

 Successfully completed: relative error between two consecutive iterates is at most TolX.

Elapsed time: 4.460953712463379
y_sol = [-0.01618971 -0.90848895 -0.32525986 -0.09447217  0.03432296  0.24184433
 14.80036436]
foo = [ 7.79323983e-09 -1.26483179e-09  2.13029667e-10  8.43600263e-10
 -1.94929158e-10  1.56334936e-11 -3.31533134e-10]

Plots

[8]:
p0 = y_sol[:-1]
tf = y_sol[-1]
tspan = list(np.linspace(0, tf, N+1))
xf, pf = np.array( f(t0, x0, p0, tspan, tmax, mass0, beta, mu) )
P  = xf[:, 0]
ex = xf[:, 1]
ey = xf[:, 2]
hx = xf[:, 3]
hy = xf[:, 4]
L  = xf[:, 5]
cL = np.cos(L)
sL = np.sin(L)
W  = 1+ex*cL+ey*sL
Z  = hx*sL-hy*cL
C  = 1+hx**2+hy**2
q  = np.zeros((N+1, 3))
q[:, 0] = P*( (1+hx**2-hy**2)*cL + 2*hx*hy*sL ) / (C*W)
q[:, 1] = P*( (1-hx**2+hy**2)*sL + 2*hx*hy*cL ) / (C*W)
q[:, 2] = 2*P*Z / (C*W)

%matplotlib widget
plt.rcParams['legend.fontsize'] = 20
plt.rcParams['figure.figsize'] = (10, 10)
fig1 = plt.figure()
ax = fig1.gca(projection='3d')
ax.set_xlim3d(-60, 60) # ax.axis('equal') not supported
ax.set_ylim3d(-60, 60)
ax.set_zlim3d(-40, 40)
u, v = np.mgrid[ 0:2*np.pi:20j, 0:np.pi:10j ]
r = 6.378 # Earth radius (in Mm)
x1 = r*np.cos(u)*np.sin(v)
x2 = r*np.sin(u)*np.sin(v)
x3 = r*np.cos(v)
ax.plot_wireframe(x1, x2, x3, color='b')
ax.plot(q[:, 0], q[:, 1], q[:, 2], label='Orbit transfer')
#ax.quiver(q[:, 0], q[:, 1], q[:, 2], u, v, w, length=0.1, normalize=True)
ax.legend()
/Users/ocots/opt/miniconda3/envs/gallery/lib/python3.7/site-packages/mpl_toolkits/mplot3d/art3d.py:304: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  self._segments3d = np.asanyarray(segments)
[8]:
<matplotlib.legend.Legend at 0x7fc796622090>
[ ]: