if( !require('phaseR') ){
    install.packages('phaseR')
}
library(phaseR)

# Parameters of the system:
# 
epsilon_1 = 1.5
sigma_1 = 0.5
alpha_1 = 1.0
epsilon_2 = 2
sigma_2 = 1.0
alpha_2 = 1.125

# The critical points of this system:
#
CP = data.frame(
    x=c(0, 0, epsilon_1/sigma_1, (epsilon_1*sigma_2 - epsilon_2*alpha_1)/(sigma_1*sigma_2 - alpha_1*alpha_2)),
    y=c(0, epsilon_2/sigma_2, 0, (epsilon_2*sigma_1 - epsilon_1*alpha_2)/(sigma_1*sigma_2 - alpha_1*alpha_2)))

my_jacobian = function(x, y) {
    # returns the Jacobian at the point (x, y)
    #
    J = matrix(c(epsilon_1 - 2*sigma_1*x - alpha_1*y, -alpha_1*x, -alpha_2*y, epsilon_2 - 2*sigma_2*y - alpha_2*x), nrow=2, ncol=2, byrow=TRUE)
}

As = mapply(my_jacobian, CP$x, CP$y, SIMPLIFY=FALSE)
for (ii in 1:length(As)){
    es = eigen(As[[ii]])
    print(sprintf('CP: x= %f; y= %f', CP$x[ii], CP$y[ii]))
    print('Jacobian=')
    print(As[[ii]])
    print('eigenvalues=')
    print(es$values)
}

my_yprime = function(t, y, parameters) {
    xt = y[1]
    yt = y[2]
    dy = rep(NA, length(y))
    dy[1] = xt*(parameters$epsilon_1 - parameters$sigma_1*xt - parameters$alpha_1*yt)
    dy[2] = yt*(parameters$epsilon_2 - parameters$sigma_2*yt - parameters$alpha_2*xt)
    list(dy)
}

diff_eq_params = list(epsilon_1=epsilon_1, sigma_1=sigma_1, alpha_1=alpha_1, epsilon_2=epsilon_2, sigma_2=sigma_2, alpha_2=alpha_2)

#postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_4_prob_3_plot.eps', onefile=FALSE, horizontal=FALSE)

t.end = 5.0
flowField(my_yprime, x.lim = c(-4, +4), y.lim = c(-4, +6), parameters = diff_eq_params, points = 30, add = FALSE, xlab='x', ylab='y')

# Plot some trajectories with initial condition near the critical points:
#
trajectory(my_yprime, y0 = c(-1.5, 1.5), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
trajectory(my_yprime, y0 = c(3.5, 4.0), t.end = t.end, parameters = diff_eq_params, col='blue', pch=8)
trajectory(my_yprime, y0 = c(-1.5, 4.5), t.end = t.end, parameters = diff_eq_params, col='red', pch=8)
trajectory(my_yprime, y0 = c(1.5, -1.5), t.end = t.end, parameters = diff_eq_params, col='green', pch=8)
points(CP$x, CP$y, pch=19)

#dev.off()