Isoclines, quiver, etc

isoclines.jl

using Plots
using Iterators
pyplot()

function dn(n1, n2; α=[1.0 0.0; 0.0 1.0], r1=1.0, r2=1.0)
    dn1 = n1*(r1-α[1,1]*n1-α[1,2]*n2)
    dn2 = n2*(r2-α[2,2]*n2-α[2,1]*n1)
    return dn1, dn2
end

function equil(α=[1.0 0.0; 0.0 1.0], r1=1.0, r2=1.0)
    return[
        (0.0, 0.0),
        (r1/α[1,1], 0.0),
        (0.0, r2/α[2,2]),
        (-(α[2,2]*r1-α[1,2]*r2)/(α[1,2]*α[2,1]-α[1,1]*α[2,2]),(α[2,1]*r1-α[1,1]*r2)/(α[1,2]*α[2,1]-α[1,1]*α[2,2]))
        ]
end

function roc(x...; α=[1.0 0.0; 0.0 1.0], r1=1.0, r2=1.0)
    dn1, dn2 = dn(x..., α=α, r1=r1, r2=r2)
    return sqrt(dn1^2+dn2^2)
end

n = linspace(0.0, 1.25, 50)
ncoarse = linspace(0.02, 1.23, 18)
ngr = collect(product(ncoarse, ncoarse))

α = [1.1 0.9; 0.8 1.02]
# Define some functions for this matrix
xroc = (x1, x2) -> roc(x1, x2, α=α)
function dnquiver(x1, x2)
    dn1, dn2 = dn(x1, x2, α=α)
    rc = xroc(x1, x2)
    return (dn1/rc)*0.02, (dn2/rc)*0.02
end
isocl1 = (x) -> -(α[1,1]/α[1,2])*x+1.0/α[1,2]
isocl2 = (x) -> -(α[2,1]/α[2,2])*x+1.0/α[2,2]

plot(n, n, xroc, c=:Blues, t=:contour, aspectratio=1, leg=false, f=true, levels=200, clims=(0,0.6), xlim=(0,Inf), ylim=(0, Inf), size=(800, 800))
plot!(n, isocl1, c=:black, lab="Species 1", leg=true, lw=2)
plot!(n, isocl2, c=:black, ls=:dash, lab="Species 2", lw=2)
quiver!(ngr, quiver=dnquiver, c=:grey)
xaxis!("Species 1")
yaxis!("Species 2")

savefig("quiver.png")