import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# Load diagnostics
# -----------------------------

df = pd.read_csv("vgp_dynamics.csv")

v = df["AngularVelocity_deg_kyr"].values
k = df["Curvature"].values
age = df["Age_ka"].values

# Remove NaNs
mask = np.isfinite(v) & np.isfinite(k)

v = v[mask]
k = k[mask]
age = age[mask]

# -----------------------------
# Scatter phase diagram
# -----------------------------

plt.figure(figsize=(7,6))

sc = plt.scatter(
    v,
    k,
    c=age,
    s=25
)

plt.xlabel("Angular Velocity (deg/kyr)")
plt.ylabel("Pole Path Curvature")

plt.title("Geomagnetic Pole Phase Diagram")

cbar = plt.colorbar(sc)
cbar.set_label("Age (ka)")

plt.tight_layout()
plt.savefig("phase_velocity_curvature.png")
plt.show()