for pibal tracking (tenerife fieldtrip), encode the equations that convert Lift Rate, Azimuth, Elevation, TimeDifference into wind speeds using trigonometry

sympy_calculations.py

"""
WIND SPEED = f(L, b_w, az₀, el₀, az₁, el₁, T₀, T₁)
                  ____________________________________________________________
                 ╱   ⎛                         ⎛  ⎛az₀   az₁⎞⎞              ⎞
                ╱    ⎜      2       2⋅T₀⋅T₁⋅cos⎜π⋅⎜─── - ───⎟⎟         2    ⎟
               ╱     ⎜    T₀                   ⎝  ⎝180   180⎠⎠       T₁     ⎟
              ╱    L⋅⎜─────────── - ────────────────────────── + ───────────⎟
             ╱       ⎜   2⎛π⋅el₀⎞        ⎛π⋅el₀⎞    ⎛π⋅el₁⎞         2⎛π⋅el₁⎞⎟
            ╱        ⎜tan ⎜─────⎟     tan⎜─────⎟⋅tan⎜─────⎟      tan ⎜─────⎟⎟
           ╱         ⎝    ⎝ 180 ⎠        ⎝ 180 ⎠    ⎝ 180 ⎠          ⎝ 180 ⎠⎠
-1.389⋅   ╱        ──────────────────────────────────────────────────────────
         ╱                                         2/3
       ╲╱                                 (L + b_w)
───────────────────────────────────────────────────────────────────────────────
                                    T₀ - T₁
Where:
- L: the lift of the baloon (g)
- b_w: the weight of the balloon (g)
- az₀: the azimuth at time 0 (deg)
- el₀: the elevation at time 0 (deg)
- az₁: the azimuth at time 0 (deg)
- el₁: the elevation at time 0 (deg)
- T₀: time 0 (mins) 
- T₁: time 1 (mins) 

TODO: find the inverse
az₀, el₀, az₁, el₁ = f(L, b_w, T₀, T₁, ws) 
"""
import sympy as sy
import numpy as np 


# INPUT DATA
# ----------
# time (in minutes)
T_0, T_1 = sy.symbols("T_0 T_1")

# lift & balloon weight
L, b_w = sy.symbols("L b_w")
# elevation
el_0, el_1 = sy.symbols("el_0 el_1")
# azimuth
az_0, az_1 = sy.symbols("az_0 az_1")

# CALCULATIONS
# ------------
time_delta = 60 * (T_1 - T_0)
# convert to rads
el_rads_0 = el_0 * (sy.pi / 180)  # mpmath.radians(el_0)
el_rads_1 = el_1 * (sy.pi / 180)  # mpmath.radians(el_1)

# convert to rads
az_rads_0 = az_0 * (sy.pi / 180)  # mpmath.radians(az_0) 
az_rads_1 = az_1 * (sy.pi / 180)  # mpmath.radians(az_1) 

# calculate vertical distance from lift rate and time taken
LR = V = (83.34 * sy.sqrt(L)) / (sy.cbrt(L + b_w))

# calculate horizontal distance
H_0 = HcotE_0 = (LR * T_0) * sy.cot(el_rads_0)
H_1 = HcotE_1 = (LR * T_1) * sy.cot(el_rads_1)

# calculate adj & opp sides of triangle w.r.t. origin (observer)
adj_0 = cosAHcotE_0 = sy.cos(az_rads_0) * H_0
adj_1 = cosAHcotE_1 = sy.cos(az_rads_1) * H_1
opp_0 = sinAHcotE_0 = sy.sin(az_rads_0) * H_0
opp_1 = sinAHcotE_1 = sy.sin(az_rads_1) * H_1

# calculate the adj. & opp. of the difference triangle (m)
#  w.r.t each other (T_0, T_1)
delta_adj = adj_1 - adj_0
delta_opp = opp_1 - opp_0

# calculate the distance between T_0, T_1 (m)
hyp = sy.sqrt((delta_opp**2) + (delta_adj**2))
wind_speed = hyp / time_delta

# evaluate/print the outcome
sy.init_printing()
print(sy.simplify(wind_speed))

# substitute in the real values
# 20190501 - 07:15
data = {
    "b_w": 3,               # g
    "L": 7.9,               # g
    "az_0": 110.90,         # deg    
    "az_1": 117.30,         # deg    
    "el_0": 34.40,          # deg    
    "el_1": 54.30,          # deg    
    "T_0": 0.5,             # min
    "T_1": 1,               # min
}
answer = (
    wind_speed                      # m/s
    .subs(b_w, data["b_w"])         # g
    .subs(L, data["L"])             # g
    .subs(az_0, data["az_0"])       # deg    
    .subs(az_1, data["az_1"])       # deg    
    .subs(el_0, data["el_0"])       # deg    
    .subs(el_1, data["el_1"])       # deg    
    .subs(T_0, data["T_0"])         # min
    .subs(T_1, data["T_1"])         # min
).evalf()

# answer should be 5.502272
assert np.isclose(float(answer), 0.29, atol=0.1)

v_test = (
    V
    .subs(b_w, data["b_w"])         # g
    .subs(L, data["L"])             # g
).evalf()
assert np.isclose(float(v_test), 105.6473)

h_test = (
    H_0
    .subs(b_w, data["b_w"])         # g
    .subs(L, data["L"])             # g
    .subs(el_0, data["el_0"])       # deg    
    .subs(T_0, data["T_0"])         # min  
).evalf()

(
    hyp
    .subs(b_w, 3)           # g
    .subs(L, 6.3)           # g
    .subs(az_0, 110.90)      # deg    
    .subs(az_1, 117.30)      # deg    
    .subs(el_0, 34.40)       # deg    
    .subs(el_1, 54.30)       # deg   
    .subs(T_0, 0.5)         # min
    .subs(T_1, 1)           # min 
).evalf()


#  adj_0 == -27.47
adj_0_test = (
    cosAHcotE_0
    .subs(b_w, data["b_w"])         # g
    .subs(L, data["L"])             # g
    .subs(az_0, data["az_0"])       # deg    
    .subs(az_1, data["az_1"])       # deg    
    .subs(el_0, data["el_0"])       # deg    
    .subs(el_1, data["el_1"])       # deg    
    .subs(T_0, data["T_0"])         # min
    .subs(T_1, data["T_1"])         # min
).evalf()
assert np.isclose(float(adj_0_test), -27.47, atol=0.1)

opp_0_test = (
    sinAHcotE_0
    .subs(b_w, data["b_w"])         # g
    .subs(L, data["L"])             # g
    .subs(az_0, data["az_0"])       # deg    
    .subs(az_1, data["az_1"])       # deg    
    .subs(el_0, data["el_0"])       # deg    
    .subs(el_1, data["el_1"])       # deg    
    .subs(T_0, data["T_0"])         # min
    .subs(T_1, data["T_1"])         # min
).evalf()
assert np.isclose(float(opp_0_test), 72.15, atol=0.1)

# invert
# sy.solve(wind_speed, el_0, el_1, az_0, az_1)
# sy.solveset(wind_speed, el_0)
sy.simplify(wind_speed)