<jupyter_start><jupyter_text>Binary System of Sirius --- Question: How will the orbit be like if Sirius B experiences a sudden accleration?---We pick a binary star system that we are most familiar with - our neighbor Sirius A and Sirius B. Most stars have a companion, the gravitational force, serving as the centripital force, keeps them moving in orbit around the center of mass. When a star accelarates, the centripetal force it needs to keep in orbit would increase, so the gravitational force is not adequate to keep the star in its original orbit, so the star would actually leave its companion. However, if we assume there is friction or radiation -- something that wastes the kineetic energy, the velocity would go down and the star would be dragged back to its companion. In this model, we would simulate the orbit when Sirius B tries to escape and get pulled back.<jupyter_code># Units 
s = UNITS.second
N = UNITS.newton
kg = UNITS.kilogram
m = UNITS.meter
km = UNITS.kilometer

r_0 = 2.992e9 * km
T_0 = 50 * 365 * 24 * 60**2 * s
init = State(xA = -r_0 *2/3,
             yA = 0 *km,
             vxA = 0 *km/s,
             vyA = 4 / 3*pi*r_0/T_0,
             xB = r_0*1/3,
             yB = 0 * km,
             vxB = 0 *km/s,
             vyB = - 2 * pi * r_0 /3 / T_0)<jupyter_output><empty_output><jupyter_text>The mass of Sirius A is 2.1 times the mass of the sun, rounded to 2.The mass of Sirius B is 0.98 times the mass of the sun, rounded to 1.Therefore the center of mass is at 1/3 of the distance from Sirius A to B.<jupyter_code>m_sun = 1.989*10e30*kg
r_sun = 6.955*10e5 *km
systemA = System(init = init,
                 G = 6.674e-5 * N / kg**2 * km**2,
                 m = 2*m_sun,
                 r = 1.71 * r_sun,
                 t_0 = 0 *s
                 )

systemB = System (init = init,
                 r= 0.008*r_sun,
                 G = 6.674e-5 * N / kg**2 * km**2,
                 m = m_sun,
                 t_0 = 0*s)<jupyter_output><empty_output><jupyter_text>To accout for the deceleration of stars, we need to consider the drag forceF_drag = 1/2 * rho * mu^2 * Cd * A F_darg is the drag force, which is by definition the force component in the direction of the flow velocity;  rho is the is the mass density of the fluid;  mu is  the flow velocity relative to the object;  A is the reference area;  Cd is the drag coefficient, mostly depends on the Reynolds Number.<jupyter_code>def drag(state, velocity, system):
    rho = 1.67 *10e-12 * kg/km**3
    #Based on research, the average density in outer space is 1 hydrogen atom per 1 cubic cm.
    mu = velocity
    A = 2 * pi * system.r**2
 
    Cd = 2/3*rho**2
    drag_force = rho * mu**2 *Cd * A/2 * km**8/kg**2 /s**2
    
    return drag_force
    
drag (init, 10e7, systemA)
def universal_gravity(state, systemA, systemB):    
    """"
    Calculating gravitational force between Sirius A and Sirius B
    Taking the center of mass as origin
    xA, yA: posistion of Sirius A 
    xB, yB: position of Sirius B
    vxA, vyA: velocity of Sirius A
    vxB, vyB: velocity of Sirius B
    """
    xA, yA, vxA, yxA, xB, yB, vxB, vyB = state
    
    posA = Vector (xA, yA)
    posB = Vector (xB, yB)
    
    rA = posA.mag
    rB = posB.mag
    directionA = posA.hat()
    directionB = posB.hat()
    forceA = systemA.G * systemA.m *systemB.m / (rA+rB)**2 * directionA
    forceB = systemB.G * systemA.m *systemB.m / (rA+rB)**2 * directionB

    return rA, rB, forceA, forceB
    
universal_gravity(init, systemA, systemB)
def centri_force(state, systemA, systemB):
    """
    Calculating the current cetripetal force needed, related to current velocity and radius.
    
    """
    xA, yA, vxA, vyA, xB, yB, vxB, vyB = state
    
    velA = Vector(vxA, vyA)
    velB = Vector(vxB, vyB)
    posA = Vector (xA, yA)
    posB = Vector (xB, yB)
    
    vA = velA.mag
    vB = velB.mag
    rA = posA.mag
    rB = posB.mag
    
    cenA = systemA.m * vA**2 / rA
    cenB = systemB.m * vB**2 / rB
    
    cenA = cenA * (-posA.hat())
    cenB = cenB * (-posB.hat())
    
    return cenA, cenB
def comparing_force():
    """
    Updating distance and velocity to keep balance by comparing the gravitational force and the centripetal force required.
    """
    
    
#def event_func(state):

def acceleration_force(state, systemA, systemB ):
    
    xA, yA, vxA, vyA, xB, yB, vxB, vyB = state
    rA, rB,  grav_A, grav_B = universal_gravity(state, systemA, systemB)
    vA = sqrt(grav_A * rA / systemA.m)
    vB = sqrt(grav_B * rB / systemB.m)
    
   
    
    vel_A = Vector (vxA, vyA)
    vel_B = Vector (vxB, vyB)
    
    delta_vA = sqrt(rA * drag(state, vA, systemA)/systemA.m)
    delta_vB = sqrt(rB * drag(state, vB, systemB)/systemB.m)
    
    vA = vA-delta_vA
    vB = vB-delta_vB
    
    vA = vA * vel_A.hat()
    vB = vB * vel_B.hat()
    
    return vA, vB
def slope_func(state, t, systemA, systemB):
    xA, yA, vxA, vyA, xB, yB, vxB, vyB = state
    vA, vB = acceleration_force(state, systemA, systemB)
    
    
    return vA.x, vA.y,<jupyter_output>C:\Miniconda3\lib\site-packages\ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in sqrt
  """
C:\Miniconda3\lib\site-packages\pint\quantity.py:1236: RuntimeWarning: invalid value encountered in sqrt
  out = uf(*mobjs)
{'fraction_non_alphanumeric': 0.10066100094428707, 'fraction_numerical': 0.023040604343720492, 'mean_word_length': 2.4727868852459016, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, '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': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '303835', 'license': 'MIT', 'n_tokens_mistral': 1996, 'n_tokens_neox': 1875, 'n_words': 659, 'path': 'code/Project3-Binary.ipynb', 'repo_name': 'Alexis-Xinyi/ModSimPy', 'size': '5295'}