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

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] = -4*sin(xt)
    list(dy)
}

diff_eq_params = list()

## Parts (a-c):
##
##postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_3_prob_23_plots_pt_ac.eps', onefile=FALSE, horizontal=FALSE)

ts = seq(0, 10, length.out=500)

As = c(0.25, 0.5, 1.0, 1.5, 2.0) # initial conditions
Xs = c() # the numerical solutions
amplitude_estimates = c()
period_estimates = c()
for( a in As ){
    res = ode(y=c(a, 0), times=ts, fun=my_yprime, parms=diff_eq_params)
    xs = res[, 2] # the numerical function x(t)
    Xs = cbind(Xs, xs)
    amplitude_estimates = c(amplitude_estimates, my_amplitude(xs) )
    period_estimates = c(period_estimates, my_period(ts, xs))
}
colnames(Xs) = As

par(mfrow=c(1, 3))
matplot(ts, Xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)')
grid()

plot(As, amplitude_estimates, type='b', lty=1, lwd=2, pch=19, xlab='A', ylab='amplitude (est)')
grid()

plot(As, period_estimates, type='b', lty=1, lwd=2, pch=19, xlab='A', ylab='period (est)')
grid()

par(mfrow=c(1, 1))

##dev.off()

# Part (d):
#
#postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_3_prob_23_plot_pt_d.eps', onefile=FALSE, horizontal=FALSE)

a=4
res = ode(y=c(a, 0), times=ts, fun=my_yprime, parms=diff_eq_params)
xs = res[, 2] # the numerical function x(t)
matplot(ts, xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)')
grid()

#dev.off()