import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf

# Parameters
D = 2.0  # Diffusion coefficient in m^2/year
uplift_events = [(0, 14), (298, 3.93), (1000, 7.27)]  # (time, height) tuples
times = [10, 300, 400, 1005]  # Times at which to plot the scarp evolution
labels = ['First Event (1505 CE)', 'Second Event (1803 CE)', 'Reported Year (2017 CE)', 'Future Event (2803 CE)']
x = np.linspace(-100, 100, 200)  # Spatial domain

# Function to calculate the height profile after an uplift event
def scarp_height(x, t, uplift_time, uplift_height, D):
    elapsed_time = t - uplift_time
    if elapsed_time < 0:
        return np.zeros_like(x)
    else:
        return uplift_height / 2 * (1 - erf(x / (2 * np.sqrt(D * elapsed_time))))

# Plot
plt.figure(figsize=(10, 4))
for t, label in zip(times, labels):
    total_height = np.zeros_like(x)
    for uplift_time, uplift_height in uplift_events:
        total_height += scarp_height(x, t, uplift_time, uplift_height, D)
    # plt.plot(x, total_height, label=f't = {t} years')
    plt.plot(x, total_height, label=label)

plt.xlabel('Distance (m)')
plt.ylabel('Height (m)')
plt.title('Evolution of Kaladungi Fault Scarp Over Time')
plt.legend()
plt.grid(True)
plt.show()