"""
Weertman's marine ice-sheet instability (1974)
==============================================

Johannes Weertman studied the junction between an ice sheet grounded below
sea level and its floating shelf. The ice at the grounding line has the
flotation thickness, set by the water depth, and the shelf spreads at a rate
that rises as a power of that thickness. A grounding line on a bed that
deepens inland therefore cannot be steady: a small retreat thickens the ice
at the junction, which discharges more ice, which drives further retreat.
The example compares the snowfall feeding the ice sheet with a measure of
the discharge, the spreading of
10 km of cold shelf times the flotation thickness.
"""

import matplotlib.pyplot as plt
import numpy as np

from earthsciencekit.constants import GRAVITY, ICE_DENSITY, SECONDS_PER_YEAR
from earthsciencekit.glaciology import flotation_thickness, ice_shelf_strain_rate

x = np.linspace(300e3, 900e3, 601)  # m, distance from the divide
accumulation = 0.3 / SECONDS_PER_YEAR
supply = accumulation * x
fig, ax = plt.subplots()
for label, depth, sign in [("bed deepening seaward", 300 + 1e-3 * (x - 300e3), 1), ("bed deepening inland", 900 - 1e-3 * (x - 300e3), -1)]:
    thickness = flotation_thickness(depth, ICE_DENSITY, 1028)
    discharge = ice_shelf_strain_rate(thickness, 1e-25, 3, ICE_DENSITY, 1028, GRAVITY) * 10e3 * thickness
    crossing = np.argmin(np.abs(discharge - supply))
    imbalance_slope = np.gradient(discharge - supply, x)[crossing]
    assert np.sign(imbalance_slope) == sign  # positive: advance raises discharge (stable)
    line = ax.semilogy(x / 1e3, discharge * SECONDS_PER_YEAR / 1e6, label=f"discharge, {label}")[0]
    face = "white" if sign < 0 else line.get_color()
    ax.plot(x[crossing] / 1e3, discharge[crossing] * SECONDS_PER_YEAR / 1e6, "o", ms=9, zorder=3, color=line.get_color(), mfc=face)
ax.plot(x / 1e3, supply * SECONDS_PER_YEAR / 1e6, "k--", label="snowfall upstream")
ax.set(xlabel="grounding-line position [km]", ylabel="flux [km^2/year]", title="Weertman: open circle, an unstable grounding line")
ax.legend()

plt.show()