Lab6

Web VPython 3.2
from visual import *

G = 6.7e-11



# sun and earth

sun = sphere(pos = vec(0,0,0), radius = 7e9, color = color.yellow, make_trail = True, mass = 5.972e24)

earth = sphere(pos = vec(1.49e11,0,0), radius = 6.4e9, color = color.blue, make_trail = True, mass = 5.972e24)




sun.v = vec(0,0,0)
sun.p = sun.mass * sun.v

earth.v = vec(0,2.978e4,0)
earth.p = earth.mass * earth.v

old = earth.v
T = 0

t = 0
dt = 3600

while t < 2*365*24*60*60: 
    rate(4000)
    
    r = earth.pos - sun.pos
    r_mag = mag(r)
    FG = -G * sun.mass * earth.mass / r_mag**2 * r.hat
    
    
    earth.p += FG * dt
    earth.v = earth.p / earth.mass
    
    sun.p -= FG * 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 =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

    
    
    
    
    '''
    r = earth.pos - sun.pos
    r_mag = mag(r)
    FG = -G * sun.mass * earth.mass / r_mag**2 * r.hat
    
    
    earth.v += FG / earth.mass * dt
    earth.pos += earth.v * dt
    
    earth.p = earth.mass * earth.v
    sun.p = sun.mass * sun.v
    
    t += dt
    
    if mag(earth.pos - vec(1.49e11,0,0)) < 1e8 and T == 0:
        T = t / (3600 * 24)
        print (f"Period of Orbid: {T} days")
    '''