Skrypt liczący promień Ziemi

earth.py

#!/usr/bin/env python
from mpmath import *
import sys

lim_d1 = (mpf(3050), mpf(3200))
lim_d2 = (mpf(12200), mpf(12800))
lim_h = (mpf(1.5), mpf(3.5))
lim_H = (mpf(69.5), mpf(73))
ratio = (mpf(50)/mpf(1180), mpf(120)/mpf(1110))

mp.dps = 100


def horizon_elev(R, h):
    return asin(R/(R+h))

def elev(R, src_h, distance, tgt_h):
    elev = atan((R+tgt_h)*sin(distance/R)/((R+src_h) - (R+tgt_h)*cos(distance/R)))
    if elev < 0:
        elev += pi
    return elev

def angles(R, h, H, d1, d2):
    horiz_elev = horizon_elev(R, h)
    bridge_elev = elev(R, h, d2, H)
    island_elev = elev(R, h, d1, 0)
    return horiz_elev-island_elev, bridge_elev-horiz_elev

def find_r_for_params(params):
    d1, d2, h, H, tgt_ratio = params
    Rmin = 2400000
    Rmax = 100000000

    while Rmax - Rmin > 1000:
        R = (Rmin + Rmax) / 2
        a, b = angles(mpf(R), h, H, d1, d2)
        ratio = a/b

        if ratio < tgt_ratio:
            Rmin = R
        elif ratio >= tgt_ratio:
            Rmax = R

    return R

def iter_params():
    for i in range(0,32):
        id1 = i % 2
        id2 = (i >> 1) % 2
        ih = (i >> 2) % 2
        iH = (i >> 3) % 2
        iratio = (i >> 4) % 2
        yield (lim_d1[id1], lim_d2[id2], lim_h[ih], lim_H[iH], ratio[iratio])

if __name__ == "__main__":
    min_radius = 1e9
    max_radius = 0
    for params in iter_params():
        r = find_r_for_params(params)
        if r < min_radius:
            min_radius = r
        if r > max_radius:
            max_radius = r

    print("Min radius:", min_radius)
    print("Max radius:", max_radius)