Lab6

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)

    # Graph


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

    # Update time:
    t += dt

    # 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