g = 32
h = 3
r = 1/5

L = 350 # how far is the fence
H = 10 # how heigh is the fence

u = 125 # initial velocity

#u = 175 # initial velocity (play with some velocities until we find some angles that get the ball over the fence)

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

t = seq( 0, 8, length.out=500 )

A = 30 * (pi/180)
x_1 = (u*cos(A)/r)*( 1 - exp(-r*t) )
y_1 = (u*sin(A)/r)*( 1 - exp(-r*t) ) + (g/r^2)*( 1 - exp(-r*t) ) - g*t/r + h

A = 45 * (pi/180)
x_2 = (u*cos(A)/r)*( 1 - exp(-r*t) ) 
y_2 = (u*sin(A)/r)*( 1 - exp(-r*t) ) + (g/r^2)*( 1 - exp(-r*t) ) - g*t/r + h

A = 60 * (pi/180)
x_3 = (u*cos(A)/r)*( 1 - exp(-r*t) )
y_3 = (u*sin(A)/r)*( 1 - exp(-r*t) ) + (g/r^2)*( 1 - exp(-r*t) ) - g*t/r + h 

min_x = min( c( x_1, x_2, x_3 ) )
max_x = max( c( x_1, x_2, x_3 ) )
min_y = min( c( y_1, y_2, y_3 ) )
max_y = max( c( y_1, y_2, y_3 ) )

min_x = 0
max_x = 600
min_y = -200
max_y = 200

plot( x_1, y_1, col='blue', xlim=c(min_x, max_x), ylim=c(min_y, max_y), type='l', xlab='x', ylab='y', main='Problem 31' )
lines( x_2, y_2, col='green', type='l', xlab='x', ylab='y', main='Problem 31' )
lines( x_3, y_3, col='red', type='l', xlab='x', ylab='y', main='Problem 31' )
#abline(v=L, col='black')
#abline(h=H, col='black')
legend( 'bottomleft', c('A=30 deg', 'A=45 deg', 'A=60 deg'), lwd=2, lty=1, col=c('blue', 'green', 'red') )
grid()

#dev.off()



# Part (d):
#
# For a fixed velocity plot the value of y as a function of angle A at x=L
#
u = 175 # for this value of u there exist many angles where the ball clears the wall

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

# We decrease the value of u until we get a very small region where the ball will clear the wall
#
u = 145.3 # at this point the region of valid angles has shrunk to a point

as = seq( 0, pi/3, length.out=100 )
x = L
ys = x * tan(as) + ( g*x ) / ( r*u*cos(as) ) + ( g / r^2 ) * log( 1 - ( r*x )/( u*cos(as) ) ) + h
plot( as * (180/pi), ys, type='l', xlab='A (degrees)', ylab='y(A)' )
abline(h=H, col='red')
grid()

#dev.off()