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

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

## Parameters:
##
a = 1
alpha = 0.5
c = 0.75
gamma = 0.25

x_bar_linear = c/gamma # the linear approximation average
y_bar_linear = a/alpha
my_yprime = function(t, y, parameters){
    xt = y[1]
    yt = y[2]
    dy = rep(NA, length(y))
    dy[1] = a*xt - alpha*xt*yt
    dy[2] = -c*yt + gamma*xt*yt
    list(dy)
}


# Parts (b):
#
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] # the numerical function x(t)
ys = res[, 3] # the numerical function y(t)

T = my_period(ts, xs-mean(xs))
T_linear = 2*pi/sqrt(a*c)
print(sprintf('T_linear= %f; T= %f', T_linear, T))

my_trapezoidal = function(ts, ys){
    ##
    ## A poor i.e. quick numerical implementation of a quadrature routine
    ##
    n = length(ys)

    delta_t = mean(diff(ts)) # expect the "x" grid to be uniform
    area = ( ys[1] + 2 * sum(ys[2:(n-1)]) + ys[n] ) * delta_t / 2
    return(area)
}

mask = ts<T
x_bar = my_trapezoidal(ts[mask], xs[mask])/T
y_bar = my_trapezoidal(ts[mask], ys[mask])/T
print(sprintf('x_bar_linear= %f; x_bar= %f', x_bar_linear, x_bar) )
print(sprintf('y_bar_linear= %f; y_bar= %f', y_bar_linear, y_bar))

# Part(c):
#
x0 = c/gamma # the critical point
y0 = a/alpha

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

As = seq(1.01, 2, length.out=10) # move the initial conditions(x, y) from the critical point(x0, y0) by multiplying by this factor
x_bar_estimates = c()
y_bar_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)
    ys = res[, 3] # the numerical function y(t)
    
    x_bar = my_trapezoidal(ts[mask], xs[mask])/T
    y_bar = my_trapezoidal(ts[mask], ys[mask])/T

    x_bar_estimates = c(x_bar_estimates, x_bar)
    y_bar_estimates = c(y_bar_estimates, y_bar)
}

par(mfrow=c(1, 2))

plot(As, x_bar_estimates/x_bar_linear, type='b', lty=1, lwd=2, pch=19,
     xlab='expansion factor from critical point',
     ylab='x_bar/x_bar_linear',
     main='x estimate over linear estimate')
grid()

plot(As, y_bar_estimates/y_bar_linear, type='b', lty=1, lwd=2, pch=19,
     xlab='expansion factor from critical point',
     ylab='y_bar/y_bar_linear', main='y estimate over linear estimate')
grid()

par(mfrow=c(1, 1))

#dev.off()