g = 32.2 w = 180 y0 = 5000 y0p = 0 # # [1] "when jump ends: t= 266.182365; y= 2.881754; v= -15.000000" # t = seq( 0, 20, length.out=500 ) # integrate to t=275 to get to the end of the sump library(deSolve) parameters = c( w=w, g=g ) state = c( y=y0, v=y0p ) dydx = function(t, state, parameters){ with(as.list(c(state, parameters)),{ # rate of change dy = v if( t<10 ){ dv = -g - ( 0.75 * g / w ) * v }else{ dv = -g - ( 12 * g / w ) * v } list(c(dy, dv)) }) } res = ode(y=state, times=t, fun=dydx, parms=parameters) ts = res[,'time'] ys = res[,'y'] vs = res[,'v'] open_index = which.min(abs(ts-10)) print( sprintf('when parachute opened: t= %f; y= %f; v= %f', ts[open_index], ys[open_index], vs[open_index]) ) print( sprintf('distance fallen= %f', y0 - ys[open_index] ) ) print( sprintf('limiting velocity= %f', vs[length(vs)]) ) end_index = which.min(abs(ys)) # when the jump ends y(t) = 0 print( sprintf('when jump ends= t= %f; y= %f; v= %f', ts[end_index], ys[end_index], vs[end_index]) ) #postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_3_prob_23_plot.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1,2)) plot( ts, ys, col='black', type='l', xlab='t', ylab='y(t)' ) grid() plot( ts, vs, type='l', xlab='t', ylab='v(t)' ) grid() par(mfrow=c(1,1)) #dev.off()