import numpy as np

# Einheitsvektoren
unit_x = np.array([1, 0, 0])
unit_y = np.array([0, 1, 0])
unit_z = np.array([0, 0, 1])

def get_axis(axis):
    "Axis are numbered from 1 to 3 from x to z."
    match axis:
        case 1:
            ax = unit_x
        case 2:
            ax = unit_y
        case 3:
            ax = unit_z
        case _:
            ax = unit_x
    return ax

def proj(vec, axis: int =1):
    """Simple vector projection onto an axis."""    
    ax = get_axis(axis)
    return np.dot(vec, ax) * ax

def rotate(v, angle=90, axis=1):
    "Rotate a vector with an angle around a axis with the right hand rule."
    angle = angle/180 * np.pi
    k = get_axis(axis)

    return (
        v * np.cos(angle)
        + np.cross(k, v) * np.sin(angle)
        + k * np.dot(k, v) * (1 - np.cos(angle))
    )
        
def agl(a, b):
    "Get the angle between two vectors. This is always between 0 and 180 degree."
    return np.round(np.acos(np.dot(a, b)/(np.linalg.norm(a) * np.linalg.norm(b)))/(2 * np.pi) * 360)

def get_angles(source, target):
    """Main function to get the phi and theta angles for a source and a target vector. Both vectors must lie on the front half sphere.
       Phi is from 0 to 180 where 0 means left when you look at the mirrors. The hardware is bounded between 45 and 135 degree. Thus the here provided angle needs to be subtracted by 45 and then doubled.
       Theta is from 0 to 90 where 0 means up."""
    source_planar = source - proj(source, 3)
    target_planar = target - proj(target, 3)

    source_phi = agl(source_planar, unit_x)
    target_phi = agl(target_planar, unit_x)

    source_theta = agl(rotate(source, 90 - source_phi, 3), unit_z)
    target_theta = agl(rotate(target, 90 - target_phi, 3), unit_z)

    phi = None
    theta = None

    theta_diff = None
    phi_diff = agl(source_planar, target_planar)
    if source_phi < target_phi:
        rota = rotate(source_planar, phi_diff, 3)
        theta_diff = agl(rota, target)
        phi = source_phi + phi_diff/2
    else:
        rota = rotate(target_planar, phi_diff, 3)
        theta_diff = agl(rota, source)
        phi = target_phi + phi_diff/2

    if source_theta < target_theta:
        theta = target_theta + theta_diff/2
    else:
        theta = source_theta + theta_diff/2

    return (phi, theta)