# Assess sensitivity of mescososm model to different alpha values, 

# Define differential equation with first order kinetics 
# and a power law residence time distribution
@sam_model("firstOrder", "powerLaw")

# Specify time span, background, initial condition, and parameters
tspan = (0.0, 1800) # model run time
h(p,t) = 310.0 # background specific conductance
u0 = [325.0] # concentration in channel after spike
τ₋ = 0.1
τ₊ = 900.0
k = 0.0
s = 27.0
sHZ = s/3
sCH = s - sHZ

# Solve model for two different alpha values
α = 1.5
qPerUnitHZ = qCalc_powerLaw(1, τ₋, τ₊, α)
p = [τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH]
prob2 = remake(prob, p = p)
sol1 = solve(prob2, MethodOfSteps(RK4()), saveat = 1)

α = 1.7
qPerUnitHZ = qCalc_powerLaw(1, τ₋, τ₊, α)
p = [τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH]
prob4 = remake(prob, p = p)
sol2 = solve(prob4, MethodOfSteps(RK4()), saveat = 1)

# Plot model output
plot(sol1, label = "alpha = 1.5", ylabel = "Nitrate Concentration (ppb N)", lw = 3, legend = :topleft)
plot!(sol2, label = "alpha = 1.7", xlabel = "Seconds", lw = 3)
