import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform

# ----------------------------------------
# Load embedded attractor state
# ----------------------------------------

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

omega = df["AngularVelocity_deg_kyr"].values
energy = df["ExcursionEnergy"].values
age = df["Age_ka"].values

# Normalize
omega = (omega - np.mean(omega)) / np.std(omega)

# ----------------------------------------
# Takens embedding
# ----------------------------------------

tau = 3
m = 3

N = len(omega) - (m - 1) * tau
X = np.zeros((N, m))

for i in range(m):
    X[:, i] = omega[i*tau : i*tau + N]

age_emb = age[:N]
energy_emb = energy[:N]

# ----------------------------------------
# Distance matrix
# ----------------------------------------

D = squareform(pdist(X))

# Recurrence threshold
eps = np.percentile(D, 10)

R = D < eps

# ----------------------------------------
# Plot recurrence map
# ----------------------------------------

plt.figure(figsize=(6,6))
plt.imshow(R, origin="lower", cmap="binary", aspect="auto")

plt.title("Geomagnetic Attractor Recurrence Plot")

plt.xlabel("Time index")
plt.ylabel("Time index")

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