import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

# ------------------------------------------------------------
# LOAD IERS FINALS FILE
# ------------------------------------------------------------

df = pd.read_csv("finals.all.csv", sep=";")

# Keep only final solutions
df = df[df["Type"] == "final"].copy()

# Build datetime from fields
df["date"] = pd.to_datetime(dict(year=df.Year, month=df.Month, day=df.Day))

# Extract pole coordinates (mas)
x = df["x_pole"].astype(float).values
y = df["y_pole"].astype(float).values
t = df["date"].values

# Remove NaNs
mask = np.isfinite(x) & np.isfinite(y)
x = x[mask]
y = y[mask]
t = t[mask]

# ------------------------------------------------------------
# PARAMETERS
# ------------------------------------------------------------

window = 365 * 3      # 3-year rolling window
step   = 30           # evaluate every ~month

angles = []
times  = []

# ------------------------------------------------------------
# ROLLING PCA ANISOTROPY EXTRACTION
# ------------------------------------------------------------

for i in range(0, len(x) - window, step):
    Xw = np.vstack([x[i:i+window], y[i:i+window]]).T

    pca = PCA(n_components=2)
    pca.fit(Xw)

    vec = pca.components_[0]

    angle = np.degrees(np.arctan2(vec[1], vec[0]))
    angle = angle % 360

    angles.append(angle)
    times.append(t[i + window//2])

angles = np.array(angles)
times  = np.array(times)

# ------------------------------------------------------------
# PLOT
# ------------------------------------------------------------

plt.figure(figsize=(11,5))

plt.plot(times, angles, color="black", lw=1.6)

# Reference longitudes
plt.axhline(285, color='green', linestyle='--', lw=2,
            label='Forcing Axis (~75°W)')
plt.axhline(31,  color='blue', linestyle='--', lw=2,
            label='Inertial Plane (~31°E)')
plt.axhline(170, color='red', linestyle='--', lw=2,
            label='Excursion Plane (~170°E)')

plt.ylabel("Anisotropy Bearing (°E)")
plt.title("Polar Motion Anisotropy Axis vs Time (IERS Finals Daily)")
plt.legend()
plt.tight_layout()
plt.show()
