[Python] Problem Updating Angular Velocity and Angles in Binary Star Orbit Simulation (Python,...

I'm creating a simulation for a binary star system in python using matplotlib. I want it to change based on the users inputs through sliders I've provided. When I first created this without any real equations for binary star systems (w = sqrt(G*M/r) and theta += w * dt) it worked perfectly. When the masses of each star were changed/the separation between them both, the individual distances of the stars would change too. But when I try to change the speed as well, the total separation of the stars doesn't stay constant. They should always be travelling opposite to one another. I can't find the exact problem and I have a feeling it can be a combination of multiple. I'm also trying to avoid using elliptical orbits.

My code:

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()


I've tried changing the constants due to unit conversions. Thought about inputing some real values in and comparing speed values but thought it would take too much time. I also thought about the rA and rB variables and how they're defined but I'm not really sure how to go from there.

Continue reading...