# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# Import numpy and matplotlib
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline

# Create model grid
x = np.linspace(0,50,51)

# +
# Define initial elevation along hillslope
zi = np.zeros_like(x)
for i in np.arange(1,len(zi)):
    if i <= (len(zi)/2.):
        zi[i] = zi[i-1] + 2
    else:
        zi[i] = zi[i-1] - 2
        
#Define initial bedrock height along hillslope
# wi = np.zeros_like(x)
# for i in np.arange(1,len(wi)):
#     if i <= (len(wi)/2.):
#         wi[i] = wi[i-1] + 9.5
#     else:
#         wi[i] = wi[i-1] - 9.5
wi = zi - 0.2
        
# Plot initial elevation of model to make sure it looks good
plt.plot(x,zi,'k')
plt.plot(x,wi,'r')


# -

def evolve(x,zi,wi,k,aw,bw,t,dt,u,bd):
    # Calculate dx
    dx = x[1]-x[0]
    # Create array so store elevation at each timestep
    z_out = np.zeros((int((t/dt)+1),len(x)))
    # Create array to store bedrockelevation at each timestep
    w_out = np.zeros((int((t/dt)+1),len(x)))
    # Store initial Elevation as timestep 0 in array
    z_out[0] = zi
    # Store initial bedrock Elevation as timestep 0 in array
    w_out[0] = wi
    # Create array of times at different timesteps
    t = np.linspace(dt,t,t/dt)
    t = np.concatenate(([0],t))
    # Set Elevation at boundaries
    z_out[:,0] = zi[0]
    z_out[:,-1] = zi[0]
    # Set Bedrock Elevation at boundaries
    w_out[:,0] = wi[0]
    w_out[:,-1] = wi[0]
    #Create loop to iterate through timesteps for diffusion
    for t in np.arange(1,len(t)):
        # Create loop to iterate through model grid
        for l in np.arange(0,len(x)):
            th = z_out[t-1][l] - w_out[t-1][l]
            w_r = aw*np.exp(-bw*th)
            w_out[t][l] = w_out[t-1][l] - w_r*dt + u*dt   
        for l in np.arange(1,len(x)-1):
            q_l = k*((z_out[t-1][l-1]-z_out[t-1][l])/dx)
            q_r = k*((z_out[t-1][l+1]-z_out[t-1][l])/dx)
            z_out[t][l] = z_out[t-1][l] + q_r*dt + q_l*dt + u*dt
        for l in [0,len(x)-1]:
            z_out[t][l] = z_out[t][l] - bd*dt
    return z_out,w_out


z_out,w_out = evolve(x,zi,wi,0.05,0.5,5,10000,1,0,-1)

plt.plot(z_out[5000])
plt.plot(w_out[5000])

dt = 0.5
t = 100
np.linspace(dt,t,t/dt)







StarcoderdataJupyterScriptsDedupFiltered