from sympy import *
t,mu=symbols('t,mu')
x=Function('x')(t)
y=Function('y')(t)
r=Function('r')(t)
theta=Function('theta')(t)
u_r=Matrix([cos(theta),sin(theta)])
init_printing()
u_r


u_theta=Matrix([-sin(theta),cos(theta)])
u_theta

"""
Explanation: <h2>Un poco de mecánica celeste</h2>

Los cálculos, como es ya habitual, se los dejaremos a SymPy, pero esta vez
usaremos <a href="http://ipython.org/">IPython</a> como interfaz de SymPy.
IPython permite  trabajar con notebooks, que son documentos que se visualizan con un navegador y permiten mezclar código Python (u otros lenguajes), con código html y LateX. La presente notebook se aloja en el siguiente <a href="https://github.com/fdmazzone/Ecuaciones_Diferenciales/blob/master/scripts/Mec%C3%A1nicaCeleste.ipynb">enlace</a>.
Los mismos cálculos se efectuaron con SAGE, y se pueden encontrar en el script 2cuerpos.sage dentro de la carpeta scripts en
el repositorio de <a href="https://github.com/fdmazzone/Ecuaciones_Diferenciales"> GitHub</a> que mantiene materiales de este curso. Ver el siguiente 
<a href="https://github.com/fdmazzone/Ecuaciones_Diferenciales}{https//github.com/fdmazzone/Ecuaciones_Diferenciales"> enlace </a>
Supondremos que el movimiento del planeta se realiza en un plano. Esta suposición no entraña una pérdida de generalidad, dado que es posible demostrar que es cierta. Esta demostración es sencilla, pero nos alejaría de nuestro objetivo, por tanto la omitiremos. 
Vamos a usar coordenadas polares $(r,\theta)$ y los versores $\overrightarrow{u}r:=(\cos\theta,\sin\theta)$ y $\overrightarrow{u}{\theta}:=(-\sin\theta, \cos\theta)$. Notar que
$\overrightarrow{u}r \perp \overrightarrow{u}{\theta}$ y por consiguiente $\mathcal{B}:={\overrightarrow{u}r , \overrightarrow{u}{\theta}}$ forma una base del espacio euclideano 2-dimensional. Usaremos este hecho
para representar distintos vectores como combinación lineal de vectores de la base.
End of explanation
"""


pos=r*u_r
a1,a2=symbols('a1,a2')
sol=solve(a1*u_r+a2*u_theta-pos.diff(t,2),[a1,a2])
sol

sol[a1]

sol[a2]

"""
Explanation: Vamos a representar  el vector aceleración $\overrightarrow{a}=\frac{d^2\overrightarrow{u}}{dt^2}$
End of explanation
"""


h=symbols('h')
ed=(sol[a1]).subs(theta.diff(t),h/r**2)
ed

"""
Explanation: El vector aceleración debe ser igual a la fuerza por unidad de masa $-\mu \overrightarrow{r}/r^3$. Notemos que esta fuerza es
<a href="http://es.wikipedia.org/wiki/Campo_central">central</a>, es decir tiene componente nula
respecto al vector $\overrightarrow{u}_{\theta}$. Por consiguiente se debe satisfacer que
$$\frac{1}{r}\frac{d}{dt} \left( r^2\dot{\theta} \right)=0\Longleftrightarrow \exists h\in\mathbb{R}: \boxed{ r^2\dot{\theta}=h}. $$
Hemos derivado la <a href=http://es.wikipedia.org/wiki/Leyes_de_Kepler>Segunda Ley de Kepler</a>: El radio vector barre áreas iguales en tiempos iguales.
<img src="https://raw.githubusercontent.com/fdmazzone/Ecuaciones_Diferenciales/master/imagenes/2leykepler.jpeg">
En la dirección radial $\overrightarrow{u}_r$ la componente de la fuerza es $-\mu/r^2$. Es decir se satisface la ecuación
$$
\ddot{r}-r\dot{\theta}^2=-\frac{\mu}{r^2}
$$
Notar que esta ecuación entraña dos incognitas $r$ y $\theta$, pero $\dot{\theta}$ puede ser remplazado por $h/r^2$ por la segunda Ley de Kepler.
Declaremos la variable $h$ que juega un rol importante y reemplacemos $\dot{\theta}$ en la ecuación
End of explanation
"""


z=Function('z')(theta)
r=1/z
ed2=r.diff(t,2)+mu/r**2-h**2/r**3
ed2

"""
Explanation: Conseguimos una ecuación no lineal de segundo orden para $r$. De los métodos que hemos visto, ninguno se aplica a esta ecuación.
El truco mágico consiste en considerar la nueva variable dependiente $z=1/r$ y la nueva variable independiente
$\theta$.
End of explanation
"""


theta2diff=(h*z**2).diff(t).subs(theta.diff(t),h*z**2)
ed3=ed2.subs(theta.diff(t,2),theta2diff)
ed3=ed3.subs(theta.diff(t),h*z**2)
ed3

(ed3/h**2/z**2).expand()

"""
Explanation: En la ecuación resultante, nuevamente aparece $\dot{\theta}$ y además ahora aparece $\ddot{\theta}$. Tenemos que reemplazar
$\dot{\theta}$ por $hz^2$ y $\ddot{\theta}$ por $\tfrac{d}{dt}hz^2$.
End of explanation
"""


