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

# -----------------------------
# Energy index
# -----------------------------

energy = v * k

energy = np.nan_to_num(energy)

df["ExcursionEnergy"] = energy

df.to_csv("vgp_energy_index.csv", index=False)

# -----------------------------
# Plot
# -----------------------------

plt.figure()

plt.plot(age, energy)

plt.gca().invert_xaxis()

plt.xlabel("Age (ka)")
plt.ylabel("Excursion Energy Index")

plt.title("Geomagnetic Excursion Energy")

plt.tight_layout()
plt.savefig("excursion_energy_index.png")

plt.show()

# -----------------------------
# Identify strongest events
# -----------------------------

top = df.sort_values("ExcursionEnergy", ascending=False).head(10)

print("\nTop candidate instability events:\n")
print(top[["Age_ka","ExcursionEnergy"]])