g = 9.8
m = 0.15
y0 = 30
y0p = 20

t = seq( 0, 5.5, length.out=1000 )

# Position and velocity for problem 20:
#
v_20 = - g * t + y0p
y_20 = -(g/2) * t^2 + y0p * t + y0

# Position and velocity for problem 21:
#
#
t_21 = t
v_21 = -30 * m * g + ( y0p + 30 * m * g ) * exp( -t/(30*m) )
y_21 = -30 * m * g * t + 30 * m * ( y0p + 30 * m * g ) * ( 1 - exp( -t/(30*m) ) ) + y0

y_21_max_height = max(y_21) 
print(sprintf('Prob 21: maximum height %f', y_21_max_height))
mi = which.min(abs(y_21))
t_hit_ground = t_21[mi]
print( sprintf('Prob 21: t when hit ground %f', t_hit_ground) )

# Position and velocity for problem 22:
#
library(deSolve)

parameters = c( m=m, g=g )
state = c( y=y0, v=y0p )
dydx = function(t, state, parameters){
  with(as.list(c(state, parameters)),{
    # rate of change
    dy = v
    dv = -g - (sign(v)*v^2)/(1325*m)
    list(c(dy, dv))
  })
}

res = ode(y=state, times=t, fun=dydx, parms=parameters)

t_22 = res[,'time']
y_22 = res[,'y']
v_22 = res[,'v']

y_22_max_height = max(y_22)
print(sprintf('Prob 22: maximum height %f', y_22_max_height))
mi = which.min(abs(y_22))
t_hit_ground = t_22[mi]
print( sprintf('Prob 22: t when hit ground %f', t_hit_ground) )

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

par(mfrow=c(1,2))
plot( t, y_20, col='black', type='l', xlab='t', ylab='y(t)' )
lines( t_21, y_21, col='blue' )
lines( t_22, y_22, col='green' )
grid()

plot( t, v_20, type='l', xlab='t', ylab='v(t)' )
lines( t_21, v_21, col='blue' )
lines( t_22, v_22, col='green' )
grid()

legend( 'topright', c('No drag', '|v| drag', 'v^2 drag'), col=c('black', 'blue', 'green'), lty=1, lwd=2 )
par(mfrow=c(1,1))

#dev.off()