Notebook source code:
examples/solarsail/solarsail-simple-version.ipynb
Run the notebook yourself on binder
Solar sail (maximum orbit change)¶
Maximize the dispacement of the solar sail in the tangent manifold of its orbital elements over one orbit:
\[\delta I' = \mathbf{G}(I,M) \, u,\quad u(M) \in U,\]
\[\delta I(0) = 0, \quad \delta I(2\pi) \parallel \mathbf{d},\]
\[\left( \delta I (2\pi) \, | \, \mathbf{d} \right) \longrightarrow \max.\]
Fixed initial Keplerian orbit and direction of the displacement \(\mathbf{d}\).
Initializations¶
[1]:
import numpy as np
import scipy
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [15, 15]
plt.rcParams['figure.dpi'] = 200 # 200 e.g. is really fine, but slower
from mpl_toolkits.mplot3d import Axes3D
import time
from nutopy import nle
from nutopy import tools
from nutopy import ocp
from scipy.linalg import null_space
from gve import *
from os import system
# Definition of optical parameters
rho = 0.6 # Specular reflection coefficient
a = 1. - rho # Absorption coefficient
s = 0.94 # Diffuse reflection coefficient
eps_f = 0.05 # Emissivity front coeff 0.05
eps_b = 0.55 # Emissivity back coeff 0.55
Bf = 0.79 # Front Lambertian coeff 0.79
Bb = 0.55 # Back Lambertian coeff 0.55
eps = (Bf * eps_f - Bb * eps_b) / (eps_b + eps_f)
#opt_p = np.array([rho, s, eps_f, eps_b, Bf, Bb, eps])
b = np.array([1. - rho * s, 2. * rho * s, Bf * rho * (1. - s) + (1. - rho) * eps ])
# Parameters of the problem
cont = 1. # 0 for triangular cone, 1 for real drop
mu = 1.
sDir = np.array([0., 0., -1.])
# Initial guess of pI taken from Matlab Convex Programming
pI = np.array([-0.1492, -0.2744, -0.1004, -0.2594, 1.0000])
M1_guess = 0.8762
M2_guess = 5.3809
I1_guess = np.array([-0.6552, 0.1258, 0.4852, -0.0240, 0.0263])
I2_guess = np.array([-0.6543, 0.1304, 0.4841, -0.0274, 0.0267])
z = np.hstack((pI, M1_guess, M2_guess, I1_guess, I2_guess))
def deg2rad(x):
y = x * np.pi / 180.
return y
def rad2deg(x):
y = x * 180. / np.pi
return y
# Initial state
I = np.array([deg2rad(10), deg2rad(20), deg2rad(30), 1, 0.6]) # Initial state
d = np.array([ 0., 0., 0., 0., 1.]) # Direction of the displacement
I0 = np.array([0., 0., 0., 0., 0]) # delta I (0)
M0 = 0
Mf = 2*np.pi
# Real cone force definition
cBeta = (- b[0] * b[2] - 2. * b[1] * b[2] + \
np.sqrt(b[0] ** 2 * b[2] ** 2 - 4. * b[0] * b[2] ** 2 * b[1] + 8. * b[0] ** 2 * b[1] ** 2 + 4. * b[1] ** 3 * b[0])) / \
(4. * b[0] * b[1] + 2. * b[1] ** 2);
sBeta = np.sqrt(1 - cBeta ** 2);
fs = b[0] * cBeta + (b[1] * cBeta ** 2 + b[2] * cBeta) * cBeta;
fperp = (b[1] * cBeta ** 2 + b[2] * cBeta) * sBeta;
fCone = np.array([fs, fperp])
pars = np.hstack((mu, I, fCone, b))
# Orthogonal vector family
d_orth = scipy.linalg.null_space([d]) # Orthogonal space of b
d_orth = d_orth.transpose()
[2]:
#!pygmentize hfun.f90
Hamiltonian¶
[3]:
!python -m numpy.f2py -c hfun.f90 hfun.pyf -m hfun > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d.f90 hfun_d.pyf -m hfun_d > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d_d.f90 hfun_d_d.pyf -m hfun_d_d > /dev/null 2>&1
!python -m numpy.f2py -c hfunalpha_da.f90 -m hfunalpha_da > /dev/null 2>&1
[4]:
from hfun import hfun
from hfun_d import hfun_d
from hfun_d_d import hfun_d_d
from hfunalpha_da import hfunalpha_da as dhfunalpha
hfun_u = lambda M, q, p, pars, cont : hfun.hfun_u(M, q, p, pars, cont)
ufun = lambda M, q, p, cont : hfun.control(M, q, p, pars, cont)
theta2alpha = lambda theta, b : hfun.theta2alpha(theta, b)
gveeci = lambda M, I, mu : hfun.gveeci(M, I, mu)
dhfun_u = lambda M, dM, q, dq, p, dp, pars, cont, dcont : hfun_d.hfun_u_d(M, dM, q, dq, p, dp, pars, cont, dcont)
d2hfun_u = lambda M, dM, d2M, q, dq, d2q, p, dp, d2p, pars, cont, dcont, d2cont : hfun_d_d.hfun_u_d_d(M, dM, d2M, q, dq, d2q, p, dp, d2p, pars, cont, dcont, d2cont)
[5]:
hfun_u(M2_guess, I, 2*pI, pars, cont)
[5]:
-1.641289891501918e-05
[6]:
dhfun_u(M2_guess, 4, I1_guess, np.array([0,0.0,0.0,0.0,0.0]), pI, np.array([0.0,0.0,0.0,0.0,0.0]), pars, cont,0)
[6]:
(-8.20644945750959e-06, 0.1205097240305964)
[7]:
def f_aux(M,I, pdeltaI):
[fx, fy, fz] = gveeci(M, I, mu)
pIG = np.array([0.0, 0.0, 0.0])
pIG[0] = np.dot(fx[0:5],pdeltaI)
pIG[1] = np.dot(fy[0:5],pdeltaI)
pIG[2] = np.dot(fz[0:5],pdeltaI)
theta = np.arccos(np.dot(pIG,sDir)/np.linalg.norm(pIG))
alpha = theta2alpha(theta, b)
return alpha
[8]:
def dhfun_0(M, dM, q, dq, p, dp, pars, cont, dcont):
h = 0.0
dh = 0.0
return h, dh
def d2hfun_0(M, dM, d2M, q, dq, d2q, p, dp, d2p, pars, cont, dcont, d2cont):
h = 0.0
dh = 0.0
d2h = 0.0
return h, dh, d2h
def hfun_0(M, q, p, pars, cont):
h = 0.0
return h
[9]:
hfun_u = tools.tensorize(dhfun_u, d2hfun_u, tvars=(1, 2, 3, 5), full=True)(hfun_u)
hfun_0 = tools.tensorize(dhfun_0, d2hfun_0, tvars=(1, 2, 3, 5), full=True)(hfun_0)
ufun = tools.vectorize(vvars=(1,2,3))(ufun)
Hu = ocp.Hamiltonian(hfun_u)
H0 = ocp.Hamiltonian(hfun_0)
fu = ocp.Flow(Hu)
f0 = ocp.Flow(H0)
Shooting function¶
[10]:
#Initial condition
deltaI0 = np.zeros(5)
M0 = 0.
Mf = 2. * np.pi
[11]:
def dshoot(z, dz, pars, d, d_orth, cont):
#def dshoot(z, dz, pars, d, d_orth, cont, next=False):
I0 = np.zeros(5)
M0 = 0.
Mf = 2. * np.pi
pI = z[0 : 5]
M1 = z[5]
M2 = z[6]
I1 = z[7 : 12]
I2 = z[12 : 17]
dpI = dz[0 : 5]
dM1 = dz[5]
dM2 = dz[6]
dI1 = dz[7 : 12]
dI2 = dz[12 : 17]
dM0 = 0.
dcont = 0.
dI0 = 0.
#----------------------------------
(I1sol, dI1sol), (pI1, dpI1) = fu(M0, I0, (pI, dpI), (M1, dM1), pars, cont)
(I2sol, dI2sol), (pI2, dpI2) = f0((M1, dM1), (I1, dI1), (pI, dpI), (M2, dM2), pars, cont)
(If, dIf), (pIf, dpIf) = fu((M2, dM2), (I2, dI2), (pI, dpI), Mf, pars, cont)
#----------------------------------
s = np.zeros(17)
s[0] = pIf[0] * d[0] + pIf[1] * d[1] + pIf[2] * d[2] + pIf[3] * d[3] + pIf[4] * d[4] - np.linalg.norm(d)**2
s[1 : 5] = np.dot(d_orth, If)
s[7 : 12] = I1 - I1sol
s[12 : 17] = I2 - I2sol
ds = np.zeros(17)
ds[0] = dpIf[0] * d[0] + dpIf[1] * d[1] + dpIf[2] * d[2] + dpIf[3] * d[3] + dpIf[4] * d[4]
ds[1 : 5] = np.dot(d_orth, dIf)
ds[7 : 12] = dI1 - dI1sol
ds[12 : 17] = dI2 - dI2sol
s[5], ds[5] = Hu((M1, dM1), (I1, dI1), (pI, dpI), pars, cont) #- H0((M1, dM1), (I1, dI1), (pI, dpI), pars, cont)
s[6], ds[6] = Hu((M2, dM2), (I2, dI2), (pI, dpI), pars, cont) #- H0((M2, dM2), (I2, dI2), (pI, dpI), pars, cont)
# on impose pas la continuité de pI car il est constant
# il faut imposer la continuité de l'etat deltaI entre les switchs entre 7 et 17
#if not next: return s, ds
#else: return s, ds, (Mf, If, pIf, None)
return s, ds
#@tools.vectorize(vvars=(3,))
#@tools.vectorize(vvars=(4,), next=True)
@tools.tensorize(dshoot, tvars=(1,), full=True)
def shoot(z, pars, d, d_orth, cont):
#def shoot(z, pars, d, d_orth, cont, next=False):
I0 = np.zeros(5)
M0 = 0.
Mf = 2. * np.pi
pI = z[0 : 5]
M1 = z[5]
M2 = z[6]
I1 = z[7 : 12]
I2 = z[12 : 17]
#----------------------------------
I1sol, pI1 = fu(M0, I0, pI, M1, pars, cont)
I2sol, pI2 = f0(M1, I1, pI, M2, pars, cont)
If, pIf = fu(M2, I2, pI, Mf, pars, cont)
print(I1sol,I2sol, If)
#----------------------------------
s = np.zeros(17)
s[0] = pIf[0] * d[0] + pIf[1] * d[1] + pIf[2] * d[2] + pIf[3] * d[3] + pIf[4] * d[4] - np.linalg.norm(d)**2
s[1 : 5] = np.dot(d_orth, If)
# theta1 = atan2(pI1*sPerp/(pI1*sDir))
# theta2 = atan2(pI2*sPerp/(pI2*sDir))
# s[5] = theta1 - cone_alpha - deg2rad(90.)
# s[6] = theta2 - cone_alpha - deg2rad(90.)
s[5] = Hu(M1, I1, pI, pars, cont)
s[6] = Hu(M2, I2, pI, pars, cont)
s[7 : 12] = I1 - I1sol
s[12 : 17] = I2 - I2sol
# on impose pas la continuité de pI car il est constant
# il faut imposer la continuité de l'etat deltaI entre les switchs
#if not next: return s
#else: return s, (Mf, If, pIf, None)
return s
[12]:
@tools.vectorize(vvars =(1,))
def control(M, I, pI, cont):
sol = np.array([0., 0., 0.])
if type(M)==float or type(M)==int:
if Hu(M, I, pI, pars, cont) < 0.:
sol = np.array([0., 0., 0.])
else:
sol = ufun(M, I, pI, cont)
else:
for ii in range(len(M)):
if Hu(M[ii], I, pI, pars, cont) < 0.:
sol = np.vstack((sol, np.array([0., 0., 0.])))
else:
sol = np.vstack((sol, ufun(M[ii], I, pI, cont)))
sol = np.delete(sol, 0, 0)
return sol
Solve¶
[13]:
dfoo = lambda z, dz, cont: shoot((z, dz), pars, d, d_orth, (cont, 0.))
foo = lambda z, cont: shoot(z, pars, d, d_orth, cont)
#foo = tools.tensorize(dfoo, full=True)(foo)
#foo = tools.tensorize(dfoo, tvars=(1,), full=True)(foo)
[14]:
nleopt = nle.Options(SolverMethod='hybrj', Display='on', TolX=1e-8)
#et = time.time(); sol = nle.solve(foo, z, df=foo, options=nleopt, args=cont);
et = time.time(); sol = nle.solve(foo, z, options=nleopt, args=cont);
#et = time.time(); sol = nle.solve(foo, z, options=nleopt);
et = time.time() - et
print('Elapsed time:', et)
z_sol = sol.x;
print('z_sol =', z_sol)
print('foo =', foo(z_sol, cont))
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
Calls |f(x)| |x|
1 5.670902558167515e-01 5.679741568944841e+00
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703461 0.14090175 0.82166179 -0.04146046 0.02516804] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006258 -0.0932223 -0.14265169 -0.0221303 0.09145796]
[-1.02703469 0.14090177 0.82166186 -0.04146048 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006244 -0.09322228 -0.14265159 -0.02213027 0.09145796]
[-1.02703474 0.14090177 0.82166193 -0.04146045 0.02516804] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006235 -0.09322228 -0.14265149 -0.02213029 0.09145795]
[-1.02703467 0.14090176 0.82166187 -0.04146042 0.02516805] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006243 -0.09322228 -0.14265161 -0.02213024 0.09145796]
[-1.02703467 0.14090176 0.82166186 -0.04146046 0.02516804] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.1426516 -0.02213027 0.09145796]
[-1.02703469 0.14090177 0.82166186 -0.04146048 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.2000624 -0.09322229 -0.14265157 -0.02213026 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.65519999 0.1258 0.4852 -0.024 0.0263 ] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.12580001 0.4852 -0.024 0.0263 ] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.48520001 -0.024 0.0263 ] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.02399999 0.0263 ] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.02630001] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006246 -0.09322229 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322227 -0.1426516 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.14265158 -0.02213028 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.1426516 -0.02213026 0.09145795]
[-1.02703468 0.14090176 0.82166185 -0.04146047 0.02516803] [-0.6552 0.1258 0.4852 -0.024 0.0263] [ 0.20006245 -0.09322229 -0.1426516 -0.02213028 0.09145797]
[-1.09052814 0.21585281 0.85341875 -0.03763988 0.03884065] [-1.09087596 0.20867219 0.85386362 -0.03454811 0.03928921] [ 0.03748703 0.00426459 -0.03306304 -0.00499942 0.12804089]
2 5.153381521200057e-02 5.777224618971843e+00
[-1.08222595 0.20376867 0.85455454 -0.03209697 0.03875793] [-1.08390747 0.20652489 0.85567477 -0.0330346 0.0387136 ] [-0.0180777 -0.00148415 0.01483277 0.00163212 0.12360018]
3 2.392243221026179e-02 5.789399146522509e+00
[-1.08417522 0.20695523 0.85424342 -0.03314101 0.03892934] [-1.08416573 0.20700725 0.85424067 -0.03311168 0.03894241] [-9.06131950e-05 -1.56138241e-06 7.06452324e-05 1.70825924e-05
1.24966652e-01]
4 1.484563069365445e-04 5.784666203493183e+00
[-1.08418939 0.20700145 0.85424367 -0.0331261 0.03894175] [-1.08418384 0.20699912 0.85423975 -0.03312276 0.03894282] [ 1.57250904e-05 -8.26933608e-06 -1.12564100e-06 -8.75763867e-06
1.24974819e-01]
5 2.158898188935950e-05 5.784754695761182e+00
[-1.08419262 0.20700755 0.8542465 -0.03312352 0.03894365] [-1.08419394 0.20700726 0.85424744 -0.03312329 0.03894373] [-3.64839749e-07 -8.90324879e-07 8.28349788e-07 -9.50895724e-07
1.24975958e-01]
6 2.319768562870607e-06 5.784768319366612e+00
[-1.08419293 0.20700832 0.85424685 -0.03312304 0.03894394] [-1.0841931 0.20700826 0.85424697 -0.03312302 0.03894395] [ 2.46094328e-07 -3.40444006e-08 -1.80594269e-07 -7.94559305e-08
1.24976185e-01]
7 3.897942659182793e-07 5.784769176219384e+00
[-1.08419292 0.20700833 0.85424687 -0.033123 0.03894396] [-1.08419292 0.20700833 0.85424688 -0.033123 0.03894396] [ 1.14084066e-08 -3.31875872e-09 -5.15658186e-09 -1.68575496e-09
1.24976168e-01]
8 1.389707715245373e-08 5.784769347693150e+00
[-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396] [-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396] [-1.41021240e-10 6.42491147e-11 -2.31293457e-12 3.80533834e-10
1.24976168e-01]
9 5.585999511799720e-10 5.784769358763763e+00
Results of the nle solver method:
xsol = [-0.15486778 -0.23740693 -0.11720524 -0.24128161 1. 0.99543827
5.23631576 -1.08419292 0.20700834 0.85424687 -0.033123 0.03894396
-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396]
f(xsol) = [ 0.00000000e+00 6.42491147e-11 -2.31293457e-12 3.80533834e-10
1.41021240e-10 -8.01714806e-11 -7.61704311e-11 -6.84392543e-11
4.87117013e-11 3.71664921e-11 -3.26306357e-10 -1.26744899e-10
0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00]
nfev = 9
njev = 1
status = 1
success = True
Successfully completed: relative error between two consecutive iterates is at most TolX.
Elapsed time: 0.8945662975311279
z_sol = [-0.15486778 -0.23740693 -0.11720524 -0.24128161 1. 0.99543827
5.23631576 -1.08419292 0.20700834 0.85424687 -0.033123 0.03894396
-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396]
[-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396] [-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396] [-1.41021240e-10 6.42491147e-11 -2.31293457e-12 3.80533834e-10
1.24976168e-01]
foo = [ 0.00000000e+00 6.42491147e-11 -2.31293457e-12 3.80533834e-10
1.41021240e-10 -8.01714806e-11 -7.61704311e-11 -6.84392543e-11
4.87117013e-11 3.71664921e-11 -3.26306357e-10 -1.26744899e-10
0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00]
[15]:
pI = z_sol[0 : 5]
M1 = z_sol[5]
M2 = z_sol[6]
M_sol = np.array([M1,M2])
I1 = z_sol[7 : 12]
I2 = z_sol[12 : 17]
u = control(M0, I0, pI, cont)
Function of continuous integration¶
[16]:
def plotSolutionBSB(pI, M1, M2, Mf, sDir, cont, nOrb):
N = 500
M_vect = np.array([])
g1 = np.array([])
g2 = np.array([])
g3 = np.array([])
a = np.array([])
e = np.array([])
dvect = np.array([])
I = np.array([0., 0., 0., 0., 0])
I0 = np.array([0., 0., 0., 0., 0])
M0 = 0
M0_2_M1 = np.linspace(M0, M1, N+1)
M1_2_M2 = np.linspace(M1, M2, N+1)
M2_2_Mf = np.linspace(M2, Mf, N+1)
u_vect = np.array([0., 0., 0.])
for ii in range(nOrb):
Mspan1 = list((ii)*2*np.pi*np.ones(N+1)+M0_2_M1)
Mspan2 = list((ii)*2*np.pi*np.ones(N+1)+M1_2_M2)
Mspanf = list((ii)*2*np.pi*np.ones(N+1)+M2_2_Mf)
I1sol, pI1sol = fu(M0, I0, pI, Mspan1, pars, cont) # on [ 0, t1]
I2sol, pI2sol = f0(Mspan1[-1], I1sol[-1], pI, Mspan2, pars, cont) # on [t1, t2]
Ifsol, pIfsol = fu(Mspan2[-1], I2sol[-1], pI, Mspanf, pars, cont) # on [t2, tf]
M_vect = np.hstack((M_vect, Mspan1, Mspan2, Mspanf))
qd = np.zeros(5)
pd = np.zeros(5)
for kk in range(N+1):
val = I1sol[kk]
g1 = np.hstack((g1, val[0]))
g2 = np.hstack((g2, val[1]))
g3 = np.hstack((g3, val[2]))
a = np.hstack((a, val[3]))
e = np.hstack((e, val[4]))
alpha = f_aux(Mspan1[kk], pars[1:6], pI1sol[kk])
dvect = np.hstack((dvect, dhfunalpha(Mspan1[kk],0.0,I1sol[kk], qd,pI,pd,alpha,1.0,pars,cont)[1]))
for kk in range(N+1):
val = I2sol[kk]
g1 = np.hstack((g1, val[0]))
g2 = np.hstack((g2, val[1]))
g3 = np.hstack((g3, val[2]))
a = np.hstack((a, val[3]))
e = np.hstack((e, val[4]))
dvect = np.hstack((dvect,0))
for kk in range(N+1):
val = Ifsol[kk]
g1 = np.hstack((g1, val[0]))
g2 = np.hstack((g2, val[1]))
g3 = np.hstack((g3, val[2]))
a = np.hstack((a, val[3]))
e = np.hstack((e, val[4]))
alpha = f_aux(Mspanf[kk], pars[1:6], pIfsol[kk])
dvect = np.hstack((dvect, dhfunalpha(Mspanf[kk],0.0,Ifsol[kk], qd,pI,pd,alpha,1.0,pars,cont)[1]))
I0 = Ifsol[-1]
M0 = Mspanf[-1]
u_vect = control(M_vect, I, pI, cont)
u_s = np.dot(u_vect,sDir)
u_1 = u_vect[:,0]
u_2 = u_vect[:,1]
u_3 = u_vect[:,2]
ang1 = np.arccos(u_s/np.linalg.norm(u_vect,axis=1))
where_are_NaNs = np.isnan(ang1)
ang1[where_are_NaNs] = 0
ang1 = rad2deg(ang1)
#for i in range(np.size(ang1)-1):
#if(ang1[i]==0 and ang1[i+1]!=0):
#print(deg2rad(ang1[i+1]))
#ang2 = np.sqrt(np.divide(u_3 ** 2, u_2 ** 2))
#where_are_NaNs = np.isnan(ang2)
#ang2[where_are_NaNs] = 0
ang2 = np.arctan2(u_1 ** 2, u_2 ** 2)
ang2 = rad2deg(ang2)
#M_vect = M_vect/(2*np.pi)
M_vect = M_vect / 2 / np.pi * 360
fig = plt.figure(); plt.tight_layout()
ax1 = fig.add_subplot(421);
ax1.plot(M_vect, g1);
ax2 = fig.add_subplot(422);
ax2.plot(M_vect, g2);
ax3 = fig.add_subplot(423);
ax3.plot(M_vect, g3);
ax4 = fig.add_subplot(424);
ax4.plot(M_vect, a);
ax5 = fig.add_subplot(425);
ax5.plot(M_vect, e);
ax6 = fig.add_subplot(426);
ax6.plot(M_vect, ang1, label='cone angle');
ax6.plot(M_vect, ang2, label='clock angle');
ax7 = fig.add_subplot(427);
ax7.plot(M_vect, dvect);
label_name = 'Satellite longitude [deg]'
ax1.set_xlabel(label_name); ax1.set_ylabel('$\gamma_1 \epsilon^{-1}$'); ax1.axhline(0, color='k')
ax2.set_xlabel(label_name); ax2.set_ylabel('$\gamma_2 \epsilon^{-1}$'); ax2.axhline(0, color='k')
ax3.set_xlabel(label_name); ax3.set_ylabel('$\gamma_3 \epsilon^{-1}$'); ax3.axhline(0, color='k')
ax4.set_xlabel(label_name); ax4.set_ylabel('$a \epsilon^{-1}$'); ax4.axhline(0, color='k')
ax5.set_xlabel(label_name); ax5.set_ylabel('$e \epsilon^{-1}$'); ax5.axhline(0, color='k')
ax6.set_xlabel(label_name); ax6.set_ylabel('$Control \ angle $');
ax6.legend()
ax7.set_xlabel(label_name); ax7.set_ylabel('dhfunalpha');
plt.savefig('figures/plot.png')
Plots¶
[17]:
pI = z_sol[0 : 5]
M1 = z_sol[5]
M2 = z_sol[6]
I1 = z_sol[7 : 12]
I2 = z_sol[12 : 17]
print(z_sol)
If, pI = fu(M2, I2, pI, 2*np.pi, pars, cont)
[-0.15486778 -0.23740693 -0.11720524 -0.24128161 1. 0.99543827
5.23631576 -1.08419292 0.20700834 0.85424687 -0.033123 0.03894396
-1.08419292 0.20700834 0.85424687 -0.033123 0.03894396]
One orbit¶
Trajectory¶
[18]:
plotSolutionBSB(pI, M1, M2, Mf, sDir, cont, 1)
/net/home/c/ct705827/.conda/envs/ct-gallery/lib/python3.7/site-packages/ipykernel_launcher.py:71: RuntimeWarning: invalid value encountered in true_divide
Multiple orbits¶
[19]:
plotSolutionBSB(pI, M1, M2, Mf, sDir, cont, 3)
/net/home/c/ct705827/.conda/envs/ct-gallery/lib/python3.7/site-packages/ipykernel_launcher.py:71: RuntimeWarning: invalid value encountered in true_divide
[ ]: