LAB6_FINISHED

Web VPython 3.2
from vpython import *

# --------------------- PLANETARY MOTION PROGRAM --------------------------

# Universal gravitational constant
G = 6.7e-11

# Create sphere objects for sun-Earth simulation:
sun = sphere(pos=vec(0, 0, 0), radius=7e9, color=color.yellow, make_trail=True, mass=1.989e30)
earth = sphere(pos=vec(1.49e11, 0, 0), radius=6.4e9, color=color.blue, make_trail=True, mass=5.972e24)

# Initialize:
sun.v = vec(0, 0, 0)  # Initial velocity approximating the Sun as stationary
sun.p = sun.mass * sun.v  # Initial momentum
earth.v = vec(0, 2.978e4, 0)  # Initial velocity of Earth
earth.p = earth.mass * earth.v  # Initial momentum of Earth
old = earth.v  # Old value of velocity to compare to updated value.
T = 0  # Initialize periodicity variable
t = 0  # Initialize time variable
dt = 3600  # Time step (1 hour)

    
# Velocity Components Graph
velocity_graph = graph(width=600, height=600, title="Velocity Components vs Time", xtitle="Time (s)", ytitle="Velocity (m/s)")
velocity_curve_x_earth = gcurve(color=color.blue, label="Earth x-velocity")
velocity_curve_y_earth = gcurve(color=color.red, label="Earth y-velocity")
velocity_curve_x_sun = gcurve(color=color.green, label="Sun x-velocity")
velocity_curve_y_sun = gcurve(color=color.orange, label="Sun y-velocity")

# Momentum Components Graph
momentum_graph = graph(width=600, height=600, title="Momentum Components vs Time", xtitle="Time (s)", ytitle="Momentum (kg m/s)")
momentum_curve_x_earth = gcurve(color=color.blue, label="Earth x-momentum")
momentum_curve_y_earth = gcurve(color=color.red, label="Earth y-momentum")
momentum_curve_x_sun = gcurve(color=color.green, label="Sun x-momentum")
momentum_curve_y_sun = gcurve(color=color.orange, label="Sun y-momentum")

# Total Momentum Graph
total_momentum_graph = graph(width=600, height=600, title="Total Momentum Components vs Time", xtitle="Time (s)", ytitle="Total Momentum (kg m/s)")
total_momentum_curve_x = gcurve(color=color.blue, label="Total x-momentum")
total_momentum_curve_y = gcurve(color=color.red, label="Total y-momentum")



# While Loop
while t < 2 * 365 * 24 * 3600:
    rate(4000)  # Simulation speed control

    # Calculate gravitational force:
    r_vec = earth.pos - sun.pos
    r_mag = mag(r_vec)
    F_gravity = -G * sun.mass * earth.mass / r_mag**2 * r_vec.hat

    # Update Earths momentum using P = mv:
    earth.p += F_gravity * dt 
    earth.v = earth.p / earth.mass  
    
    # Update Sun's momentum and velocity using P = mv:
    sun.p -= F_gravity * dt  
    sun.v = sun.p / sun.mass  
    
    # Update position:
    earth.pos += earth.v * dt

    # Update momentum for Sun (which is negligible but still updated):
    sun.p = sun.mass * sun.v
    

    # Update time:
    t += dt

    # Plot the velocity components:
    velocity_curve_x_earth.plot(pos=(t, earth.v.x))
    velocity_curve_y_earth.plot(pos=(t, earth.v.y))
    velocity_curve_x_sun.plot(pos=(t, sun.v.x))
    velocity_curve_y_sun.plot(pos=(t, sun.v.y))

    # Plot the momentum components:
    momentum_curve_x_earth.plot(pos=(t, earth.p.x))
    momentum_curve_y_earth.plot(pos=(t, earth.p.y))
    momentum_curve_x_sun.plot(pos=(t, sun.p.x))
    momentum_curve_y_sun.plot(pos=(t, sun.p.y))

    # Plot the total momentum:
    total_momentum_curve_x.plot(pos=(t, earth.p.x + sun.p.x))
    total_momentum_curve_y.plot(pos=(t, earth.p.y + sun.p.y))

    # Checking for periodicity to calculate the orbital period:
    if mag(earth.pos - vec(1.49e11, 0, 0)) < 1e8 and T == 0:  # Close to initial position again
        T = t / (3600 * 24)  # Convert time to days
        print(f"Orbital period: {T} days")  # Output orbital period