#!/usr/bin/env python3

import numpy as np

PERMUTATIONS = 100000


def latlon_to_vec(lat_deg, lon_deg):
    lat = np.radians(lat_deg)
    lon = np.radians(lon_deg)
    return np.array([
        np.cos(lat) * np.cos(lon),
        np.cos(lat) * np.sin(lon),
        np.sin(lat)
    ])


def angular_sep(v1, v2):
    dot = np.clip(np.dot(v1, v2), -1.0, 1.0)
    return np.degrees(np.arccos(dot))


def random_unit_vector():
    v = np.random.normal(size=3)
    return v / np.linalg.norm(v)


if __name__ == "__main__":

    # Observed dipole regression axis
    obs_lat = 2.68139
    obs_lon = -74.13507

    obs_vec = latlon_to_vec(obs_lat, obs_lon)

    # Candidate Euler axes (equatorial)
    candidates = {
        "59W Euler": (0.0, -59.0),
        "121E Euler": (0.0, 121.0),
        "90E Hemisphere": (0.0, 90.0),
    }

    print("Observed axis:", obs_lat, obs_lon)
    print()

    for name, (lat, lon) in candidates.items():

        cand_vec = latlon_to_vec(lat, lon)
        observed_angle = angular_sep(obs_vec, cand_vec)

        # Monte Carlo
        count = 0
        for _ in range(PERMUTATIONS):
            rand_vec = random_unit_vector()
            angle = angular_sep(rand_vec, cand_vec)
            if angle <= observed_angle:
                count += 1

        p_value = count / PERMUTATIONS

        print(f"{name}:")
        print(" Angular separation:", observed_angle)
        print(" p-value:", p_value)
        print()
