# eatenLimit.py
# input 0.25, 0.5, 0.75
# The Lotka–Volterra predator–prey model with prey limits

import math
import matplotlib.pyplot as plt 
from frange import *


# birth/death rate of prey (rabbits)
alpha = 2 
beta =  2 

# birth/death rate of predators (foxes)
delta = 1
gamma = 1

dt  = 0.001
N_t = 20000   # no of steps

K = 20  # carrying capacity


ts = list(frange(0, N_t*dt, dt))
preds = [0]*N_t
prey = [0]*N_t

# population density
pop = float(input("Pop density? "))
prey[0] = pop
preds[0] = pop

# Step equations forward in time
for i in range(N_t-1):
  prey[i+1] = prey[i] + dt*alpha*prey[i]*(1 - prey[i]/K) - \
                        dt*beta*prey[i]*preds[i]

  preds[i+1] = preds[i] + dt*delta*prey[i]*preds[i] - \
                          dt*gamma*preds[i]

plt.plot(ts, preds, label="preds")
plt.plot(ts, prey, label="prey")
plt.legend()
plt.xlabel('time')
plt.ylabel('pop density')
plt.title("Limited Pred/Prey Pop. ("+str(pop)+")")
plt.show()