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

my_yprime = function(t, y, parameters) {
    xt = y[1]
    yt = y[2]
    zt = y[3]
    dy = rep(NA, length(y))
    dy[1] = parameters$sigma * (-xt + yt)
    dy[2] = parameters$r * xt - yt - xt*zt
    dy[3] = -parameters$b * zt + xt*yt
    list(dy)
}

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

diff_eq_params = list(sigma=10, b=8/3, r=21)

rm1 = diff_eq_params$r - 1
b_sqrt_rm1 = sqrt(diff_eq_params$b*rm1)

CP = data.frame(
    x=c(-b_sqrt_rm1, 0, b_sqrt_rm1),
    y=c(-b_sqrt_rm1, 0, b_sqrt_rm1),
    z=c(rm1, 0, rm1)
)

ts = seq(0, 30, length.out=2000)

res_1 = ode(y=c(3, 8, 0), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_1 = res_1[, 2] # the numerical function x(t)

res_2 = ode(y=c(5, 5, 5), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_2 = res_2[, 2]

res_3 = ode(y=c(5, 5, 10), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_3 = res_3[, 2]

Xs = cbind(xs_1, xs_2, xs_3)

matplot(ts, Xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)')
abline(h=CP$x, lwd=3, col='gray')
grid()
legend('topright', c('(3, 8, 0)', '(5, 5, 5)', '(5, 5, 10)'), lty=c(1, 1, 1), lwd=2, col=c('black','red','green') )

#dev.off()


# Part(b-c):
#
diff_eq_params = list(sigma=10, b=8/3, r=22) # pick a new value for x

rm1 = diff_eq_params$r - 1
b_sqrt_rm1 = sqrt(diff_eq_params$b*rm1)

CP = data.frame(
    x=c(-b_sqrt_rm1, 0, b_sqrt_rm1),
    y=c(-b_sqrt_rm1, 0, b_sqrt_rm1), z=c(rm1, 0, rm1)
)

ts = seq(0, 100, length.out=2000)

res_1 = ode(y=c(3, 8, 0), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_1 = res_1[, 2] # the numerical function x(t)

res_2 = ode(y=c(5, 5, 5), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_2 = res_2[, 2]

res_3 = ode(y=c(5, 5, 10), times=ts, fun=my_yprime, parms=diff_eq_params)
xs_3 = res_3[, 2]

Xs = cbind(xs_1, xs_2, xs_3)

matplot(ts, Xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)', main=sprintf('r= %.2f', diff_eq_params$r))
abline(h=CP$x, lwd=3, col='gray')
grid()
legend('topright', c('(3, 8, 0)', '(5, 5, 5)', '(5, 5, 10)'), lty=c(1, 1, 1), lwd=2, col=c('black','red','green') )