#!/usr/bin/env python3

import numpy as np


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_separation(v1, v2):
    dot = np.clip(np.dot(v1, v2), -1.0, 1.0)
    return np.degrees(np.arccos(dot))


# Extracted axis from PCA
axis_lat = -6.8078
axis_lon = 103.4002

axis_vec = latlon_to_vec(axis_lat, axis_lon)

candidates = {
    "0° Meridian": (0.0, 0.0),
    "90°E Meridian": (0.0, 90.0),
    "31°E (TPW axis)": (0.0, 31.0),
    "170°E (Excursion)": (0.0, 170.0),
    "Sampling axis (~90°W)": (24.4033, -89.0999),
}

print("Extracted axis:", axis_lat, axis_lon)
print()

for name, (lat, lon) in candidates.items():
    vec = latlon_to_vec(lat, lon)
    angle = angular_separation(axis_vec, vec)
    print(f"{name}: {angle:.2f}°")
