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 - xt^2 ) * yt - xt 
    list(dy)
}

# Part (a):
#
ts = seq(0, 30, length.out=1000)

mus_requested = c( 0.2, 1.0, 5.0 ) # the special values of mu we want to compute 
for( m in mus_requested ){
    diff_eq_params = list(mu=m)
    res = ode(y=c(-3, 2), times=ts, fun=my_yprime, parms=diff_eq_params) # guess the initial conditions the book uses for its plot
    us = res[, 2] # the numerical function u(t)

    plot(ts, us, type='l', xlab='t', ylab='u(t)')
    grid()

    A = my_amplitude(us)
    T = my_period(ts, us)
    print(sprintf('mu= %f; Amplitude (est)= %f; Period (est)= %f', diff_eq_params$mu, A, T))
}

# Part (b-c):
#
mus = seq(0.2, 6.0, length.out=50)
u_period_estimates = c()
for( m in c(mus_requested, mus) ){ # put the special samples "at the front"
    diff_eq_params = list(mu=m)
    res = ode(y=c(-3, 2), times=ts, fun=my_yprime, parms=diff_eq_params)
    us = res[, 2] # the numerical function u(t)
    
    T = my_period(ts, us)
    
    u_period_estimates = c(u_period_estimates, T)
}

## Extract the special period 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_14_plot.eps', onefile=FALSE, horizontal=FALSE)

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()

##dev.off()