Problem Updating Angular Velocity and Angles in Binary Star Orbit Simulation (Python, Matplotlib)

Question

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

from scipy.constants import pi

dt = 0.1
numsteps = 10000
pi = pi ## imported by from scipy.constants import 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")

# 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.

Answer 1

try changing:

# Update speed and angles for next positions
        wA = sqrt(G*M/rA)
        wB = sqrt(G*M/rB)

with :

# Update speed and angles for next positions
        wA = sqrt(G*M/rA*rA)
        wB = sqrt(G*M/rB*rB)

and let us know