import numpy as np
import matplotlib.pyplot as plt

# =========================
# DATA: Dogger Bank RSL curve
# =========================
time_ka = np.array([
    30, 26, 21, 19, 16, 14.6, 13.5, 12.8, 11.7,
    10.5, 9.5, 9.0, 8.7, 8.2, 7.8, 7.5, 7.0, 6.0, 0
])

rsl_m = np.array([
    -115, -120, -120, -115, -100, -90, -70, -65, -60,
    -50, -40, -35, -30, -25, -15, -5, 0, 1.5, 0
])

# Uncertainty (only for 6 ka highstand)
rsl_min = rsl_m.copy()
rsl_max = rsl_m.copy()

# Mid-Holocene highstand range
rsl_min[time_ka == 6.0] = 1.0
rsl_max[time_ka == 6.0] = 2.0

# =========================
# SORT (descending time → ascending time)
# =========================
idx = np.argsort(time_ka)
time_ka = time_ka[idx]
rsl_m = rsl_m[idx]
rsl_min = rsl_min[idx]
rsl_max = rsl_max[idx]

# =========================
# INTERPOLATION (smooth curve)
# =========================
time_smooth = np.linspace(time_ka.min(), time_ka.max(), 500)
rsl_smooth = np.interp(time_smooth, time_ka, rsl_m)

# =========================
# PLOT
# =========================
plt.figure()

# Main curve
plt.plot(time_ka, rsl_m, marker='o', label="RSL (data)")
plt.plot(time_smooth, rsl_smooth, linestyle='--', label="Interpolated")

# Uncertainty band
plt.fill_between(time_ka, rsl_min, rsl_max, alpha=0.2)

# Axis formatting
plt.gca().invert_xaxis()  # 30 ka → present
plt.xlabel("Time (ka BP)")
plt.ylabel("Relative Sea Level (m)")
plt.title("Dogger Bank Relative Sea Level (30 ka → Present)")

# Reference line (modern sea level)
plt.axhline(0)

plt.legend()
plt.tight_layout()

# Save
plt.savefig("dogger_rsl_curve.png", dpi=150)

plt.show()

# =========================
# OPTIONAL: EXPORT CSV
# =========================
import pandas as pd

df = pd.DataFrame({
    "time_ka": time_ka,
    "rsl_m": rsl_m,
    "rsl_min": rsl_min,
    "rsl_max": rsl_max
})

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

print("Saved: dogger_rsl_curve.png, dogger_rsl_curve.csv")