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

diff_eq_params = list(gamma=0.8)

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

my_jacobian = function(x, y, parameters=diff_eq_params) {
    # returns the Jacobian at the point (x, y)
    #
    J = matrix(c(0, -1, -(x-0.15)*(x-2) - x*(x-2) - x*(x-0.15), -parameters$gamma), 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] = -yt
    dy[2] = -parameters$gamma * yt - xt * ( xt - 0.15 ) * ( xt - 2 ) 
    list(dy)
}

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

par(mfrow=c(1, 2))

# gamma = 0.8:
# 
diff_eq_params = list(gamma=0.8)

t.end = 7.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', main=sprintf('gamma= %.2f', diff_eq_params$gamma))

h = 0.01
trajectory(my_yprime, y0 = c(2-h, h), t.end = 2*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)


# gamma = 1.5:
#
diff_eq_params = list(gamma=1.5)

t.end = 7.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', main=sprintf('gamma= %.2f', diff_eq_params$gamma))

trajectory(my_yprime, y0 = c(2-h, h), t.end = 4*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)

par(mfrow=c(1, 1))

#dev.off()


if( FALSE ){ 
    # Part (c): find the gamma where the trajectory from (2, 0) goes to (0, 0)
    #
    diff_eq_params = list(gamma=1.3)

    t.end = 7.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', main=sprintf('gamma= %.2f', diff_eq_params$gamma))

    trajectory(my_yprime, y0 = c(2-h, h), t.end = 4*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)
}