# Erdradius
r = 6371

from math import sin, cos, tan, radians, asin, acos, atan, atan2, degrees
import numpy as np

# Mit "nord" ist im Folgenden die nördliche Breite in Grad gemeint. 
# Mit "ost" ist im Folgenden die östliche Länge in Grad gemeint. 

def xyz_von_sphaerisch(nord, ost):
    z = r * sin(radians(nord))
    r_z = r * cos(radians(nord))
    x = r_z * cos(radians(ost))
    y = r_z * sin(radians(ost))
    return x, y, z

def sphaerisch_von_xyz(x, y, z):
    nord = asin(z/r)
    ost = atan2(y, x)
    return nord, ost

def daten_zur_orthodrome(A_nord, A_ost, B_nord, B_ost):
    A = np.array(xyz_von_sphaerisch(A_nord, A_ost))
    B = np.array(xyz_von_sphaerisch(B_nord, B_ost))
    print('Ausgaben in km bzw. Grad.')    
    print(f'Tunnellänge: {np.linalg.norm(B-A)}')
    alpha = acos(np.dot(A, B)/(np.linalg.norm(A)* np.linalg.norm(B)))
    d = alpha * r
    print(f'sphärische Distanz: {d}')
    print(f'Mittelpunktswinkel: {degrees(alpha)}')
    
    # Der Rest liefert nur manchmal die richtigen Werte (etwa im Beispiel Zürich-JFK);
    # schon im Beispiel JFK-Zürich versagt er. 
    # Vermutlich müsste man diverse Fallunterscheidungen betrachten.
    print('Rest nur manchmal richtig.')
    normalenvektor = np.cross(B, A)
    nordpol = np.array([0, 0, 1])
    beta = acos(np.dot(normalenvektor, nordpol)/(np.linalg.norm(normalenvektor))* np.linalg.norm(nordpol))
    print(f'Der nördlichster Punkt auf der Orthodrome hat die geographische Breite {degrees(beta)} N.')
    v = np.cross(normalenvektor, nordpol)
    w = np.cross(v, normalenvektor)
    w = r / np.linalg.norm(w) * w
    noerdlichst_nord, noerdlichst_ost = sphaerisch_von_xyz(w[0], w[1], w[2])
    print(f'Koordinaten des nördlichsten Punktes: {degrees(noerdlichst_nord)} N, {degrees(noerdlichst_ost)} E')
    return d

# Zürich
zuerich_nord=47.5
zuerich_ost=8.5

# JFK New York
jfk_nord = 40.6
jfk_ost = -73.8

daten_zur_orthodrome(zuerich_nord, zuerich_ost, jfk_nord, jfk_ost)