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

# The critical points of this system:
#
CP = data.frame(
    x=c(0, 3),
    y=c(0, 2)
)

my_jacobian = function(x, y, a=1, b=1) {
    # returns the Jacobian at the point (x, y)
    #
    J = matrix(c(a*(1-y/2), -a*x/2, b*y/3, b*(-1 + x/3)), 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] = parameters$a*xt*(1 - yt/2)
    dy[2] = parameters$b*yt*(-1.0 + xt/3)
    list(dy)
}

# Part (a):
#
#postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_5_prob_9_plot_part_a.eps', onefile=FALSE, horizontal=FALSE)

diff_eq_params = list(a=1, b=1)
t.end = 10.0
flowField(my_yprime, x.lim = c(-0.1, 7), y.lim = c(-0.1, 5), parameters = diff_eq_params, points = 30, add = FALSE, xlab='x', ylab='y')
trajectory(my_yprime, y0 = c(5, 2), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)

#dev.off()


# Part (b):
#
#postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_5_prob_9_plot_part_b.eps', onefile=FALSE, horizontal=FALSE)

par(mfrow=c(1, 2))

diff_eq_params = list(a=3, b=1)
t.end = 10.0
flowField(my_yprime, x.lim = c(-0.1, 7), y.lim = c(-0.1, 5), parameters = diff_eq_params, points = 25, add = FALSE, xlab='x', ylab='y', main='a=3; b=1')
trajectory(my_yprime, y0 = c(5, 2), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)

diff_eq_params = list(a=1/3, b=1)
t.end = 15.0
flowField(my_yprime, x.lim = c(-0.1, 7), y.lim = c(-0.1, 5), parameters = diff_eq_params, points = 25, add = FALSE, xlab='x', ylab='y', main='a=1/3; b=1')
trajectory(my_yprime, y0 = c(5, 2), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)

par(mfrow=c(1, 1))

#dev.off()


# Part (c):
#
#postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_5_prob_9_plot_part_c.eps', onefile=FALSE, horizontal=FALSE)

par(mfrow=c(1, 2))

diff_eq_params = list(a=1, b=3)
t.end = 10.0
flowField(my_yprime, x.lim = c(-0.1, 7), y.lim = c(-0.1, 5), parameters = diff_eq_params, points = 25, add = FALSE, xlab='x', ylab='y', main='a=1; b=3')
trajectory(my_yprime, y0 = c(5, 2), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)

diff_eq_params = list(a=1, b=1/3)
t.end = 15.0
flowField(my_yprime, x.lim = c(-0.1, 7), y.lim = c(-0.1, 5), parameters = diff_eq_params, points = 25, add = FALSE, xlab='x', ylab='y', main='a=1; b=1/3')
trajectory(my_yprime, y0 = c(5, 2), t.end = t.end, parameters = diff_eq_params, col='black', pch=8)
points(CP$x, CP$y, pch=19)

par(mfrow=c(1, 1))

#dev.off()