<jupyter_start><jupyter_text>QuestionHow does mass and initial velocity affect the shape of the orbit of a binary star?A binary star is a system containing two stars in close proximity.  In this system, both stars orbit around a common center or one star orbits around the other.  We wanted to create a model to understand how the masses of the two stars affected the shape of the stars orbit.  We also wanted to understand the shape of the binary star’s orbit when the star’s velocities are not constant.<jupyter_code>plt.figure(figsize = (15,15))
myImage = mpimg.imread("Untitled-1.png")
myImage = plt.imshow(myImage)
plt.axis('off')
plt.show<jupyter_output><empty_output><jupyter_text>Creating the ModelThe code below calculates the the initial velocity and distance from the origin of one star in a binary star system. To do so, Kepler's Third Law, Newton's Law of Universal Gravitation, and centripetal acceleration are used. The model itself also assumes that no other masses exert a force on the system. Calculating Initial Distance and Velocity from OriginModifying Kepler's Third Law gives the formula:$ D= \sqrt[3]{\frac{p^2 G (m_{1}+m_{2})} {4 \pi^2}} $This allows for the distance between stars to be calculated by hand picking masses for each sun and a desired period. From here, the radius from the origin to mass one can be calculated by the formula:$ r = \frac{m_{2} D} {m_{1}+m_{2}} $Then, to calculate the initial velocity, the Newton's Law of Universal Gravitation can be used to calculate the force exerted on each star. Dividing the given force by mass one gives the velocity of mass one, which then can be plugged into the formula:$ v = \sqrt{r a} $Repeating this process with the masses switched allows the same values for the second mass to be calculated.<jupyter_code>def init_r(p, m1, m2):
    """distance of m1 from origin
    
    p: orbital period
    m1: mass 1
    m2: mass 2
    
    returns: radius to origin of m1
    """
    G = 6.674e-11 * N / kg**2 * m**2
    
    d=((p**2)*G*(m1+m2)/(4*pi**2))**(1/3)
    r = m2*d/(m1+m2)
    return r

def init_v(p, m1, m2):
    """initial velocity of m1
    
    p: orbital period
    m1: mass 1
    m2: mass 2
    
    returns: initial velocity of m1
    """
    G = 6.674e-11 * N / kg**2 * m**2
    
    d=((p**2)*G*(m1+m2)/(4*pi**2))**(1/3)
    F = G*m1*m2/d**2
    a = F/m1
    r = d/2
    v = ((r)*a)**(1/2)
    return v<jupyter_output><empty_output><jupyter_text>Making the SystemThe make_sys function takes inputs for the mass of both stars as well as the orbital period, eccentricity, and end time. It then creates a system with those values as well as a State that hold the initial position and velocity of each star.<jupyter_code>def make_sys(sun1, sun2, p, e, t_end):
    """Computes gravitational force.
    
    sun1: mass of sun 1
    sun2: mass of sun 2
    re: distance from earth to origin
    ve: initial velocity of earth
    p: period of binary star orbit
    e: factor to multiply initial velocities of sun by
    
    returns: force vector
    """
    init = State(x1=init_r(p, sun1, sun2),y1=0,
                 vx1=0,vy1=-init_v(p, sun1, sun2)*e,
                 x2=-init_r(p, sun2, sun1),y2=0,
                 vx2=0, vy2=init_v(p, sun2, sun1)*e)
    
    system = System(init=init,
                    m1=sun1,
                    m2=sun2,
                    t_0=0 * s,
                    t_end=t_end)
    return system
sun1=1.989e30 * kg
sun2=1.989e30 * kg

p = 12.154e7 * s
e = 1
t_end = 12.154e8 * s

system = make_sys(sun1, sun2, p, e, t_end)<jupyter_output><empty_output><jupyter_text>Updating the StateThe fg function calculates a force vector based on the Law of Universal Gravitation to be used in the slope function. This makes the assumption that each planet is a point mass.<jupyter_code>def fg(r, m1, m2):
    """Computes gravitational force.
    
    r: distance between masses
    m1: mass 1
    m2: mass 2
    
    returns: force vector
    """
    G = 6.674e-11 * N / kg**2 * m**2

    direction = r.hat()
    radius = r.mag
    
    force = direction * G * m1 * m2 / radius**2
    return force<jupyter_output><empty_output><jupyter_text>The slope function below returns the derivatives of each value in the state, which means that it returns velocities and accelerations of each star in each direction. The derivative for position is velocity; these are already stored in the State. To calculate acceleration, the derivative of velocity, Newton's Law of Universal Gravitation:$ F = {\frac{G m_{1} m_{2}} {r}} $   As well as Newton's Second Law of Motion were used:                $ F = m a $By substituting for force, we get:$ m_{1} a = {\frac{G m_{1} m_{2}} {r}} $Then, divide by $m_{1}$ on both sides to get acceleration.$a = {\frac{G m_{2}} {r}} $The function below accomplishes this with the use of the fg function seen above.<jupyter_code>def slope_func(state, t, system):
    """Compute derivatives of the state.
    
    state: position, velocity
    t: time
    system: System object containing `g`
    
    returns: derivatives of y and v
    """
    x1, y1, vx1, vy1, x2, y2, vx2, vy2 = state
    unpack(system)    
    
    xy1 = Vector(x1,y1)
    xy2 = Vector(x2,y2)
    
    r1 = xy1.hat()*(xy1.mag+xy2.mag)
    r2 = xy2.hat()*(xy1.mag+xy2.mag)
    
    force1 = fg(r1, m1, m2)
    force2 = fg(r2, m1, m2)
    
    a1 = force1 / m1 * -1
    a2 = force2 / m2 * -1
    
    return (vx1, vy1, a1.x, a1.y, vx2, vy2, a2.x, a2.y)
