g = 32
u = 125
h = 3

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

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

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

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

A = 60 * (pi/180) 
x_3 = u*cos(A)*t
y_3 = -(g*t^2)/2 + u*sin(A)*t + 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 30' )
lines( x_2, y_2, col='green', type='l', xlab='x', ylab='y', main='Problem 30' )
lines( x_3, y_3, col='red', type='l', xlab='x', ylab='y', main='Problem 30' )
legend( 'bottomleft', c('A=30 deg', 'A=45 deg', 'A=60 deg'), lwd=2, lty=1, col=c('blue', 'green', 'red') )
grid()
 
#dev.off()


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

g = 32
L = 350
H = 10
u = 110

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

#dev.off()

library(stats)
f = function(a){
    -(g/2) * ( L / (u*cos(a)) )^2 + L * tan(a) + h - H
}
res = uniroot( f, c(30, 40)*(pi/180) )
print(res)
res = uniroot( f, c(50, 60)*(pi/180) )
print(res)


# Part (f):
#
# Method 1: Decrease the value of u until we get a very small region where the ball will clear the wall
# 
g = 32
L = 350
H = 10

u = 106.9

as = seq( 0, pi/3, length.out=100 )
Hs = -(g/2) * ( L / (u*cos(as)) )^2 + L * tan(as) + h
plot( as * (180/pi), Hs, type='l', xlab='A (degrees)', ylab='LHS', main=sprintf('u= %f', u) )
abline(h=H, col='red')
grid()


# Method 2: Put u and A on different sides of the equation
# 
as = seq( 0, pi/3, length.out=1000 )
rhs = ( 2*cos(as)^2 / ( g * L^2 ) ) * ( L * tan(as) + h - H ) 
plot( as, rhs, type='l', xlab='A (radian)', ylab='RHS', col='black' )
grid()

am = which.max(rhs)
print( c( as[am], rhs[am], 1/sqrt(rhs[am]) ) )