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

source('../Chapter3/utils.R')

my_yprime = function(t, y, parameters) {
    xt = y[1]
    yt = y[2]
    dy = rep(NA, length(y))
    dy[1] = xt * ( 2.4 - 0.2 * xt - 2 * yt / (xt + 6) )
    dy[2] = yt * ( -0.25 + xt / (xt + 6) )
    list(dy)
}


my_jacobian = function(x, y, a=2.4) {
    ## returns the Jacobian at the point(x, y)
    ##
    J = matrix(c(a - 0.4*x - 2*y/(x+6) + 2*x*y/(x+6)^2, -2*x/(x+6), y/(x+6) - x*y/(x+6)^2, -0.25 + x/(x+6)), nrow=2, ncol=2, byrow=TRUE)
}

CP = data.frame(x=c(0, 12, 2), y=c(0, 0, 8))

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(round(es$values, 4))
}

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

t.end = 30.0

diff_eq_params = list()

flowField(my_yprime, x.lim = c(0, 13), y.lim = c(0, 13), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y')
trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
points(CP$x, CP$y, pch=19)

#dev.off()