Skrypt liczący promień Ziemi
#!/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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment