Проблема обновления угловой скорости и углов в моделировании орбиты двойной звезды (Python, Matplotlib) ⇐ Python

[Anonymous]

Я создаю симуляцию двойной звездной системы на Python, используя matplotlib. Я хочу, чтобы он менялся в зависимости от действий пользователей с помощью предоставленных мной ползунков. Когда я впервые создал это без каких-либо реальных уравнений для двойных звездных систем (w = sqrt(G*M/r) и тета += w * dt), это сработало отлично. Когда массы каждой звезды изменились/расстояние между ними, индивидуальные расстояния до звезд тоже изменились. Но когда я пытаюсь изменить скорость, общее расстояние между звездами не остается постоянным. Они всегда должны двигаться навстречу друг другу. Я не могу найти точную проблему и чувствую, что это может быть комбинация нескольких. Я также стараюсь избегать использования эллиптических орбит.
Мой код:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider
import scipy
from math import sqrt

dt = 0.1
numsteps = 10000
pi = scipy.constants.pi
G = 4.30091e-3 # AU^3 * M_sun^-1 * yr^-2
wA =3.0
wB =3.0

thetaA = 0
thetaB = thetaA+pi #to put the other store on the opposite end of starA

#initialise variables
r = 5
mA = 50 #mass in solar mass
mB = 50
M = mA+mB

x_valA,y_valA = [],[]
x_valB,y_valB = [],[]

# Create the animation, fig represents the object/canvas and ax means it is the area being plotted on
fig, ax = plt.subplots()
# the graph, sets limits for axis
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')

ax.set_facecolor("black")
fig.patch.set_facecolor("black")

# Create the stars and COM plot
starA, = ax.plot([], [], 'o', color='blue', markersize=10, label='Star A')
starB, = ax.plot([], [], 'o', color='red', markersize=10, label='Star B')
COM = ax.plot([0],[0], '+', color='white', markersize=5, label='COM')
ax.legend()

def orbit(r,mA,mB,M):
global x_valA, y_valA, x_valB, y_valB, thetaA, thetaB, G, dt

# Reset variables
x_valA, y_valA = [], []
x_valB, y_valB = [], []

M = mA+mB
rA = r*(mB/M)
rB = r*(mA/M)

#initial positions
positionA = np.array([rA * np.cos(thetaA), rA * np.sin(thetaA)]) # Star A initial position
positionB = np.array([rB * np.cos(thetaB), rB * np.sin(thetaB)]) # Star B initial position

# SIMULATION LOOP
for _ in range(numsteps):
# Store positions for both stars
x_valA.append(positionA[0])
y_valA.append(positionA[1])
x_valB.append(positionB[0])
y_valB.append(positionB[1])

# Update speed and angles for next positions
wA = sqrt(G*M/rA)
wB = sqrt(G*M/rB)
thetaA += wA * dt #update angle for starA
thetaB += wB * dt #update angle for starB

# Calculate new positions based on updated angles
positionA = np.array([rA * np.cos(thetaA), rA * np.sin(thetaA)])
positionB = np.array([rB * np.cos(thetaB), rB * np.sin(thetaB)])

#initialising the data
def init():
starA.set_data([], [])
starB.set_data([],[])
return starA, starB

def update(frame):
starA.set_data([x_valA[frame]], [y_valA[frame]]) # Pass as lists
starB.set_data([x_valB[frame]], [y_valB[frame]])
return starA, starB

ani = FuncAnimation(fig, update, frames=numsteps, init_func=init, blit=True, interval=50)
plt.title("Binary Star System")

def up(val):
global r, mA, mB, M
r = seperation_slider.val
mA = mA_slider.val
mB = mB_slider.val
M = mA + mB
orbit(r, mA, mB, M) # updated values into function
ani.event_source.stop() # Stop the current animation
ani.event_source.start() # Restart the animation with updated orbit
fig.canvas.draw_idle()

seperation_slider = Slider(ax=plt.axes([0.1, 0.00, 0.10, 0.04]), label='Radius', valmin=1, valmax=15, valinit=r, facecolor='w')
mA_slider = Slider(ax=plt.axes([0.45, 0.00, 0.15, 0.04]), label="Mass A", valmin=0.1, valmax=100, valinit=mA, facecolor='b')
mB_slider = Slider(ax=plt.axes([0.80, 0.00, 0.15, 0.04]), label="Mass B", valmin=0.1, valmax=100, valinit=mB, facecolor='r')

seperation_slider.on_changed(up)
mA_slider.on_changed(up)
mB_slider.on_changed(up)

orbit(r,mA,mB,M)

plt.show()
< /code>
Я попытался изменить константы из -за преобразования единиц. Подумал о введении некоторых реальных значений и сравнении значений скорости, но подумал, что это займет слишком много времени. Там.

Подробнее здесь: https://stackoverflow.com/questions/793 ... ulation-py