Построение орбиты в Python, но когда x <0, траектория внезапно становится линейной

Question

Я строю в простой 2D-плоскости путь, пройденный телом, движущимся мимо другого гравитационно-привлекательного тела.

В каждом цикле время увеличивается на секунду, и вычисляется новое положение тела, ираспечатанный.Затем я вставляю результаты в электронную таблицу и строю график.

Все выглядит хорошо, пока компонент тела x не станет отрицательным - и затем траектория станет линейной и улетит в верхний левый угол.

Это все для моего собственного сольного развлечения, я не студент.Поэтому, немного почесав голову, я, наконец, вздохнул, попросив кого-нибудь о помощи.Я, наверное, пропустил что-то очевидное.Я подозреваю, что моей тригонометрии чего-то не хватает.

Я использую Python 2.7.10

import sys
import os
import math

mu = 4.0*(10**14)
massAst = 1
earthRadius = 6371000.
alt = 100000.
r = earthRadius+ alt
rTheta = 270.
rAngtoX = math.radians(rTheta)
tInc = 1 ## increment time by 1 seconds - one update of pos&velocity per second of travel
calcTime = 1100 ## simulation runtime (86400 seconds = 1 day) 10 mins
t = 1 ## integral of time t to use in the calcs in the loop.
printbell = 120 ## print results now
printclock = 0
hourCount = 0

## Initialise velocity vectors for Asteroid:
uAstX = 1500.
uAstY = 0.
vAstX = 0.
vAstY = 0.
## Displacement
dAstX = r*math.cos(rAngtoX)
dAstY = r*math.sin(rAngtoX)

for i in range(0, calcTime):
    acc = -1*(mu/r**2)

    accX = acc*math.cos(rAngtoX)
    accY = acc*math.sin(rAngtoX)

    vAstX = uAstX + accX*t ## new value for velocity in X direction
    vAstY = uAstY + accY*t ## and in Y

    deltaDAstX = uAstX*t + 0.5*accX*(t**2) ## change in position over this time interval
    deltaDAstY = uAstY*t + 0.5*accY*(t**2)

    dAstX = dAstX + deltaDAstX
    dAstY = dAstY + deltaDAstY

    uAstX = vAstX
    uAstY = vAstY
    ## Now calculate new angle and range
    ## tan(theta) = dAstY/dAstX, so:
    rAngtoX = math.atan(dAstY/dAstX) ##+(2*3.141592654)
    ##print 'theta:', math.degrees(rAngtoX)
    r = dAstY/math.sin(rAngtoX)
    ## if i == print
    print dAstX, '  ', dAstY

Answer 1

Я не знаком с математикой, но я взял ваш код, изменил его для запуска на моей машине, подготовил данные с использованием библиотеки seaborn и придумал следующее:

Это код, который я использовал:

import math
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt


def calc(calcTime):
    mu = 4.0*(10**14)
    earthRadius = 6371000.0
    alt = 100000.0
    r = earthRadius+ alt
    rTheta = 270.0
    rAngtoX = math.radians(rTheta)
    t = 1               # integral of time t to use in the calcs in the loop.

    # Initialise velocity vectors for Asteroid:
    uAstX = 1500.0
    uAstY = 0.0
    # Displacement
    dAstX = r*math.cos(rAngtoX)
    dAstY = r*math.sin(rAngtoX)

    point_list = []
    for i in range(0, calcTime):
        acc = -1*(mu/r**2)

        accX = acc*math.cos(rAngtoX)
        accY = acc*math.sin(rAngtoX)

        vAstX = uAstX + accX*t                  # new value for velocity in X direction
        vAstY = uAstY + accY*t                  # and in Y

        deltaDAstX = uAstX*t + 0.5*accX*(t**2)      # change in pos over time interval
        deltaDAstY = uAstY*t + 0.5*accY*(t**2)

        dAstX = dAstX + deltaDAstX
        dAstY = dAstY + deltaDAstY

        uAstX = vAstX
        uAstY = vAstY
        # Now calculate new angle and range
        # tan(theta) = dAstY/dAstX, so:
        rAngtoX = math.atan(dAstY/dAstX)            #+(2*3.141592654)
        # print 'theta:', math.degrees(rAngtoX)
        r = dAstY/math.sin(rAngtoX)
        # if i == print
        if i % 5 == 0:
            print('{:05d} | {:15.2f} | {:15.2f}'.format(i, dAstX, dAstY))
            point_list.append((i, dAstX, dAstY))

    df = pd.DataFrame(data=point_list, columns=['i', 'x', 'y'])
    return df


if __name__ == '__main__':
    df = calc(950)
    sns.scatterplot(data=df, x='x', y='y')
    plt.show()

---

Мой анализ: промежутки между точками с каждым разом увеличиваются (примечание: я только заговариваю каждый5-й пункт, чтобы сделать изображение более читабельным, на мой взгляд).С точки зрения физики это указывало бы на то, что объект набирает скорость.

Возможно ли, что ваши расчеты верны и что объект покидает орбиту, потому что он набрал достаточно скорости (или энергии), чтобы уйтигравитационное поле центральной массы (иначе Земли)?

Как я уже сказал, я не знаком с конкретной математикой, но для меня имеет смысл, что объект может вырваться из орбиты сскорость набирает в пол оборота.

Answer 2

Когда dAstX приближается к нулю, dAstY/dAstX приближается к делению на ноль ... Это вызовет всевозможные проблемы (по крайней мере, проблемы округления).

Я бы порекомендовалразделяя x / y компоненты расстояния / скорости / ускорения отдельно.Конечно, расстояние между объектами важно, но его можно рассчитать, используя r=sqrt(dAstX**2 + dAstY**2).