#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Feb 20 14:04:31 2020 @author: Prof Roy """ import numpy as np from matplotlib import pyplot as plt import sys ############ ASSIGNMENT #1 ####################### # Question 1: We use the theory of q1: for each value of (xA,v0,h) we compute # the values of angle alpha which guarantee that A and B are connected for # given values of h. # Since 3 values of v0 are given, we set up an array for v0. # We add a fourth value which corresponds to v0min, below which there # is no solution. The fifth value is the optimal path. # On the graph, we plot 4 trajectories and the rectangular obstacle. # set the parameter values g = 9.8 # [m/s^2] xA = -50 # [m] h = 33 # [m] d = 10 # [m] epsilon= 0.0001 # a small parameter color = np.array(['b', 'r', 'o' , 'g', 'y']) # color scheme for each v0 value ymax = 90 # [m] v0min= np.sqrt(g*h*(1+np.sqrt(1.0+xA**2/h**2))) v0 = np.array([40,35,30,v0min]) # [m/s] print('v0min=',v0min,'[m/s]') fig1 = plt.figure(1) fig1.set_size_inches(10, 6) ax = plt.axes(xlim=(xA, 0), ylim=(0, ymax)) Npt = 100 # number to point to define trajectory plt.xlabel('x [m]') plt.ylabel('y [m]') plt.title('Projectile Motion') rect = plt.Rectangle((-d,0), d, h, facecolor= 'red', edgecolor= 'black') ax.add_patch(rect) def Projectile(X, v, alpha): # returns an array of [x, y] values from t=0 to tfinal # starting at A(X,0) with initial speed v and launch angle # alpha t0 = 0 tfinal = -X/(v*np.cos(alpha)) # final time to reach point B(0,h) starting from # point A(xA,0) t = np.linspace(t0, tfinal, Npt) x = X + v*t*np.cos(alpha) y = v*t*np.sin(alpha)-0.5*g*t**2 return [x, y] k=-1 for v in v0: k += 1 A= 2*h + g*xA**2/v**2 # A,B, C are 3 coefficients of quadratic equation B= 2*xA # to solve for the launch angles alpha. C= g*xA**2/v**2 discrim = B**2 -4*A*C # what follows depends on the value of discrim: # you will get 2 solutions, 1 solution or no solution. # We also test whether the path avoids the rectangular obstacle. # This is useful for question 2. if discrim < -epsilon: print('v = ',v, '[m/s], No solution!') if discrim > epsilon: roots= np.array([0.5*(-B+np.sqrt(discrim))/C,0.5*(-B-np.sqrt(discrim))/C]) for i in [0,1]: # two possible launch angles alpha = np.arctan(roots[i]) t_hit= -(d+xA)/np.cos(alpha)/v y_hit= -0.5*g*t_hit**2 + v*t_hit*np.sin(alpha) print('v = ',v, '[m/s], alpha = ',alpha*180/np.pi,'degrees','yhit= ',y_hit) if y_hit < h: print('HIT!') traj = Projectile(xA, v, alpha) plt.plot(traj[0],traj[1], color=color[k], label="v0[m/s]= {}".format(v) ) if np.abs(discrim) < epsilon: # this case is to address the case when discrim # is nearly zero: in that case, there is a double # root alpha= np.arctan(-0.5*B/C) t_hit= -(d+xA)/np.cos(alpha)/v y_hit= -0.5*g*t_hit**2 + v*t_hit*np.sin(alpha) print('v = ',v, '[m/s], alpha = ',alpha*180/np.pi,'degrees','yhit= ',y_hit) if y_hit < h: print('HIT!') traj = Projectile(xA, v, alpha) plt.plot(traj[0],traj[1], color=color[k], label="v0[m/s]= {00}".format(v)) ############ QUESTION 2 ############################### print('Part 2: Case 1, rectangular obstacle') ########### CASE 1: RECTANGULAR BUILDING ################### # # Q.1 shows that, for a given distance xA, there is a minimum # launch speed v0min necessary to reach point B. For v0< v0min, there is no solution. # We can find v0min= sqrt(g*h*(1+sqrt(1.0+xA^2/h^2)) # # We can use a scheme based on two nested loops: # First we create a loop of xA values, we increase xA as the loop proceeds. # # Inside this loop, we create an array of v0 values from v0min to v0max. # # We then create a v0 loop, by decreasing v0 as this second loop proceeds: # Inside this loop, we determine the path connecting A(xA,0) to B(0,h) by finding # the two launch angles. This path is either acceptable or not. # A path is acceptable is no point P(x,y) t=0..tfinal is inside the region defined by # the building: D< x <0 0< y h, the # path is acceptable. # # If the path is acceptable, we test if the corresponding v0 value is lower that our previously # recorded minimum. If it is, this is our new minimum. # # Eventually, as you proceed along these 2 loops, there are no longer smaller values of v0. # You can get closer to the obstacle, but it takes a larger speed v0 to reach B. # The scheme converges very fast to the solution. You can refine your solution by adding more values to # the loops. xAmax = -30 # [m] the most distant point A to O (start with sufficiently large xAmax) xAmin = -1.1*d # [m] the closest point A to origin O h = 33 # [m] height of obstacle d= 10 # [m] 1/2 width of obstacle Nxvalues = 400 # number of xA values Nvvalues = 400 # number of v0 values X= np.linspace(xAmax,xAmin,Nxvalues) # xA array Vmin= 100; # Vmin is the sought value of optimal (minimal) value of v0 for xA in X: v0min= np.sqrt(g*h*(1+np.sqrt(1.0+xA**2/h**2))) v0max= 1.2*v0min Vvalues= np.linspace(-v0max, -v0min, Nvvalues, endpoint=False)# array of v0 values Vvalues = -Vvalues for V0 in reversed(Vvalues): # v0 loop from v0max to v0min, excluding v0min # (at v0min, discrim is a very small number, which # be <0 due to machine precision) A =2*h + g*xA**2/V0**2 B = 2*xA C= g*xA**2/V0**2 discrim=B**2-4*A*C alpha1 = np.arctan(0.5*(-B+np.sqrt(discrim))/C) # launch angle 1 t1_hit= -(d+xA)/np.cos(alpha1)/V0 # corresponding time to reach x=-d y1_hit= -0.5*g*t1_hit**2 + V0*t1_hit*np.sin(alpha1) # height at x=-d alpha2 = np.arctan(0.5*(-B-np.sqrt(discrim))/C) # launch angle 2 t2_hit= -(d+xA)/np.cos(alpha2)/V0 # time to reach x=-d y2_hit= -0.5*g*t2_hit**2 + V0*t2_hit*np.sin(alpha2) # height at x=-d if (y1_hit > h): # test if path #1 is acceptable if(V0 < Vmin): # if acceptable, test if V0 is smaller than Vmin Vmin=V0 # if smaller, then this is the new minimum Xmin=xA # record the corresponding value of xA AlphaMin= alpha1 # record the corresponding value of alpha print('Rectangle Case: xA=',xA,'[m], v0=',V0,'[m/s], alpha=',alpha1*180/np.pi,'degrees') if (y2_hit > h): # repeat for the second path if(Vmin>V0): Vmin=V0 Xmin=xA AlphaMin=alpha2 print('Rectangle Case: xA=',xA,'[m], v0=',V0,'[m/s], alpha=',alpha2*180/np.pi,'degrees') # The last printed value is the optimal path. Plot the solution on the previous plot print('This is the optimal path:') print('Rectangle Case: xA=',Xmin,'[m], v0=',Vmin,'[m/s], alpha=',AlphaMin*180/np.pi,'degrees') k += 1 traj = Projectile(Xmin, Vmin, AlphaMin) plt.plot(traj[0],traj[1], color=color[k] , label="v0[m/s]= {00}".format(Vmin)) leg = ax.legend(loc="upper left", ncol=1, shadow=False, fancybox=False) for t in leg.texts: # truncate label text to 13 characters t.set_text(t.get_text()[:13]) plt.savefig('Proj1Figure1.png') plt.show() option = int(input("Continue? enter 1: ")) if option != 1: sys.exit() ########### CASE 2: SPHERICAL BUILDING ################### # # same procedure: however to check whether a trajectory is acceptable or not, we # need to compute the entire trajectory to connect A to B, given xA, v0 and, alpha. # We write a function to check the distance of all points of the trajectory to the center # of the circular building: if x**2 + (x-h/2)**2 > (h/2)**2 for all (x,y) on the path, # then the trajectory. # As in the first case, we create two nested loops of xA values and v0 values: # For each xA and v0, we compute the two trajectories from A to B. We check if they are # acceptable. If they are, we test if the value of v0 is smaller. If so, we record v0. # # Eventually, we reach the optimal solution: the solution connecting reaching B with smallest # launch speed. print('Part 2: Case 1, circular obstacle') xA = -40 # [m] h = 33 # [m] Npt= 400 # nbr of points on trajectory Nvalues = 400 X= np.linspace(xA,-1.1*d,Nvalues) Vmin= 100; # python program to check if all # values in array DistanceArray are greater # than value using all() function def check(DistanceArray, value): return(all(x > value for x in DistanceArray)) for x in X: v0min= np.sqrt(g*h*(1+np.sqrt(1.0+x**2/h**2))) v0max= 1.2*v0min Vvalues= np.linspace(-v0max, -v0min, Nvalues, endpoint=False) Vvalues = -Vvalues for V0 in reversed(Vvalues): A =2*h + g*x**2/V0**2 B = 2*x C= g*x**2/V0**2 discrim=B**2-4*A*C alpha1 = np.arctan(0.5*(-B+np.sqrt(discrim))/C) traj1 = Projectile(x,V0,alpha1) distance1 = traj1[0]**2 + (traj1[1]-h/2)**2 check1= check(distance1, h**2/4) alpha2 = np.arctan(0.5*(-B-np.sqrt(discrim))/C) traj2 = Projectile(x,V0,alpha2) distance2 = traj2[0]**2 + (traj1[1]-h/2)**2 check2=check(distance2, h**2/4) # check that trajectory1 is acceptable if (check1 == True): if(V0 < Vmin): Vmin=V0 Xmin=x AlphaMin= alpha1 print('Sphere Case: xA=',x,'[m], v0=',V0,'[m/s], alpha=',alpha1*180/np.pi,'degrees') # check that trajectory2 is acceptable if (check2 == True): if(Vmin>V0): Vmin=V0 Xmin=x AlphaMin=alpha2 print('Sphere Case: xA=',x,'[m], v0=',V0,'[m/s], alpha=',alpha2*180/np.pi,'degrees') print('This is the optimal path:') print('Sphere Case: xA=',Xmin,'[m], v0=',Vmin,'[m/s], alpha=',AlphaMin*180/np.pi,'degrees') fig2 = plt.figure(2) ax2 = plt.axes(xlim=(-40, 0), ylim=(0, 40)) ax.set_aspect('equal') plt.xlabel('x [m]') plt.ylabel('y [m]') plt.title('Projectile Motion') circ = plt.Circle((0,h/2), h/2, facecolor= 'red', edgecolor= 'black') ax2.add_patch(circ) traj = Projectile(Xmin, Vmin, AlphaMin) plt.plot(traj[0],traj[1], "y", label="v0= {00} m/s".format(v)) plt.savefig('Proj1Figure2.png') plt.show()