Ejercicio: Simulación de mecánica orbital 2D (python)

Question

Ejercicio: Simulación de mecánica orbital 2D (python)

orbita
Espacio
gravedad
n-cuerpo-simulaciones

espacio de Estados

Solo un pequeño descargo de responsabilidad de antemano: nunca he estudiado astronomía ni ninguna ciencia exacta (ni siquiera TI), por lo que estoy tratando de llenar este vacío mediante la autoeducación. La astronomía es una de las áreas que más me ha llamado la atención y mi idea de la autoeducación es directa con un enfoque aplicado. Entonces, directo al grano: este es un modelo de simulación orbital en el que estoy trabajando casualmente cuando tengo tiempo/animo. Mi objetivo principal es crear un sistema solar completo en movimiento y la capacidad de planificar lanzamientos de naves espaciales a otros planetas.

¡Todos son libres de retomar este proyecto en cualquier momento y divertirse experimentando!

¡¡¡actualizar!!! (10 de noviembre)

la velocidad ahora es deltaV propio y al dar movimiento adicional ahora calcula el vector de suma de la velocidad
puede colocar tantos objetos estáticos como desee, en cada objeto unitario en movimiento, comprueba los vectores de gravedad de todas las fuentes (y comprueba las colisiones)
mejorado mucho el rendimiento de los cálculos
una solución para tener en cuenta la modificación interactiva en matplotlib. Parece que esta es la opción predeterminada solo para ipython. Python3 regular requiere esa declaración explícitamente.

Básicamente, ahora es posible "lanzar" una nave espacial desde la superficie de la Tierra y trazar una misión a la Luna haciendo correcciones del vector deltaV a través de giveMotion(). Lo siguiente en la línea es intentar implementar una variable de tiempo global para permitir el movimiento simultáneo, por ejemplo, la Luna orbita la Tierra mientras la nave espacial prueba una maniobra de asistencia por gravedad.

¡Los comentarios y sugerencias para mejorar son siempre bienvenidos!

Hecho en Python3 con la biblioteca matplotlib

import matplotlib.pyplot as plt
import math
plt.ion()

G = 6.673e-11  # gravity constant
gridArea = [0, 200, 0, 200]  # margins of the coordinate grid
gridScale = 1000000  # 1 unit of grid equals 1000000m or 1000km

plt.clf()  # clear plot area
plt.axis(gridArea)  # create new coordinate grid
plt.grid(b="on")  # place grid

class Object:
    _instances = []
    def __init__(self, name, position, radius, mass):
        self.name = name
        self.position = position
        self.radius = radius  # in grid values
        self.mass = mass
        self.placeObject()
        self.velocity = 0
        Object._instances.append(self)

    def placeObject(self):
        drawObject = plt.Circle(self.position, radius=self.radius, fill=False, color="black")
        plt.gca().add_patch(drawObject)
        plt.show()

    def giveMotion(self, deltaV, motionDirection, time):
        if self.velocity != 0:
            x_comp = math.sin(math.radians(self.motionDirection))*self.velocity
            y_comp = math.cos(math.radians(self.motionDirection))*self.velocity
            x_comp += math.sin(math.radians(motionDirection))*deltaV
            y_comp += math.cos(math.radians(motionDirection))*deltaV
            self.velocity = math.sqrt((x_comp**2)+(y_comp**2))

            if x_comp > 0 and y_comp > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity))  # update motion direction
            elif x_comp > 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 90
            elif x_comp < 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 270

        else:
            self.velocity = self.velocity + deltaV  # in m/s
            self.motionDirection = motionDirection  # degrees
        self.time = time  # in seconds
        self.vectorUpdate()

    def vectorUpdate(self):
        self.placeObject()
        data = []

        for t in range(self.time):
            motionForce = self.mass * self.velocity  # F = m * v
            x_net = 0
            y_net = 0
            for x in [y for y in Object._instances if y is not self]:
                distance = math.sqrt(((self.position[0]-x.position[0])**2) +
                             (self.position[1]-x.position[1])**2)
                gravityForce = G*(self.mass * x.mass)/((distance*gridScale)**2)

                x_pos = self.position[0] - x.position[0]
                y_pos = self.position[1] - x.position[1]

                if x_pos <= 0 and y_pos > 0:  # calculate degrees depending on the coordinate quadrant
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+90

                elif x_pos > 0 and y_pos >= 0:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))+180

                elif x_pos >= 0 and y_pos < 0:
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+270

                else:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))

                x_gF = gravityForce * math.sin(math.radians(gravityDirection))  # x component of vector
                y_gF = gravityForce * math.cos(math.radians(gravityDirection))  # y component of vector

                x_net += x_gF
                y_net += y_gF

            x_mF = motionForce * math.sin(math.radians(self.motionDirection))
            y_mF = motionForce * math.cos(math.radians(self.motionDirection))
            x_net += x_mF
            y_net += y_mF
            netForce = math.sqrt((x_net**2)+(y_net**2))

            if x_net > 0 and y_net > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce))  # update motion direction
            elif x_net > 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 90
            elif x_net < 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 270

            self.velocity = netForce/self.mass  # update velocity
            traveled = self.velocity/gridScale  # grid distance traveled per 1 sec
            self.position = (self.position[0] + math.sin(math.radians(self.motionDirection))*traveled,
                             self.position[1] + math.cos(math.radians(self.motionDirection))*traveled)  # update pos
            data.append([self.position[0], self.position[1]])

            collision = 0
            for x in [y for y in Object._instances if y is not self]:
                if (self.position[0] - x.position[0])**2 + (self.position[1] - x.position[1])**2 <= x.radius**2:
                    collision = 1
                    break
            if collision != 0:
                print("Collision!")
                break

        plt.plot([x[0] for x in data], [x[1] for x in data])

Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(100.0, 100.0), radius=1.737, mass = 7.347e22)  # position not to real scale
Craft = Object(name="SpaceCraft", position=(49.0, 40.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)
Craft.giveMotion(deltaV=2000.0, motionDirection=90, time=60000)
plt.show(block=True)

Cómo funciona

Todo se reduce a dos cosas:

Creación de objetos como Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)con parámetros de posición en la cuadrícula (1 unidad de cuadrícula es 1000 km por defecto, pero esto también se puede cambiar), radio en unidades de cuadrícula y masa en kg.
Dando al objeto algo deltaV como Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)obviamente requiere Craft = Object(...)ser creado en primer lugar como se mencionó en el punto anterior. Los parámetros aquí están deltaVen m/s (tenga en cuenta que por ahora la aceleración es instantánea), motionDirectiones la dirección de deltaV en grados (desde la posición actual, imagine un círculo de 360 grados alrededor del objeto, por lo que la dirección es un punto en ese círculo) y finalmente el parámetro timees cuántos segundos después de que la trayectoria de empuje deltaV del objeto sea monitoreada. Las subsiguientes giveMotion()comienzan desde la última posición de la anterior giveMotion().

