#!/usr/bin/env python3

# here, we integrate Eq. (1) using a simple explicit Euler method 
# where we add random noise at each timestep
# we're also going to change the parameter b linearly as time progresses

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
plt.rcParams.update({'font.size':14})

# DEFINE FUNCTIONS
# function to calculate the stable fixed point
def stable_fixed_point(r,a,b):
    return a/2 + np.sqrt(a**2/4-b) 

# defining the model: i.e. f(N,r,a,b)
def model(N,params):
    r,a,b = params
    ndot = r*N*(-N**2+a*N-b)
    return ndot

# DEFINE PARAMETERS: this is what you might need to change
# define constants r and a, and the start and end values of b
r = 1
a = 8 
b_start = 14 
b_end = 15

# set up start and end times of the integration
tstart = 0
tend = 250

# timestep (recommend not touching this)
dt = 0.05

# amplitude of random noise 
amplitude = 0.6 

# INTEGRATE AND PLOT  
t = np.linspace(tstart,tend,int((tend-tstart)/dt))
b_ramp = np.linspace(b_start,b_end,len(t))
N = np.zeros(len(t))
N[0] = stable_fixed_point(r,a,b_ramp[1]) 

for i in range(1,len(t)):
    N[i] = N[i-1] + model(N[i-1],[r,a,[b_ramp[i]]])*dt + np.random.normal(0,amplitude*dt)

plt.plot(t,N)
plt.xlabel('Time')
plt.ylabel('N')
plt.show()