phy_542.py

import sys
import math

g = 9.8
r = 6400e3

# problem 5.42 from Savchenko

def e0(v):
    # total energy
    return v * v / 2 - g * r

def xm(v):
    # maximum distance from the Earth centre
    return -g*r*r / e0(v)

def gamma(x):
    # a helper function
    return 2*g*r*r/x

def tme(v):
    # half of the travel time
    x = xm(v)
    assert x >= r
    return (1/gamma(x)**0.5) * (x*math.acos((r/x)**0.5) + (r*(x-r))**0.5)

def main(tdays=7):
    nmax = 1000
    epsmax = 1e-6

    vmax = (2*g*r) **0.5
    vmin = (g*r) ** 0.5

    t = tdays * 24 * 3600

    for i in range(nmax):
        v = 0.5 * (vmin + vmax)
        ti = tme(v)

        eps = abs(t - ti)
        if  eps <= epsmax:
            break

        if ti < t:
            vmin = v
        else:
            vmax = v

    print(v, eps)
    return v
    
print(main(tdays=30./2.) - main(tdays=7./2.))