# L-29 MCS 471 Fri 28 Oct 2022 : predpreysim.jl

# Simulation of a predator-prey model with the modified Euler method.

A = 1.1    # growth rate of the prey
B = -0.5   # correction to prey from encounter with predator
C = -0.75  # decay rate of the predator
D = 0.25   # correction to predator from encounter with prey

using Printf

"""
    modeulerpredprey(T::Float64,n::Int,
                     y1t0::Float64,y2t0::Float64,
                     verbose::Bool=true)

applies the modified Euler method to the predator-prey model
with time span [0, T], n steps, and initial values (y1t0, y2t0),
respectively for y1 and y2.
"""
function modeulerpredprey(T::Float64,n::Int,
                          y1t0::Float64,y2t0::Float64,
                          verbose::Bool=true)
    h = T/n
    if verbose
        println("  i   t         y1          y2 ")
    end
    (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))
    y1 = y1t0
    y2 = y2t0
    (rt[1], ry1[1], ry2[1]) = (0, y1, y2)
    for i=1:n
        t = i*h
        F1p = A*y1 + B*y1*y2
        F2p = C*y2 + D*y1*y2
        (y1p, y2p) = (y1 + h*F1p, y2 + h*F2p)
        F1c = A*y1p + B*y1p*y2p
        F2c = C*y2p + D*y1p*y2p
        y1 = y1 + (h/2)*(F1p + F1c)
        y2 = y2 + (h/2)*(F2p + F2c)
        (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)
        if verbose
            stri = @sprintf("%3d", i)
            strt = @sprintf("%.2f", t)
            stry1 = @sprintf("%.6e", y1)
            stry2 = @sprintf("%.6e", y2)
            println("$stri  $strt  $stry1  $stry2")
        end
    end
    return (rt, ry1, ry2)
end

"""
Computes one trajectory of the predator-prey simulation.
"""
function main()
    (endtime, nbpts) = (20.0, 100)
    (y1valt0, y2valt0) = (3.0, 2.0)
    println("Running the modified Euler method ...")
    t, y1, y2 = modeulerpredprey(endtime,nbpts,y1valt0,y2valt0,false)
    freq = 10
    idx = 1
    step = Int(nbpts/freq)
    for k=1:freq
        stri = @sprintf("%5d", idx)
        strt = @sprintf("%5.2f", t[idx])
        stry1 = @sprintf("%13.6e", y1[idx])
        stry2 = @sprintf("%13.6e", y2[idx])
        serr = @sprintf("%.2e", abs(y1[idx] - y1[1]))
        println("$stri  $strt  $stry1  $stry2  $serr")
        idx = idx + step
    end
end

main()