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

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


## Part (c):
##

## Plots around a_0 = 2:
##

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

par(mfrow=c(1, 2))

t.end = 20.0

diff_eq_params = list(a=1.5)
flowField(my_yprime, x.lim = c(0, 4), y.lim = c(2, 9), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('a= %.2f', diff_eq_params$a) )
trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
points(2, 4*diff_eq_params$a - 1.6, pch=19)

diff_eq_params = list(a=2.1)
flowField(my_yprime, x.lim = c(0, 4), y.lim = c(2, 9), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('a= %.2f', diff_eq_params$a) )
trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
points(2, 4*diff_eq_params$a - 1.6, pch=19)

par(mfrow=c(1, 1))

##dev.off()


## Plots for increasing values for a:
##
par(mfrow=c(1, 3))

t.end = 20.0

diff_eq_params = list(a=2.5)
flowField(my_yprime, x.lim = c(0, 20), y.lim = c(3, 22), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('a= %.2f', diff_eq_params$a) )
trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
points(2, 4*diff_eq_params$a - 1.6, pch=19)

diff_eq_params = list(a=3.0)
flowField(my_yprime, x.lim = c(0, 20), y.lim = c(3, 22), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('a= %.2f', diff_eq_params$a) )
trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
points(2, 4*diff_eq_params$a - 1.6, pch=19)

diff_eq_params = list(a=3.5)
flowField(my_yprime, x.lim = c(0, 20), y.lim = c(3, 22), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('a= %.2f', diff_eq_params$a) )
trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19)
points(2, 4*diff_eq_params$a - 1.6, pch=19)

par(mfrow=c(1, 1))