orbiter.py

from vpython import sphere, vector, rate, color, arrow, canvas, cross, triangle, vertex
from math import sqrt


G = 6.67e-11  # N kg^-2 m^2
au = 1.495978707e11
day = 60 * 60 * 24
year = day * 365.25

earth_mass = 5.972e24  # kg


def area_tri(center, prev, curr):
    return 0.5 * cross((prev - center), (curr - center)).mag


if __name__ == '__main__':

    scene = canvas(width=1200, height=1200)

    ship = arrow(pos=vector(au, 0, 0), make_trail=False)
    ship.r = vector(au, 0, 0)
    ship.m = 400e3
    ship.pos = ship.r
    ship.prev = ship.r

    sun = sphere(color=color.yellow, radius=0.1 * au)
    sun.r = vector(0, 0, 0)
    sun.m = 1.989e30  # kg
    sun.pos = sun.r

    d = (ship.r - sun.r).mag
    ship.velocity = vector(0, sqrt(G * sun.m / d), 0)
    ship.axis=ship.velocity.norm() * 0.05 * au
    area = 0

    t = 0
    dt = day/2

    while t < 2 * year:
        r = ship.r - sun.r
        f = G * sun.m * ship.m * r / r.mag ** 3
        if t < 20 * day:
            f += cross(r.norm(), vector(0, 0, 1)) * 1500
        ship.velocity -= f / ship.m * dt
        ship.r, ship.prev = ship.r + ship.velocity * dt, ship.r

        ship.pos, ship.axis = ship.r, ship.velocity.norm() * 0.05 * au

        triangle(v0=vertex(pos=sun.pos, color=color.blue), v1=vertex(pos=ship.prev), v2=vertex(pos=ship.pos))

        area += area_tri(sun.pos, ship.prev, ship.pos)

        print(area)

        t += dt
        rate(60)