Preguntas:

¿Es este un algoritmo válido para calcular órbitas?
¿Cuáles son las mejoras obvias a realizar?
He estado considerando la variable "timeScale" que optimizará los cálculos, ya que podría no ser necesario volver a calcular los vectores y las posiciones para cada segundo. ¿Alguna idea sobre cómo debería implementarse o, en general, es una buena idea? (pérdida de precisión vs rendimiento mejorado)

Básicamente, mi objetivo es iniciar una discusión sobre el tema y ver a dónde conduce. Y, si es posible, aprenda (o mejor aún, enseñe) algo nuevo e interesante.

¡Siéntete libre de experimentar!

Intenta usar:

Earth = Object(name="Earth", position=(50.0, 100.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(434.0, 100.0), radius=1.737, mass = 7.347e22)
Craft = Object(name="SpaceCraft", position=(43.0, 100.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=10575.0, motionDirection=180, time=322000)
Craft.giveMotion(deltaV=400.0, motionDirection=180, time=50000)

Con dos encendidos, uno progrado en la órbita terrestre y otro retrógrado en la órbita lunar, logré una órbita lunar estable. ¿Están cerca de los valores teóricos esperados?

Ejercicio sugerido: Pruébelo en 3 encendidos: órbita terrestre estable desde la superficie de la Tierra, encendido progresivo para llegar a la Luna, encendido retrógrado para estabilizar la órbita alrededor de la Luna. Luego intente minimizar deltaV.

Nota: Planeo actualizar el código con extensos comentarios para aquellos que no están familiarizados con la sintaxis de python3.

usuario2449

¡Una muy buena idea para la autoeducación! ¿Sería posible resumir sus fórmulas para aquellos de nosotros que no estamos familiarizados con la sintaxis de Python?

espacio de Estados

Seguro supongo. Haré comentarios más extensos en el código para aquellos que quieran retomarlo y resumir la lógica general en la pregunta misma.

usuario21

Lo mejor de mi cabeza: considere usar un vector para la velocidad en lugar de tratar la velocidad y la dirección de manera diferente. Cuando dices "F = m * v", ¿quieres decir "F = m*a"? ¿Estás asumiendo que la Tierra no se mueve porque es mucho más pesada que un asteroide? Considere consultar github.com/barrycarter/bcapps/blob/master/bc-grav-sim.pl

espacio de Estados

Puedes darle movimiento a cualquier objeto, incluida la Tierra. Para fines de prueba, incluí solo objeto -> Relación con la tierra en el bucle principal. Se puede convertir fácilmente que cada objeto se relacione con todos los demás objetos que se crean. Y cada objeto puede tener su propio vector de movimiento. Razón por la que no lo hice: cálculos muy lentos incluso para 1 objeto. Espero que escalar las unidades de tiempo ayude mucho, pero todavía no estoy seguro de cómo hacerlo bien.

usuario21

ESTÁ BIEN. Un pensamiento: haga la simulación para dos objetos reales (p. ej., Tierra/Luna o Tierra/Sol) y compare sus resultados con ssd.jpl.nasa.gov/?horizons para ver la precisión. No será perfecto debido a las perturbaciones de otras fuentes, pero le dará una idea de la precisión.

espacio de Estados

Genial gracias. Solo ese enlace valió la pena para iniciar la discusión. Entonces, primer paso: intentaré volver a escribir el código para que los objetos verifiquen todos los demás objetos en busca del vector de fuerza de gravedad. Por cierto, el tema más apremiante de la escala de tiempo todavía está aquí. Y sería interesante comprobar cómo mejoran o empeoran los resultados cambiando el parámetro de escala de tiempo.

espacio de Estados

Ah, y F=m*v es realmente eso. Se supone que al dar movimiento al objeto en primer lugar, se detiene por completo y por eso la velocidad es igual a la aceleración. Como la velocidad instantánea desde 0 -> cualquier valor, y continúa indefinidamente con esa velocidad si ninguna otra fuerza actúa sobre el objeto.

Respuestas (1)

Ejercicio: Simulación de mecánica orbital 2D (python)

¡Una muy buena idea para la autoeducación! ¿Sería posible resumir sus fórmulas para aquellos de nosotros que no estamos familiarizados con la sintaxis de Python?
Seguro supongo. Haré comentarios más extensos en el código para aquellos que quieran retomarlo y resumir la lógica general en la pregunta misma.
Lo mejor de mi cabeza: considere usar un vector para la velocidad en lugar de tratar la velocidad y la dirección de manera diferente. Cuando dices "F = m * v", ¿quieres decir "F = m*a"? ¿Estás asumiendo que la Tierra no se mueve porque es mucho más pesada que un asteroide? Considere consultar github.com/barrycarter/bcapps/blob/master/bc-grav-sim.pl
Puedes darle movimiento a cualquier objeto, incluida la Tierra. Para fines de prueba, incluí solo objeto -> Relación con la tierra en el bucle principal. Se puede convertir fácilmente que cada objeto se relacione con todos los demás objetos que se crean. Y cada objeto puede tener su propio vector de movimiento. Razón por la que no lo hice: cálculos muy lentos incluso para 1 objeto. Espero que escalar las unidades de tiempo ayude mucho, pero todavía no estoy seguro de cómo hacerlo bien.
ESTÁ BIEN. Un pensamiento: haga la simulación para dos objetos reales (p. ej., Tierra/Luna o Tierra/Sol) y compare sus resultados con ssd.jpl.nasa.gov/?horizons para ver la precisión. No será perfecto debido a las perturbaciones de otras fuentes, pero le dará una idea de la precisión.
Genial gracias. Solo ese enlace valió la pena para iniciar la discusión. Entonces, primer paso: intentaré volver a escribir el código para que los objetos verifiquen todos los demás objetos en busca del vector de fuerza de gravedad. Por cierto, el tema más apremiante de la escala de tiempo todavía está aquí. Y sería interesante comprobar cómo mejoran o empeoran los resultados cambiando el parámetro de escala de tiempo.
Ah, y F=m*v es realmente eso. Se supone que al dar movimiento al objeto en primer lugar, se detiene por completo y por eso la velocidad es igual a la aceleración. Como la velocidad instantánea desde 0 -> cualquier valor, y continúa indefinidamente con esa velocidad si ninguna otra fuerza actúa sobre el objeto.

david ketcheson · Answer 1

Que los dos cuerpos involucrados tengan masas $m_1, m_2$ . Comience con la segunda ley de Newton

F = metro a

$F = ma$ dónde

a

$a$ es aceleración. La fuerza gravitacional sobre el cuerpo 2 del cuerpo 1 está dada por

F_{21} = \frac{GRAMO {metro}_{1} {metro}_{2}}{| r_{21} |^{3}} r_{21}

$F_{21} = \frac{G m_1 m_2}{|r_{21}|^3}r_{21}$

dónde $r_{21}$ es el vector de posición relativa para los dos cuerpos en cuestión. La fuerza sobre el cuerpo 1 del cuerpo dos es, por supuesto $F_{12}=-F_{21}$ desde $r_{12}=-r_{21}$ . Si denotamos las posiciones por $(x_1,y_1)$ y $(x_2,y_2)$ , después

r_{21} = (\begin{matrix} X_{1} - X_{2} \\ y_{1} - y_{2} \end{matrix}) .

$r_{21} = \begin{pmatrix} x_1-x_2 \\ y_1-y_2 \end{pmatrix}.$

y

| r | = \sqrt{(X_{1} - X_{2})^{2} + (y_{1} - y_{2})^{2}} .

$|r| = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2}.$ Entonces la ley de Newton

a = F / m

$a=F/m$ se convierte

\begin{aligned} X_{1}^{″} (t) & = \frac{GRAMO {metro}_{2} (X_{2} - X_{1})}{| r |^{3}} \\ y_{1}^{″} (t) & = \frac{GRAMO {metro}_{2} (y_{2} - y_{1})}{| r |^{3}} \\ X_{2}^{″} (t) & = \frac{GRAMO {metro}_{1} (X_{1} - X_{2})}{| r |^{3}} \\ y_{2}^{″} (t) & = \frac{GRAMO {metro}_{1} (y_{1} - y_{2})}{| r |^{3}} . \end{aligned}

$\begin{align} x_1''(t) & = \frac{G m_2 (x_2-x_1)}{|r|^3} \\ y_1''(t) & = \frac{G m_2 (y_2-y_1)}{|r|^3} \\ x_2''(t) & = \frac{G m_1 (x_1-x_2)}{|r|^3} \\ y_2''(t) & = \frac{G m_1 (y_1-y_2)}{|r|^3}. \end{align}$

Junto con las posiciones y velocidades iniciales, este sistema de ecuaciones diferenciales ordinarias (EDO) comprende un problema de valor inicial. El enfoque habitual es escribir esto como un sistema de primer orden de 8 ecuaciones y aplicar un método de Runge-Kutta o de varios pasos para resolverlas.

Si aplica algo simple como Euler hacia adelante o Euler hacia atrás, verá que la Tierra gira en espiral hacia el infinito o hacia el sol, respectivamente, pero eso es un efecto de los errores numéricos. Si usa un método más preciso, como el método clásico de Runge-Kutta de cuarto orden, encontrará que se mantiene cerca de una órbita real por un tiempo, pero eventualmente se mueve en espiral hacia adentro (siempre que el paso de tiempo sea lo suficientemente pequeño) o va al infinito (si el paso de tiempo es demasiado grande). El enfoque correcto es utilizar un método simpléctico, que mantendrá a la Tierra en la órbita correcta, aunque su fase seguirá estando desfasada debido a errores numéricos.

Para el problema de 2 cuerpos, es posible derivar un sistema más simple al basar su sistema de coordenadas alrededor del centro de masa. Pero creo que la formulación anterior es más clara si esto es nuevo para usted.

Todavía digiriendo. Demasiadas palabras desconocidas para mí, pero de alguna manera siento que en algún momento llegaré. Por ahora, mi propio algoritmo es lo suficientemente bueno para que las cosas simplemente funcionen. Pero cuando conecte el movimiento simultáneo, me veré obligado a consultar la literatura y leer sobre los algoritmos adecuados. Dado que las limitaciones del hardware moderno son mucho más flexibles, puedo darme el lujo de perder el tiempo con ecuaciones simples. Sin embargo, temo que no por mucho tiempo.
De hecho, los métodos simplécticos son, con mucho, los más precisos, pero creo que es difícil para alguien sin experiencia en ciencias implementarlos. En su lugar, puede usar el método de Euler muy simple junto con la corrección de Feynman. No creo que necesites algo más complejo que eso para fines de autoeducación.

Ejercicio: Simulación de mecánica orbital 2D (python)

espacio de Estados

usuario2449

espacio de Estados

usuario21

espacio de Estados

usuario21

espacio de Estados

espacio de Estados

Respuestas (1)

david ketcheson

espacio de Estados

espacio de Estados

Chris Papageorgakis

¿Es posible tener un sistema estable de 3 cuerpos que orbite en un círculo perfecto?

Si hubiera una sonda orbitando en el borde de la Esfera de influencia de la Tierra, ¿a qué velocidad orbitaría?

¿Puedes explicar el patrón de los tamaños de las esferas de Hill de los objetos del sistema solar?

Al calcular los seis parámetros orbitales keplerianos, ¿por qué necesitamos tanto la excentricidad como el semieje mayor? ¿No te dice uno el otro?

¿Podría un dron escapar de la atmósfera de Marte?

¿Cuál es la diferencia entre la esfera de influencia y la esfera de la colina?

¿Cuál es la trayectoria de escape óptima cerca de un agujero negro?

¿Está la superficie de la Tierra "en órbita"?

¿Cómo puedo calcular los períodos orbitales en un sistema estelar binario?

¿Se podría pasar una luna de la órbita de un planeta a la de otro?