#postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_3_prob_18_plot.eps", onefile=FALSE, horizontal=FALSE) 

omega = pi/12
T0 = 60
T1 = 15
k = 0.2

ts = seq( 0, 75, length.out=500 )
T = T0 + T1 * cos( omega * ts )
pt1 = ( omega * sin( omega * ts ) ) / ( k^2 + omega^2 )
pt2 = ( k * cos( omega * ts ) ) / ( k^2 + omega^2 )
S = T0 + k * T1 * ( pt1 + pt2 )

plot( ts, T, col='black', type='l', xlab='t', ylab='' )
lines( ts, S, col='blue', type='l' )
abline(h=T0, col='black')
legend( 'bottomright', c('T(t)', 'S(t)'), lty=1, lwd=2, col=c ('black', 'blue') )
grid()

#dev.off()

print( max( S ) )

# Estimate the time lag between the maximums:
#
time_mask = (ts>=20) & (ts<=35)
time_sub = ts[time_mask] # subsets
T_sub = T[time_mask]
S_sub = S[time_mask]
T_max_index = which.max( T_sub )
print( c( T_max_index, time_sub[T_max_index], T_sub[T_max_index]) ) # location and maximum in T(t)
S_max_index = which.max( S_sub )
print( c( S_max_index, time_sub[S_max_index], S_sub[S_max_index] ) ) # location and maximum in S(t)

tau = time_sub[S_max_index] - time_sub[T_max_index]
print( tau )