MATLAB: Help!! I have no idea what this question is asking, nor how to solve it…

projectile motion

Exercise 4: 2D Model with Air Resistance-Naval Gun
Model a shell red from one of the 16 inch guns on a World War II era battleship. Start with your realistic (and accurate!) 2D projectile model with air resistance, and then use real parameter values for a round red from the naval gun. Since a range table exists (<http://www.navweaps.com/Weapons/WNUS_1650_mk7.php>) for the Mark 7, let's model this gun ring an AP Mark 8 shell, so we have some data to which we can compare the results of our model. Of particular relevance for your projectile model are the Mark 7 muzzle velocity (magnitude of the initial velocity vector), the mass of the shell, and the radius of the shell (important for calculating cross-sectional area A). There is one addition that must be made in order to make your battleship salvo model realistic. Because your projectile, as you will soon discover, is capable of attaining a relatively high altitude you will have take 1 into account the fact that the atmospheric density varies with the vertical position of the red shell. To model a realistic projectile that has been red with a WWII era naval gun, we need to build into our model an altitude-dependent air density. The following equation accurately describes how the atmospheric density decreases with altitude: %(y) = %01− cy T0α, (1) where %0 = 1.2 kg/m3, c = 6.5×10−3 K/m, T0 = 300 K (sea level temperature in Kelvins), and α= 2.5. This equation is derived from thermodynamic considerations, and is a good approximation for the behavior of the earth's atmosphere up to about 10 kilometers. Note that the density depends only on the vertical distance above the surface of the Earth. The most dicult parameter to determine is the drag coecient D. The typical shells red from the big naval guns were cylindrical with a point at the front end. Furthermore, the naval guns were ried, so that the cylindrical projectiles were spinning during their ight. Thus, the drag coecient D should be expected to be a very complicated function of velocity. However, we can determine an eectively constant value via the following procedure. The drag coecient in your working program can be varied until the range produced in your model for dierent angles closely agrees with the actual measured values found in the Range Table (for new gun muzzle velocities) at the aforementioned web site. Reproduce the table of Range values, add another column for Drag Coecient, and then for each angle specied, vary D in your model until the range is matched and record the value of D in the table. Include a column in your table for the maximum altitude achieved for each of the elevation angles.The Range Table is reproduced below for your convenience.
So far my code is:
n=350; %Number of Plots 
m=0.0585; %mass of object (kg)
d=0.06858;%diameter of object 
dt=0.01; %time increment (s)
A=pi*(d/2)^2;%cross-sectional area 
p=1.225;%density of air (kg/m^3)
D=0.5;%drag coeff. 
g=9.81;%gravity force (m/s^2)
v_0=30; %initial velocity 
theta=45; %launch angle (degrees)
time=zeros(1,n);
velocity=zeros(1,n); 
velocity(1)=v_0;
x_position=zeros(1,n); %Array of x position
x_position(1)=0;
y_position=zeros(1,n); %array of y position
y_position(1)=0;
x_velocity=zeros(1,n); %array of x velocity 
x_velocity(1)=v_0*cosd(theta);
y_velocity=zeros(1,n); %array of y velocity
y_velocity(1)=v_0*sind(theta);
x_acceleration=zeros(1,n); %array of x acceleration 
y_acceleration=zeros(1,n); %array of y acceleration 
for i=2:n
  time(i)=time(i-1)+dt;
  x_acceleration(i)=-((D*p*A)/(2*m))*velocity(i-1)*x_velocity(i-1);
  y_acceleration(i)=-((D*p*A)/(2*m))*velocity(i-1)*y_velocity(i-1)-g; 
  x_velocity(i)=x_velocity(i-1)+x_acceleration(i-1)*dt;
  y_velocity(i)=y_velocity(i-1)+y_acceleration(i-1)*dt;
  velocity(i)=sqrt((x_velocity(i-1))^2+(y_velocity(i))^2); %equation of v
  x_position(i)=x_position(i-1)+x_velocity(i-1)*dt;
  y_position(i)=y_position(i-1)+y_velocity(i-1)*dt;
end
plot(x_position,y_position);
xlabel('x position');
ylabel('y position');
title('Trajectory of a Tennis Ball');
But I don't know how to do anything else…

Best Answer

You have a couple errors in your code.
First, your calculation of the area is wrong. Area is pi*radius^2, or pi*diameter^2/4. Next, you have not computed the drag force correctly. The drag force depends on the local air density, which is a function of altitude. Unfortunately, the problem statement seems to give a function for the air density that does not work. So I recommend that you use the following:
K = 2.256e-5;     % curve fit coefficient
rhomsl = 1.225    % air density at mean sea level (Kg/m^3)
rho = rhomsl*(1-K*H)^4.2561  % this is the local air density
where H is the altitude above mean sea level (not above the ground). H is computed from
H = y + ygl    % y is your trajectory vertical coordinate (above the ground)
               % and ygl is the altitude of the ground plane above sea level (meters) (a constant)
The above must be added inside your loop. Now also add the drag force calculation:
Drag = A/2*(rho*speed^2)*(drag_coefficient)  % this is the drag force (N)
This force acts to oppose the velocity vector, so you have to compte the direction of the velocity. The flight path angle, gamma, is the angle of the velocity vector with the ground:
gamma = atan(y_velocity/x_velocity);
x_force = -Drag*cos(gamma);
y_force = -Drag*sin(gamma) - Weight;  % Weight  is mass*g
Note that the forces acts in the negative direction. Now the accelerations are computed from;
x_acceleration = x_force/mass;
y_acceleration = y_force/mass;
Finally, I suggest that you let your calculation loop run until the y position goes negative (below the ground). (You can't know in advance how many caculations this will require)
Putting this all together, you have:
H = y_position + ygl;    
rho = rhomsl*(1-K*H)^4.2561  % this is the local air density 
Drag = A/2*(rho*speed^2)*(drag_coefficient);  % Drag force  
gamma = atan(y_velocity/x_velocity);
x_force = -Drag*cos(gamma);
y_force = -Drag*sin(gamma) - Weight;
x_acceleration = x_force/mass;
y_acceleration = y_force/mass;
x_velocity = x_velocity(i-1) + x_acceleration*dt;
y_velocity = y_velocity(i-1) + y_acceleration*dt;
speed = sqrt(x_velocity^2 + y_velocity^2);
x_position = x_position(i-1) + x_velocity*dt;
y_position = y_position(i-1) + y_velocity*dt;
(Of course, I have left out many of the subscripts for brevity). You can set the ground elevation (ygl) to the desired altitude above sea level and see how this effects your trajectory. (Note that the curve fit equation for density that I gave you is only good for the troposphere, which is from sea level up to about 11,000 meters.
Related Question