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] = yt
    dy[2] = parameters$mu * ( 1 - (1/3)*yt^2 ) * yt - xt 
    list(dy)
}

# Parts (c):
#
ts = seq(0, 30, length.out=1000)

res = ode(y=c(5, 2), times=ts, fun=my_yprime, parms=list(mu=1)) # guess the initial conditions the book uses for its plot
us = res[, 2]
plot(ts, us, type='l', xlab='t', ylab='u(t)')
grid()

A = my_amplitude(us)
T = my_period(ts, us)
print(sprintf ('gammma= %f; Amplitude (est)= %f; Period (est)= % f', 1.0, A, T))


# Part (d):
#
mus_requested = c( 0.2, 0.5, 2.0, 5.0 ) # the special values of mu we want to compute 

mus = seq(0.2, 6.0, length.out=25)
u_amplitude_estimates = c()
u_period_estimates = c()
for( m in c(mus_requested, mus) ){ # put the special samples "at the front"
    res = ode(y=c(5, 2), times=ts, fun=my_yprime, parms=list(mu=m) )
    us = res[, 2] # the numerical function u(t)
 
    A = my_amplitude(us)
    T = my_period(ts, us)

    u_amplitude_estimates = c(u_amplitude_estimates, A)
    u_period_estimates = c(u_period_estimates, T)
}

## Extract the special amplitude and period estimates:
##
u_amplitude_estimates_requested = u_amplitude_estimates[1:length(mus_requested)]
u_amplitude_estimates = u_amplitude_estimates[(length(mus_requested)+1):length(u_amplitude_estimates)]

u_period_estimates_requested = u_period_estimates[1:length(mus_requested)]
u_period_estimates = u_period_estimates[(length(mus_requested)+1):length(u_period_estimates)]

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

par(mfrow=c(1, 2))

plot(mus, u_amplitude_estimates, type='l', lty=1, lwd=2, pch=19, xlab='mu', ylab='amplitude estimate')
lines(mus_requested, u_amplitude_estimates_requested, type='p', pch=19, col='blue')
grid()

plot(mus, u_period_estimates, type='l', lty=1, lwd=2, pch=19, xlab='mu', ylab='period estimate')
lines(mus_requested, u_period_estimates_requested, type='p', pch=19, col='blue')
grid()

par(mfrow=c(1, 1))

##dev.off()