#
# Lotka-Volterra Equations
#
# Euler integration
# September 2012, Hoffmann
#

# libraries
import numpy
import matplotlib
import pylab

# model parameters
m = 2 # model selection
a = [0.00280, 0.00100, 0.00280][m]   
b = [0.00032, 0.00032, 0.00032][m]   
c = [0.00800, 0.00800, 0.00200][m] 
d = [0.00080, 0.00080, 0.00080][m] 
initX = [20, 20, 20][m]
initY = [50, 50, 50][m]

# difference equations
def deltaX(x,y,a,b):
    return x*(a - b*y)

def deltaY(x,y,c,d):
    return y*(-c+d*x)

# print every so distant values
def printSome(x,y,c,gap):
    if c % gap == 0:
        print("step",c,"-",x,"rabbits and",y, "foxes")

# plot values for rabbits and foxes
def drawCurve(yList,zList):
    n = len(yList)
    x = pylab.arange(n)
    pylab.plot(x, yList, 'b')
    pylab.plot(x, zList, 'r')
    pylab.show()

# main simulation
def main():
    gap = 500
    count = 0
    rabbits = minRab = initX
    foxes  = minFox = initY
    nSteps = 5000
    rabList = list()
    foxList = list()
    rDead = -1
    fDead = -1
    #
    for count in range(nSteps):
        # plot foxes and rabbits in different color
        rabList = rabList + [rabbits]
        foxList = foxList + [foxes]
        # printSome(rabbits,foxes,count,gap)
        delX = (deltaX(rabbits,foxes,a,b))
        delY = (deltaY(rabbits, foxes,c,d))
        rabbits = rabbits + delX
        foxes = foxes + delY
        if rabbits<0.001:
            rabbits = 0
            if rDead < 0:
                rDead = count
        if foxes<0.001:
            foxes = 0
            if fDead < 0:
                fDead = count
        if rabbits < minRab:
            minRab = rabbits
        if foxes < minFox:
            minFox = foxes
#        count = count+1
    print("Test case", m+1)
    print("min fox population: {0:.6f}".format(minFox), ", min rabbit pop: {0:.6f}".format(minRab))
    if rDead >= 0:
        print("rabbits died out at step",rDead)
    if fDead >= 0:
        print("foxes died out at step",fDead)
    drawCurve(rabList,foxList)

main()