mu_given_gamma = function(gamma){
    sqrt(4 - gamma^2)/2
}

u_fn = function(t, gamma){
    mu = mu_given_gamma(gamma)
    exp(-(gamma*t)/2) * ( 2 * cos( mu * t ) + (gamma/mu) * sin( mu * t ) )
}

compute_tau = function(gamma, t_min=0, t_max=75, y_thresh=0.01){
    #
    # Computes the time tau where |u(t)| < y_thresh for all t > tau
    #
    mu = mu_given_gamma(gamma)

    t = seq( t_min, t_max, length.out=50000 )
    y = exp(-(gamma*t)/2) * ( 2 * cos( mu * t ) + (gamma/mu) * sin( mu * t ) )

    # Walk backwards until we get the first FALSE:
    #
    below_t = abs(y) < y_thresh
    t_index = length(y)
    while( below_t[t_index] ){
        #print(c(gamma, t_index, below_t[t_index], y[t_index], y_thresh))
        t_index = t_index-1
    }
    tau = t[t_index+1]
    list(t=t, y=y, tau=tau)
}

approx_tau_fn = function(gamma){
    #
    # An approximation to tau see Part (e)
    #
    (2/gamma)*log(400/sqrt(4-gamma^2))
}

# Part (a):
#
y_thresh = 0.01
gamma = 0.25

ct_res = compute_tau(gamma)
tau = ct_res$tau
print(sprintf('tau= %.3f', tau))
tau = approx_tau_fn(gamma)
print(sprintf('gamma= %.2f; approximate tau= %.2f', gamma, tau))

#postscript("../../WriteUp/Graphics/Chapter3/sect_7_prob_25_plots.eps", onefile=FALSE, horizontal=FALSE)

par(mfrow=c(1, 3))

ts = seq(0, 50, length.out=500)
ys = u_fn(ts, gamma)
plot( ts, ys, type='l', main=sprintf('Problem 25: gamma: %.2f; tau= %.2f', gamma, tau), xlab='t', ylab='y(t)')
abline(v=tau, col='black')
abline(h=c(-y_thresh, y_thresh), col='blue')
grid()

# Part (b):
#
gammas = seq(0.5, 1.5, length.out=100)
taus = sapply(gammas, function(x){ compute_tau(x, gamma)$tau })
plot(gammas, taus, type='b', xlab='gamma', ylab='tau',
main=sprintf('%.2f < gamma < %.2f', min(gammas), max(gammas)))
lines(gammas, approx_tau_fn(gammas), type='l', col='red')
grid()

# Part (c):
#
gammas = seq(1.5, 1.99, length.out=100)
taus = sapply(gammas, function(x){ compute_tau(x, gamma)$tau })

min_tau_indx = which.min(taus)
tau0 = taus[min_tau_indx]
gamma0 = gammas[min_tau_indx]
ts = sprintf('gamma0= %.3f; tau0= %.3f', gamma0, tau0)
print(ts)

plot(gammas, taus, type='b', xlab='gamma', ylab='tau', main=sprintf('%.2f < gamma < %.2f', min(gammas), max(gammas)))
lines(gammas, approx_tau_fn(gammas), type='l', col='red')
abline(v=gammas[min_tau_indx], col='black')
grid()

par(mfrow=c(1, 1))

#dev.off()