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]
    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)
}

##diff_eq_params = list(sigma=10, b=8/3, r=100)
##diff_eq_params = list(sigma=10, b=8/3, r=99.98)
##diff_eq_params = list(sigma=10, b=8/3, r=99.0)
diff_eq_params = list(sigma=10, b=8/3, r=100.5)


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

res_1 = ode(y=c(5, 5, 5), times=ts, fun=my_yprime, parms=diff_eq_params)
xs = res_1[, 2] # the numerical function x(t)
ys = res_1[, 3] # the numerical function y(t)
zs = res_1[, 4] # the numerical function z(t)

mask = ts > 10.0 # restrict to the longer term behaviour
ts = ts[mask]
xs = xs[mask]
ys = ys[mask]
zs = zs[mask]

## Estimate the periods:
##
T_x = my_period(ts, xs-mean(xs))
T_y = my_period(ts, ys-mean(ys))
T_z = my_period(ts, zs-mean(zs))
print(sprintf('r= %f; T_x= %f; T_y= %f; T_z= %f', diff_eq_params$r, T_x, T_y, T_z))

## Plot the solutions:
##
par(mfrow=c(3, 1))

plot(ts, xs, type='l', xlab='t', ylab='x(t)')
grid()

plot(ts, ys, type='l', xlab='t', ylab='y(t)')
grid()

plot(ts, zs, type='l', xlab='t', ylab='z(t)')
grid()

par(mfrow=c(1, 1))