Skip to content

Instantly share code, notes, and snippets.

@fizyk20

fizyk20/earth.py Secret

Created Sep 5, 2018
Embed
What would you like to do?
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