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

# =========================
# PARAMETERS (from your fit)
# =========================
k = 0.2837249
p = 0.50465482
dz = 5.0

C = k / dz

# Topological regime (tuneable)
beta = 1.2
q = 1.3
A_c = 0.28

# =========================
# HYBRID ALPHA FUNCTION
# =========================
def alpha(A):
    if A > A_c:
        return 1 - k * (1 - A)**p
    else:
        return beta * A**q

# Vectorized version
alpha_vec = np.vectorize(alpha)

# =========================
# DISCRETE ITERATION
# =========================
def simulate_discrete(A0=1.0, steps=50):
    A = [A0]
    for _ in range(steps):
        a = alpha(A[-1])
        A.append(a * A[-1])
    return np.array(A)

# =========================
# CONTINUOUS INTEGRATION
# =========================
def dA_dz(A):
    if A > A_c:
        return -C * A * (1 - A)**p
    else:
        return -(1/dz) * (1 - beta * A**q) * A

def simulate_continuous(A0=1.0, z_max=200, dz_small=0.1):
    z_vals = np.arange(0, z_max, dz_small)
    A = np.zeros_like(z_vals)
    A[0] = A0

    for i in range(1, len(z_vals)):
        A[i] = A[i-1] + dA_dz(A[i-1]) * dz_small
        A[i] = max(A[i], 0)

    return z_vals, A

# =========================
# RUN
# =========================
A_disc = simulate_discrete()
z_cont, A_cont = simulate_continuous()

# =========================
# PLOTS
# =========================
plt.figure()
plt.plot(A_disc, label="Discrete (hybrid α)")
plt.plot(np.interp(np.arange(len(A_disc)), 
                   np.linspace(0, len(A_disc), len(A_cont)), 
                   A_cont), label="Continuous")
plt.xlabel("Step")
plt.ylabel("Area fraction A")
plt.legend()
plt.title("Dogger Bank Collapse: Hybrid Model")
plt.show()