from math import sqrt, acos

G = 6.67428e-11
earth_mass = 5.9736e24
earth_radius = 6.371e6

def potential_energy(m1, m2, r):
    return -G * m1 * m2 / r

def orbital_speed(m, r):
    return sqrt(G * m / r)

def orbital_kinetic_energy(m, m2, r):
    return m * orbital_speed(m2, r) ** 2 / 2

def free_fall_time(y0, y, m):
    return (sqrt(y0 ** 3 / (2 * G * m)) * 
            (sqrt(y / y0 * (1 - y / y0)) + acos(sqrt(y / y0))))
    


earth_surface_p = potential_energy(1.0, earth_mass, earth_radius)

ISS_orbit = 6.738e6

geostationary_orbit =  42.164e6

#print potential_energy(1.0, earthMass, ISS_orbit) - earth_surface_p
#print orbital_kinetic_energy(1.0, earthMass, ISS_orbit)

print potential_energy(1.0, earth_mass, geostationary_orbit) - earth_surface_p
print orbital_kinetic_energy(1.0, earth_mass, geostationary_orbit)

print free_fall_time(ISS_orbit, earth_radius, earth_mass) # wow, that's fast!
print free_fall_time(geostationary_orbit, earth_radius, earth_mass) / ( 60 * 60)