results, details = run_ode_solver(system, slope_func, max_step=1000000)
details<jupyter_output><empty_output><jupyter_text>The plot_results function serves to quickly and cleanly plot our results.<jupyter_code>def plot_results(results):
    """Compute derivatives of the state.
    
    results: dataframe of binary star
    
    plots what you want
    """
    plot(results.x1, results.y1, 'orange', label = 'Star One')
    plot(results.x2, results.y2, 'r', label = 'Star Two')
    title = 'Binary Star Orbit'
        
    decorate(xlabel = 'X Position (m)',
             ylabel = 'Y Position (m)',
             title = title)
plot_results(results)<jupyter_output><empty_output><jupyter_text>Note: the orbits are overlapping, so only one line can be seen.  Results and Interpretations  System of Two Binary Stars with Identical Orbits and Identical MassesHere we run our model with identical masses and constant.<jupyter_code>sun1=1.989e30 * kg
sun2=1.989e30 * kg

p = 12.154e7 * s
e = 1
t_end = p

system = make_sys(sun1, sun2, p, e, t_end)

results, details = run_ode_solver(system, slope_func, max_step=1000000)
plot_results(results)<jupyter_output><empty_output><jupyter_text>As seen in the graph, both stars follow a circular orbit both identical in shape and position. Each star is on opposite to one another while orbiting. Because the orbits are circular, it means that the stars have a constant velocity throughout. This plot acts as a control as it follows expections of a binary star.  System of Two Binary Stars with Different MassesHere we run our model with constant velocities, but different masses.<jupyter_code>sun1=1.989e30 * kg
sun2=1e30 * kg

p = 12.154e7 * s
e = 1
t_end = p*5

system = make_sys(sun1, sun2, p, e, t_end)

results, details = run_ode_solver(system, slope_func, max_step=1000000)
plot_results(results)<jupyter_output><empty_output><jupyter_text>As seen in the graph, both stars follow a mostly circular orbit of different sizes. Each star is on opposite to one another while orbiting. Because the orbits are mostly circular, it means that the stars have a constant velocity throughout. It appears that the the smaller mass will have a larger orbit that encapsulates the larger mass' orbit. This makes sense because the smaller mass will have a greater acceleration than the larger mass. We think that the orbits are not perfectly aligned because of slight inaccuracies in either the initial velocity or position.  System of Two Binary Stars with Overlapping Orbits and Identical MassesHere we run our model with identical masses, but non-constant velocities.<jupyter_code>sun1=1.989e30 * kg
sun2=1.989e30 * kg

p = 12.154e7 * s
e = .6
t_end = p

system = make_sys(sun1, sun2, p, e, t_end)

results, details = run_ode_solver(system, slope_func, max_step=1000000)
plot_results(results)<jupyter_output><empty_output><jupyter_text>As seen in the graph, both stars follow ellpitical orbits. Each star is on opposite to one another while orbiting. Because the orbits are elliptical, it means that their velocities constantly change throughout the orbit. It appears that smaller than calculated initial velocities lead to a shorter period, while larger ones lead to longer periods. This makes sense because with increased initial velocities, the orbit becomes larger. With smaller or larger than calculated initial velocities, the orbits themselves begin to misalign.  System of Two Binary Stars with Overlapping Orbits and Different MassesHere we run our model with different masses and non-constant velocities.<jupyter_code>sun1=1.989e30 * kg
sun2=1.25e30 * kg

p = 12.154e7 * s
e = 1.2
t_end = 12.154e8 * s

system = make_sys(sun1, sun2, p, e, t_end)

results, details = run_ode_solver(system, slope_func, max_step=1000000)
plot_results(results)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08667962003454231, 'fraction_numerical': 0.03626943005181347, 'mean_word_length': 3.6025832091405863, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '286387', 'license': 'MIT', 'n_tokens_mistral': 3111, 'n_tokens_neox': 2901, 'n_words': 1188, 'path': 'code/Actual Final Project (all).ipynb', 'repo_name': 'clee4/ModSimPy', 'size': '9264'}