Indice | Previo: FormalismoHamiltoniano.ProblemasSeleccionados |
En este apéndice compilamos todas aquellos algoritmos y rutinas útiles usados en el libro y que pueden aprovecharse en la solución de una amplia gama de problemas en mecánica celeste.
Del vector de estado a los elementos orbitales clásicos.
def estado_a_elementos(mu,x):
#Posición y velocidad del sistema relativo
rvec=x[:3]
vvec=x[3:]
from numpy import cross
from numpy.linalg import norm
#Momento angular relativo específico
hvec=cross(rvec,vvec)
h=norm(hvec)
#Vector excentricidad
r=norm(rvec)
evec=cross(vvec,hvec)/mu-rvec/r
#Vector nodo ascendente
nvec=cross([0,0,1],hvec)
n=norm(nvec)
#Semilatus rectum y excentricidad
p=h**2/mu
e=norm(evec)
#Orientación
from numpy import dot,arccos,pi
i=arccos(hvec[2]/h)
Wp=arccos(nvec[0]/n)
W=Wp if nvec[1]>=0 else 2*pi-Wp
wp=arccos(dot(nvec,evec)/(e*n))
w=wp if evec[2]>=0 else 2*pi-wp
fp=arccos(dot(rvec,evec)/(r*e))
f=fp if dot(rvec,vvec)>0 else 2*pi-fp
return p,e,i,W,w,f
De los elementos orbitales clásicos al vector de estado:
def elementos_a_estado(mu,elementos):
#Extrae elementos
p,e,i,W,w,f=elementos
#Calcula momento angular relativo específico
from numpy import sqrt
h=sqrt(mu*p)
#Calcula r
from numpy import cos
r=p/(1+e*cos(f))
#Posición
from numpy import cos,sin
x=r*(cos(W)*cos(w+f)-cos(i)*sin(W)*sin(w+f))
y=r*(sin(W)*cos(w+f)+cos(i)*cos(W)*sin(w+f))
z=r*sin(i)*sin(w+f)
#Velocidad
muh=mu/h
vx=muh*(-cos(W)*sin(w+f)-cos(i)*sin(W)*cos(w+f))\
-muh*e*(cos(W)*sin(w)+cos(w)*cos(i)*sin(W))
vy=muh*(-sin(W)*sin(w+f)+cos(i)*cos(W)*cos(w+f))\
+muh*e*(-sin(W)*sin(w)+cos(w)*cos(i)*cos(W))
vz=muh*(sin(i)*cos(w+f)+e*cos(w)*sin(i))
from numpy import array
return array([x,y,z,vx,vy,vz])
Método de Newton general:
def metodo_newton(f,x0=1,delta=1e-5,args=()):
#Valor inicial de la anomalía excéntrica
xn=x0
#Valor inicial del error relativo
Dn=1
#Contador de iteraciones
ni=0
while Dn>delta:
#Inicializa el valor de En
x=xn
#Nuevo valor (regla de iteración)
xn=x-f(x,*args)[0]/f(x,*args)[1]
#Valor medio
xmed=(x+xn)/2
#Criterio de convergencia
en=xn-x
Dn=abs(en/xmed)
ni+=1
return xmed,Dn,ni
Método de Laguerre-Conway:
def metodo_laguerre(f,x0=1,delta=1e-5,args=(),eta=5):
#Varifica que el valor inicial sea apropiado
disc=-1
mi=0
#Valor inicial de la anomalía excéntrica
xn=x0
#Valor inicial del error relativo
Dn=1
#Contador de iteraciones
ni=0
while Dn>delta:
#Inicializa el valor de En
x=xn
disc=-1
mi=0
while disc<0:
mi+=1
#Valor de la función y sus derivadas
y,yp,ypp=f(x,*args)
#Discriminante
disc=(eta-1)**2*yp**2-eta*(eta-1)*y*ypp
eta=eta-1 if disc<0 else eta
#Raiz del discriminante
from numpy import sqrt
raiz_disc=sqrt(disc)
#Signo en el denominador
sgn=+1 if abs(yp+raiz_disc)>abs(yp-raiz_disc) else -1
#Valor de en
en=eta*y/(yp+sgn*raiz_disc)
#Nuevo valor (regla de iteración)
xn=x-en
#Valor medio
xmed=(x+xn)/2
#Criterio de convergencia
en=xn-x
Dn=abs(en/xmed)
ni+=1
return xmed,Dn,ni+mi-1
La fuente original de este algoritmo esta disponible en línea.
def kepler_semianalitico(M,e):
from math import sin,cos,pi
#Casos extremos
if M==0 or M==2*pi or e==1:return M,0,0
Minp=M
Ecorr=0;Esgn=1.0
if M>pi:
M=2*pi-M
Ecorr=2*pi
Esgn=-1.0
#Circunferencia
if e==0:return Ecorr+Esgn*M,0,0
a=(1-e)*3/(4*e+0.5);
b=-M/(4*e+0.5);
y=(b**2/4 +a**3/27)**0.5;
x=(-0.5*b+y)**(1./3)-(0.5*b+y)**(1./3);
w=x-0.078*x**5/(1 + e);
E=M+e*(3*w-4*w**3);
#Corrección por Newton
sE=sin(E)
cE=cos(E)
f=(E-e*sE-M);
fd=1-e*cE;
f2d=e*sE;
f3d=-e*cE;
f4d=e*sE;
E=E-f/fd*(1+f*f2d/(2*fd*fd)+\
f*f*(3*f2d*f2d-fd*f3d)/(6*fd**4)+\
(10*fd*f2d*f3d-15*f2d**3-fd**2*f4d)*\
f**3/(24*fd**6))
#Corrección por Newton
f=(E-e*sE-M);
fd=1-e*cE;
f2d=e*sE;
f3d=-e*cE;
f4d=e*sE;
E=E-f/fd*(1+f*f2d/(2*fd*fd)+\
f*f*(3*f2d*f2d-fd*f3d)/(6*fd**4)+\
(10*fd*f2d*f3d-15*f2d**3-fd**2*f4d)*\
f**3/(24*fd**6))
E=Ecorr+Esgn*E
#Error relativo
Mnum=E-e*sin(E)
Dn=abs(Mnum-Minp)/Minp
return E,Dn,1
Si continuamos indefinidamente el proceso anterior es posible mostrar que la anomalía excéntrica puede calcularse como una serie de Fourier de la forma (Plummer, 1918):
\begin{equation} \label{eq:kepler_fourier} E =M + \sum_{n=1}^{\infty}\frac{e^n}{2^{n-1}}\sum_{k=0}^{\lfloor n/2\rfloor} a_{nk}\sin [(n-2k)M] \end{equation}donde:
$$ a_{nk} = (-1)^k\frac{(n-2k)^{n-1}}{(n-k)!k!} $$Puede probarse que esta serie converge uniformemente para todos los valores de $e<0.662743419349181$.$^1$
def kepler_eserie(M,e,delta=0,orden=1):
from math import sin,factorial,floor
nfac=1
En=M
Dn=1
n=0
condicion=Dn>delta if delta>0 else n<=orden
while condicion:
n+=1
E=En
prefactor=e**n/2**(n-1)
kmax=int(floor(n/2))
sgn=-1
#Los factoriales se calculan así para mayor eficiencia
nfac=nfac*n if n>0 else 1
kfac=1
nkfac=1
termino=0
for k in range(kmax+1):
sgn*=-1
kfac=kfac*k if k>0 else 1
nkfac=nkfac/(n-k+1) if k>0 else nfac
ank=sgn/(kfac*nkfac)*(n-2*k)**(n-1)
termino+=ank*sin((n-2*k)*M)
dE=prefactor*termino
En+=dE
Dn=abs(dE/En)
#La condicion depende de si se pasa o no la tolerancia
condicion=Dn>delta if delta>0 else n<orden
return En,Dn,n
Esta rutina implementa la serie en la Ec. kepler_bessel:
def kepler_bessel(M,e,delta):
from math import sin
from scipy.special import jv
Dn=1
n=1
En=M
while Dn>delta:
E=En
dE=(2./n)*jv(n,n*e)*sin(n*M)
En+=dE
Emed=(E+En)/2
Dn=abs(dE/Emed)
n+=1
return En,Dn,n
Series de Stumpff. Devuelve el valor de la series y de su primera y segunda derivada.
def serie_stumpff(t,k,N=15):
from math import factorial
sk=lambda n:t/((2*n+k+1)*(2*n+k+2))*(1-sk(n+1))\
if n<N else 0
return (1-sk(0))/factorial(k)
NOTAS AL PIE: