Commit 641e1870 authored by CAILLAU Jean-Baptiste's avatar CAILLAU Jean-Baptiste
Browse files

Merge branch 'feature/solarsail' into 'develop'

Feature/solarsail

See merge request !288
parents a53f640d 55080b74
......@@ -27,7 +27,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
......@@ -49,13 +49,13 @@
"xf_fixed = np.array([ 42.165, 0., 0., 0., 0. ]) # Final state (free final longitude)\n",
"\n",
"# tmax = 60 Newtons\n",
"#tmax = ctmax * 60.; tf = 15.2055; p0 = -np.array([ .361266, 22.2412, 7.87736, 0., 0., -5.90802 ]); N = 1000\n",
"tmax = ctmax * 60.; tf = 15.2055; p0 = -np.array([ .361266, 22.2412, 7.87736, 0., 0., -5.90802 ]); N = 1000\n",
"\n",
"# tmax = 6 Newtons\n",
"#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\n",
"\n",
"# tmax = 0.7 Newtons\n",
"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\n",
"#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\n",
"\n",
"# tmax = 0.14 Newtons\n",
"#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\n",
......@@ -73,7 +73,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
......@@ -129,7 +129,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
......@@ -227,7 +227,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
......@@ -417,7 +417,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
......@@ -435,7 +435,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
......@@ -478,7 +478,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 8,
"metadata": {},
"outputs": [
{
......@@ -488,25 +488,25 @@
"\n",
" Calls |f(x)| |x|\n",
" \n",
" 1 1.073164376904338e+00 1.210000413223070e+03\n",
" 2 1.434466432323507e-01 1.213474308936708e+03\n",
" 3 4.359026937651052e-02 1.215428931328697e+03\n",
" 4 5.821599360781245e-03 1.214866946053342e+03\n",
" 5 1.549151998131196e-03 1.214502522548834e+03\n",
" 6 5.328741165581899e-04 1.214573443587823e+03\n",
" 7 2.366668442307806e-04 1.214586014482633e+03\n",
" 8 6.447228903961488e-05 1.214593835332666e+03\n",
" 9 7.042908938241817e-06 1.214592759663435e+03\n",
" 10 1.280661104130505e-06 1.214592611776815e+03\n",
" 11 7.690550348931297e-08 1.214592632124916e+03\n",
" 12 3.965613410327104e-09 1.214592633915170e+03\n",
" 1 2.292126134511974e+00 1.523834735954001e+01\n",
" 2 4.556611736187402e-01 1.475855901478837e+01\n",
" 3 2.151564089092733e-02 1.485241942477453e+01\n",
" 4 8.256206372659084e-03 1.483036174390080e+01\n",
" 5 5.169276888558879e-03 1.483414821919238e+01\n",
" 6 5.596292416639707e-04 1.483439823807462e+01\n",
" 7 4.977547276492455e-05 1.483424598086740e+01\n",
" 8 3.628885615536121e-05 1.483411574622269e+01\n",
" 9 1.988047925176352e-05 1.483410845271090e+01\n",
" 10 4.925724021356949e-06 1.483410894452168e+01\n",
" 11 1.045417722601297e-07 1.483410881453910e+01\n",
" 12 7.952339019384582e-09 1.483410884353816e+01\n",
"\n",
" Results of the nle solver method:\n",
"\n",
" xsol = [ 6.93341643e-02 5.16187599e-01 -2.31725360e-04 -8.53653346e-01\n",
" 4.36325335e-03 -9.32349207e-05 1.21459222e+03]\n",
" f(xsol) = [ 3.65812980e-09 -1.30123645e-10 1.43819251e-09 -1.23351486e-11\n",
" -7.67939625e-11 -5.09919820e-11 -5.00196551e-10]\n",
" xsol = [-0.01618971 -0.90848895 -0.32525986 -0.09447217 0.03432296 0.24184433\n",
" 14.80036436]\n",
" f(xsol) = [ 7.79324694e-09 -1.26483227e-09 2.13028215e-10 8.43599709e-10\n",
" -1.94929269e-10 1.56336359e-11 -3.31532912e-10]\n",
" nfev = 12\n",
" njev = 1\n",
" status = 1\n",
......
import numpy as np
#############################################################################################
## Gauss Variational Equations (GVE)
def gvef(OE, mu):
#
# Gives projection of gauss' form of variational equations on r,t,n vectors
# r = r, t = theta, n = h
#
# [R, T, N] = gve(OE, mu)
# OE = [Omega, i, w, a, e, M]
#
#Omega = OE[0]
i = OE[1]
w = OE[2]
a = OE[3]
e = OE[4]
f = OE[5]
cF = np.cos(f);
sF = np.sin(f);
p = a * (1. - e**2)
r = p / (1. + e * cF)
b = a * np.sqrt(1. - e**2)
n = np.sqrt(mu / a**3)
h = n * a * b
theta = w + f
dOmega = np.array([ 0., 0., r * np.sin(theta) / h / np.sin(i)])
di = np.array([ 0., 0., r * np.cos(theta) / h])
dw = np.array([ - p * cF / h / e, (p + r) * sF / h / e, - r * np.sin(theta) * np.cos(i) / h / np.sin(i)])
da = np.array([2. * a**2 * e * sF / h, 2. * a**2 * p / h / r, 0.])
de = np.array([ p * sF / h, ((p + r) * cF + r * e) / h, 0.])
dPhi = np.array([ p * cF - 2. * r * e, - (p + r) * sF, 0.]) * b / a / h / e
R = np.array([dOmega[0], di[0], dw[0], da[0], de[0], dPhi[0]])
T = np.array([dOmega[1], di[1], dw[1], da[1], de[1], dPhi[1]])
N = np.array([dOmega[2], di[2], dw[2], da[2], de[2], dPhi[2]])
return R, T, N, f
#############################################################################################
## GVE in ECI frame
def gveecif(OE, mu):
#
# GVE with perturbation expressed in the ECI frame, i.e.,
#
# OEdot = [R, T, N] * uLVLH = [R, T, N] * Reci2lvlh * uECI = [Fx, Fy, Fz] * uECI
#
R, T, N, f = gvef(OE, mu)
Omega = OE[0]
i = OE[1]
w = OE[2]
a = OE[3]
e = OE[4]
theta = w + f
sOm = np.sin(Omega)
cOm = np.cos(Omega)
sI = np.sin(i)
cI = np.cos(i)
sTh = np.sin(theta)
cTh = np.cos(theta)
# GVE in ECI frame
Fx = R * ((- sOm * cI * sTh + cOm * cTh)) + T * ((- sOm * cI * cTh - cOm * sTh)) + N * ((sOm * sI) )
Fy = R * (( cOm * cI * sTh + sOm * cTh)) + T * (( cOm * cI * cTh - sOm * sTh)) + N * (- (cOm * sI) )
Fz = R * (sI * sTh) + T * (sI * cTh) + N * cI
return Fx, Fy, Fz
......@@ -13,7 +13,23 @@ python module hfun
double precision, intent(in) :: M, q(5), p(5), pars(11), cont
double precision, intent(out) :: u(3)
end subroutine control
subroutine kepler(ecc, M, f)
implicit none
double precision, intent(in) :: ecc, M
double precision, intent(out) :: f
end subroutine kepler
subroutine gveeci(M, I, mu, Fx, Fy, Fz)
implicit none
double precision, intent(in) :: M, I(5), mu
double precision, intent(out) :: Fx(6), Fy(6), Fz(6)
end subroutine
subroutine theta2alpha(tTheta, b, alpha)
double precision, intent(in) :: tTheta, b(3)
double precision, intent(out) :: alpha
end subroutine theta2alpha
end module hfun
......
......@@ -7,14 +7,14 @@ MODULE HFUN_D
CONTAINS
! Differentiation of hfun_u in forward (tangent) mode (with options fixinterface):
! variations of useful results: h
! with respect to varying inputs: p q
! RW status of diff variables: h:out p:in q:in
! with respect to varying inputs: cont m p q
! RW status of diff variables: h:out cont:in m:in p:in q:in
!#######################################################################
!## Hamiltonian
SUBROUTINE HFUN_U_D(m, q, qd, p, pd, pars, cont, h, hd)
SUBROUTINE HFUN_U_D(m, md, q, qd, p, pd, pars, cont, contd, h, hd)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: m, q(5), p(5), pars(11), cont
DOUBLE PRECISION, INTENT(IN) :: qd(5), pd(5)
DOUBLE PRECISION, INTENT(IN) :: md, qd(5), pd(5), contd
DOUBLE PRECISION, INTENT(OUT) :: h
DOUBLE PRECISION, INTENT(OUT) :: hd
!local variables
......@@ -35,17 +35,15 @@ CONTAINS
! Direction of the Sun
sdir = (/0.0d0, 0.0d0, -1.0d0/)
! GVEs
CALL GVEECI(m, i, mu, fx, fy, fz)
CALL GVEECI_D(m, md, i, mu, fx, fxd, fy, fyd, fz, fzd)
! Optimal control (note that cont is the continuation parameter)
fxd = 0.D0
CALL DOT_D(5, p, pd, fx, fxd, pig(1), pigd(1))
fyd = 0.D0
CALL DOT_D(5, p, pd, fy, fyd, pig(2), pigd(2))
fzd = 0.D0
CALL DOT_D(5, p, pd, fz, fzd, pig(3), pigd(3))
CALL PMPCONE_D(pig, pigd, sdir, fcone, ucone, uconed)
CALL PMPSAIL_D(pig, pigd, sdir, bsail, usail, usaild)
uconed = (1.0d0-cont)*uconed + cont*usaild
uconed = uconed*(1.0d0-cont) - ucone*contd + usaild*cont + usail*&
& contd
ucone = ucone*(1.0d0-cont) + usail*cont
! Hamiltonian
CALL DOT_D(3, pig, pigd, ucone, uconed, h, hd)
......@@ -394,6 +392,52 @@ CONTAINS
sperp = sperp/result1
u = fcone(1)*sdir + fcone(2)*sperp
END SUBROUTINE PMPCONE
! Differentiation of gveeci in forward (tangent) mode (with options fixinterface):
! variations of useful results: fx fy fz
! with respect to varying inputs: m
! ######################################################################
! ## GVE in ECI frame
SUBROUTINE GVEECI_D(m, md, i, mu, fx, fxd, fy, fyd, fz, fzd)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: m, i(5), mu
DOUBLE PRECISION, INTENT(IN) :: md
DOUBLE PRECISION, INTENT(OUT) :: fx(6), fy(6), fz(6)
DOUBLE PRECISION, INTENT(OUT) :: fxd(6), fyd(6), fzd(6)
! Local variables
DOUBLE PRECISION :: r(6), t(6), n(6), f, om, inc, w, a, e, theta, &
& som, com, si, ci, sth, cth
DOUBLE PRECISION :: rd(6), td(6), nd(6), fd, thetad, sthd, cthd
INTRINSIC SIN
INTRINSIC COS
! Orbital elements
om = i(1)
inc = i(2)
w = i(3)
a = i(4)
e = i(5)
! GVE in LVLH frame
CALL GVELVLH_D(m, md, i, mu, r, rd, t, td, n, nd, f, fd)
thetad = fd
theta = w + f
som = SIN(om)
com = COS(om)
si = SIN(inc)
ci = COS(inc)
sthd = thetad*COS(theta)
sth = SIN(theta)
cthd = -(thetad*SIN(theta))
cth = COS(theta)
! GVE in ECI frame
fxd = rd*(-(som*ci*sth)+com*cth) + r*(com*cthd-som*ci*sthd) + td*(-(&
& som*ci*cth)-com*sth) + t*(-(som*ci*cthd)-com*sthd) + som*si*nd
fx = r*(-(som*ci*sth)+com*cth) + t*(-(som*ci*cth)-com*sth) + n*(som*&
& si)
fyd = rd*(com*ci*sth+som*cth) + r*(com*ci*sthd+som*cthd) + td*(com*&
& ci*cth-som*sth) + t*(com*ci*cthd-som*sthd) - com*si*nd
fy = r*(com*ci*sth+som*cth) + t*(com*ci*cth-som*sth) + n*(-(com*si))
fzd = si*(rd*sth+r*sthd) + si*(td*cth+t*cthd) + ci*nd
fz = r*(si*sth) + t*(si*cth) + n*ci
END SUBROUTINE GVEECI_D
! ######################################################################
! ## GVE in ECI frame
SUBROUTINE GVEECI(m, i, mu, fx, fy, fz)
......@@ -426,6 +470,91 @@ CONTAINS
fy = r*(com*ci*sth+som*cth) + t*(com*ci*cth-som*sth) + n*(-(com*si))
fz = r*(si*sth) + t*(si*cth) + n*ci
END SUBROUTINE GVEECI
! Differentiation of gvelvlh in forward (tangent) mode (with options fixinterface):
! variations of useful results: f n r t
! with respect to varying inputs: m
! ######################################################################
! ## GVE in LVLH frame
SUBROUTINE GVELVLH_D(m, md, i, mu, r, rd, t, td, n, nd, f, fd)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: m, i(5), mu
DOUBLE PRECISION, INTENT(IN) :: md
DOUBLE PRECISION, INTENT(OUT) :: r(6), t(6), n(6), f
DOUBLE PRECISION, INTENT(OUT) :: rd(6), td(6), nd(6), fd
! Local variables
DOUBLE PRECISION :: inc, w, a, e, cf, sf, p, rad, b, nnorb, h, theta
DOUBLE PRECISION :: cfd, sfd, radd, thetad
INTRINSIC COS
INTRINSIC SIN
INTRINSIC SQRT
DOUBLE PRECISION :: arg1
DOUBLE PRECISION :: result1
! Orbital elements
! Omega = I(1)
inc = i(2)
w = i(3)
a = i(4)
e = i(5)
! Solving Kepler's equation
CALL KEPLER_D(e, m, md, f, fd)
! Useful variables
cfd = -(fd*SIN(f))
cf = COS(f)
sfd = fd*COS(f)
sf = SIN(f)
p = a*(1.0d0-e**2)
radd = -(p*e*cfd/(1.0d0+e*cf)**2)
rad = p/(1.0d0+e*cf)
arg1 = 1.0d0 - e**2
result1 = SQRT(arg1)
b = a*result1
arg1 = mu/a**3
nnorb = SQRT(arg1)
h = nnorb*a*b
thetad = fd
theta = w + f
! GVEs
rd(1) = 0.D0
r(1) = 0.0d0
td(1) = 0.D0
t(1) = 0.0d0
nd = 0.D0
nd(1) = (radd*SIN(theta)+rad*thetad*COS(theta))/h/SIN(inc)
n(1) = rad*SIN(theta)/h/SIN(inc)
rd(2) = 0.D0
r(2) = 0.0d0
td(2) = 0.D0
t(2) = 0.0d0
nd(2) = (radd*COS(theta)-rad*thetad*SIN(theta))/h
n(2) = rad*COS(theta)/h
rd = 0.D0
rd(3) = -(p*cfd/h/e)
r(3) = -(p*cf/h/e)
td = 0.D0
td(3) = (radd*sf+(p+rad)*sfd)/h/e
t(3) = (p+rad)*sf/h/e
nd(3) = -(COS(inc)*(radd*SIN(theta)+rad*thetad*COS(theta))/h/SIN(inc&
& ))
n(3) = -(rad*SIN(theta)*COS(inc)/h/SIN(inc))
rd(4) = 2.0d0*a**2*e*sfd/h
r(4) = 2.0d0*a**2*e*sf/h
td(4) = -(2.0d0*a**2*p*radd/h/rad**2)
t(4) = 2.0d0*a**2*p/h/rad
nd(4) = 0.D0
n(4) = 0.0d0
rd(5) = p*sfd/h
r(5) = p*sf/h
td(5) = (radd*cf+(p+rad)*cfd+e*radd)/h
t(5) = ((p+rad)*cf+rad*e)/h
nd(5) = 0.D0
n(5) = 0.0d0
rd(6) = p*cfd - 2.0d0*e*b*radd/a/h/e
r(6) = p*cf - 2.0d0*rad*e*b/a/h/e
td(6) = -(b*(radd*sf+(p+rad)*sfd)/a/h/e)
t(6) = -((p+rad)*sf*b/a/h/e)
nd(6) = 0.D0
n(6) = 0.0d0
END SUBROUTINE GVELVLH_D
! ######################################################################
! ## GVE in LVLH frame
SUBROUTINE GVELVLH(m, i, mu, r, t, n, f)
......@@ -479,6 +608,54 @@ CONTAINS
t(6) = -((p+rad)*sf*b/a/h/e)
n(6) = 0.0d0
END SUBROUTINE GVELVLH
! Differentiation of kepler in forward (tangent) mode (with options fixinterface):
! variations of useful results: f
! with respect to varying inputs: m
! ######################################################################
! ## Kepler's equation
SUBROUTINE KEPLER_D(ecc, m, md, f, fd)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: ecc, m
DOUBLE PRECISION, INTENT(IN) :: md
DOUBLE PRECISION, INTENT(OUT) :: f
DOUBLE PRECISION, INTENT(OUT) :: fd
! Local variables
INTEGER :: nmax, ii
DOUBLE PRECISION :: e, k, dk
DOUBLE PRECISION :: ed, kd, dkd
INTRINSIC SIN
INTRINSIC COS
INTRINSIC SQRT
INTRINSIC ATAN2
DOUBLE PRECISION :: result1
DOUBLE PRECISION :: arg1
DOUBLE PRECISION :: arg1d
DOUBLE PRECISION :: result2
DOUBLE PRECISION :: arg2
DOUBLE PRECISION :: arg2d
! Set parameters
nmax = 20
! Newton-Rapshon
ed = md
e = m
DO ii=1,nmax
kd = ed - ecc*ed*COS(e) - md
k = e - ecc*SIN(e) - m
dkd = ecc*ed*SIN(e)
dk = 1.0d0 - ecc*COS(e)
ed = ed - (kd*dk-k*dkd)/dk**2
e = e - k/dk
END DO
! True Anomaly
result1 = SQRT(1.0d0 + ecc)
arg1d = result1*ed*COS(e/2.0d0)/2.0d0
arg1 = result1*SIN(e/2.0d0)
result2 = SQRT(1.0d0 - ecc)
arg2d = -(result2*ed*SIN(e/2.0d0)/2.0d0)
arg2 = result2*COS(e/2.0d0)
fd = 2.0d0*(arg1d*arg2-arg2d*arg1)/(arg1**2+arg2**2)
f = 2.0d0*ATAN2(arg1, arg2)
END SUBROUTINE KEPLER_D
! ######################################################################
! ## Kepler's equation
SUBROUTINE KEPLER(ecc, m, f)
......@@ -548,4 +725,3 @@ CONTAINS
END DO
END SUBROUTINE DOT
END MODULE HFUN_D
......@@ -4,8 +4,8 @@ python module hfun_d
module hfun_d
subroutine hfun_u_d(m, q, qd, p, pd, pars, cont, h, hd)
double precision, intent(in) :: m, q(5), p(5), pars(11), cont, qd(5), pd(5)
subroutine hfun_u_d(m, md, q, qd, p, pd, pars, cont, contd, h, hd)
double precision, intent(in) :: m, q(5), p(5), pars(11), cont, md, qd(5), pd(5), contd
double precision, intent(out) :: h, hd
end subroutine hfun_u_d
......
This diff is collapsed.
......@@ -4,10 +4,10 @@ python module hfun_d_d
module hfun_d_d
subroutine hfun_u_d_d(m, q, qd0, qd, p, pd0, pd, pars, cont, h, hd, hdd)
subroutine hfun_u_d_d(m, md0, md, q, qd0, qd, p, pd0, pd, pars, cont, contd0, contd, h, hd, hdd)
double precision, intent(in) :: m, q(5), p(5), pars(11), cont
double precision, intent(in) :: qd0(5), pd0(5)
double precision, intent(in) :: qd(5), pd(5)
double precision, intent(in) :: md0, qd0(5), pd0(5), contd0
double precision, intent(in) :: md, qd(5), pd(5), contd
double precision, intent(out) :: h
double precision, intent(out) :: hd
double precision, intent(out) :: hdd
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
from os import system
tap_dir = '/Users/ocots/Boulot/recherche/Logiciels/dev/hampath/hampath308/install/tapenade3.8'
#tap_dir = '/Users/ocots/Boulot/recherche/Logiciels/dev/hampath/hampath308/install/tapenade3.8'
tap_dir = '/home/ldellelc/Documents/Lavoro/Research/Matlab/hampath/install/tapenade3.8'
TAP = tap_dir+'/bin/tapenade'
tap_home = tap_dir
......@@ -8,5 +9,5 @@ build_dir='./'
#OPTT = '-tangent -fixinterface -inputlanguage fortran90 -outputlanguage fortran90 -nooptim activity -O ' + build_dir
OPTT = '-tangent -fixinterface -inputlanguage fortran90 -outputlanguage fortran90 -O ' + build_dir
system('TAPENADE_HOME=' + tap_home + ' ' + TAP + ' ' + OPTT + ' ' + '-tgtfuncname _d' + ' ' + '-root hfun.hfun_u -vars "q, p" -outvars "H" hfun.f90')
system('TAPENADE_HOME=' + tap_home + ' ' + TAP + ' ' + OPTT + ' ' + '-tgtfuncname _d' + ' ' + '-root hfun_d.hfun_u_d -vars "q, p" -outvars "hd" hfun_d.f90')
system('TAPENADE_HOME=' + tap_home + ' ' + TAP + ' ' + OPTT + ' ' + '-tgtfuncname _d' + ' ' + '-root hfun.hfun_u -vars "M, q, p, cont" -outvars "H" hfun.f90')
system('TAPENADE_HOME=' + tap_home + ' ' + TAP + ' ' + OPTT + ' ' + '-tgtfuncname _d' + ' ' + '-root hfun_d.hfun_u_d -vars "M, q, p, cont" -outvars "hd" hfun_d.f90')
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment