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 - 0.5*xt*yt
    dy[2] = -0.75*yt + 0.25*xt*yt
    list(dy)
}

# Parts (c):
#
ts = seq(0, 10, length.out=500)

res = ode(y=c(5, 2), times=ts, fun=my_yprime, parms=list()) # guess the initial conditions the book uses for its plot
xs = res[, 2] - mean(res[, 2]) # the numerical function x(t) (demeaned)
ys = res[, 3] - mean(res[, 3]) # the numerical function y(t) (demeaned)
T_x = my_period(ts, xs)
T_y = my_period(ts, ys)


# Part (d):
#
a = 1
alpha = 0.5
c = 0.75
gamma = 0.25

x0 = c/gamma # the critical point
y0 = a/alpha

small_amplitude_period = 2*pi/sqrt(a*c)
print(sprintf('small_amplitude period= %f', small_amplitude_period) )

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

As = seq(1.01, 1.6, length.out=10) # move the initial conditions (x, y) from the critical point (x0, y0) by multiplying by this factor
x_period_estimates = c()
y_period_estimates = c()
for( a in As ){
    res = ode(y=a*c(x0, y0), times=ts, fun=my_yprime, parms=list())
    xs = res[, 2] # the numerical function x(t)
    xs = xs - mean(xs) # (demeaned)
    ys = res[, 3] # the numerical function y(t)
    ys = ys - mean(ys) # (demeaned)

    T_x = my_period(ts, xs)
    T_y = my_period(ts, ys)

    x_period_estimates = c(x_period_estimates, T_x)
    y_period_estimates = c(y_period_estimates, T_y)
}

plot(As, x_period_estimates, type='l', col='blue', lty=1, lwd=2, pch=19, xlab='expansion factor from critical point', ylab='')
lines(As, y_period_estimates, type='l', col='red', lty=1, lwd=2, pch=19, xlab='expansion factor from critical point', ylab='')
abline(h=small_amplitude_period, col='green')
grid()

#dev.off()