Este código fue desarrollado y utilizado para la materia Evolución Tidal de Sistemas Planetarios, de la FaMAF (UNC). 🇦🇷
Si principal función es la evolución de un sistema orbital de 3 cuerpos, incluyendo efectos tidales, gravitatorios y relativistas; utilizando un sistema de referencia de Jacobi.
Al comienzo del archivo main.F90 se encuentra un bloque para intorucir los parámetros de los 3 cuerpos, y de la integración.
Los parametros iniciales de los cuerpos utilizados son:
- m: Masa m. [M⊙]
- a: Semieje mayor a. [UA]
- e: Eccentricidad e.
- s: Spin Ωs. [Rad Día⁻¹]
- o: Oblicuidad ε. [Rad]
- vp: Longitud del pericentro . [Rad]
- r: Radio R. [UA]
- alp: Alpha de momento de inercia α (0.4 ⇒ esfera).
- Q: Parámetro tidal (Tidal Quality Factor).
Los parametros iniciales de la integración son:
- t0: Tiempo inicial t0. [Día]
- dt_min: Paso de tiempo para integradores de paso fijo, o mínimo dtmin. [Día]
- tf: Tiempo final tf. [Día]
- beta: Tasa de aprendizaje, en caso que se utilize un integrador de paso adaptativo β.
- e_tol: Error absoluto máximo, en caso que se utilize un integrador de paso adaptativo ϵtol (≡ |yreal - ypred|).
- n_points: Cantiad aproximada de datos de salida.
- logsp: Si la salida es log- (True) o equi- (False) espaciada.
- filename: Nombre del archivo de salida.
Se puede modificar el integrador utilizado, como también las fuerzas involucradas.
Se debe utilar uno de los 4 métodos disponibles: integ_caller, rk_half_step_caller, embedded_caller, o BStoer_caller.
Esto se realiza en la líneas 149-153 del archivo main.F90.
!!! Execute an integration method (uncomment/edit one of theese)
! call integ_caller (t, y, dt_adap, dydtidall, Runge_Kutta4, dt, ynew)
! call rk_half_step_caller (t, y, dt_adap, dydtidall, Runge_Kutta5, 5, e_tol, beta, dt_min, dt, ynew)
call embedded_caller (t, y, dt_adap, dydtidall, Dormand_Prince8_7, e_tol, beta, dt_min, dt, ynew)
! call BStoer_caller (t, y, dt_adap, dydtidall, e_tol, dt_min, dt, ynew)
-
integ_caller: Se puede modificar el integrador (
Runge_Kutta4
) por cualquier otro que no sea embebido, entre los listados en el archivo integrators.F90. El paso dtadap será constante, el igual al valor mínimo dtmin. -
rk_half_step_caller: Igual a integ_caller, pero al cambiar el integrador, también se debe introducir su orden (ej.
Runge_Kutta5
-> O(5)). En este caso el integrador intentará utilizar un paso de tiempo adaptativo dtadap adecuado, según la toleranca de error ϵtol introducida (ver integrators.F90). En este caso el error calculado será ϵcalc ≡ |1 - (yaux/ypred)|/(2ord - 1). -
embedded_caller: Similar a rk_half_step_caller (debido a que también calcula un paso de tiempo óptimo), pero solo se pueden introducir integradores embebidos (ver integrators.F90). En este caso el error calculado será ϵcalc ≡ |1 - (yaux/ypred)|.
-
BStoer_caller:, Se utiliza el integrador Bulirsch_Stoer (ver bstoer.F90).
Esto se realiza en la líneas 358-362 del archivo tidall.F90.
call dydtidal (y0, y1, 1, y01t, y10t)
call dydtidal (y0, y2, 2, y02t, y20t)
call dydtgrav (y1, 1, y2, 2, y12g, y21g)
call dydtrela (y1, 1, y1r)
call dydtrela (y2, 2, y2r)
En caso de querer despreciar (o simplemente no calcular) alguna fuerza en especial (tidal, gravitatoria o efectos relativistas), se debe comentar la línea correspondiente.
Compilar y ejecutar:
gfortran
Plotear:
python
numpy
matplotlib
pandas
Recompilar y ejecutar (en caso de realizar modificaciones al código que cambien las dependencias):
pip
python
- fortdepend (paquete de python utilzado para generar nuevo archivo de dependencias)
Si se usa el código tal y como está, solo es necesario compilarlo y ejecutarlo. En caso de realizarle modificaciones que alteren las dependencias, es necesario instalar fortdepend.
$ make
$ ./tidal
En caso de que no exista el archivo tidal.dep, es necesario tener instalado fortdepend para generarlo.
$ make install
$ make clean
En este caso, también se borra el archivo de dependencias tidal.dep, por lo que luego será necesario regenerarlo.
Se adjunta un pequeño código plot.py en python
para realizar algunos plots resultantes.
$ python plot.py
Distribucído bajo la Licencia MIT. (Archivo LICENSE)
Emmanuel Gianuzzi - [email protected]