Project 7: From the Earth to the Moon, part II

 When there is air resistance, one must use smaller timesteps to get reasonable results. I have updated the values of the timesteps you should use in the assignment below. Please check it before doing your analysis.

Thursday, April 26, at 6:00 PM
No pseudocode. Instead, please write detailed equations showing Heun's method. Write down on paper exactly how you will go from current values of height and velocity to the values at the next timestep. Take it in two pieces. First, write down the equations you will use to get the new velocity. Then, write down the equations you will use to get the new height. Use Heun's method in each case.
Monday, April 30, at noon
The finished code and analysis

Your job in this project: calculate the height and speed of a projectile fired vertically from a cannon, just like last week. This time, however, you are take into account two complicating factors:

1. gravitational force changes with distance from the Earth
2. air resistance
You will also have to compare Euler's method to Heun's method.

The situation is similar to that of the previous week: a gigantic cannon shoots a projectile vertically at high speed. The initial conditions are

• a starting height of H = 0 meters above the ground
• a starting speed of V = 5,000 meters per second
• projectile mass M = 10 metric tons
• spherical projectile radius R = 2 meters

What are the height and velocity of the projectile as a function of time?

Case A: realistic gravity only (no air resistance)

Write a Scilab routine which calculates position and velocity as a function of time for this problem. It should look like this:

```  function cannonball_a(start_velocity, timestep, method, output_file)

where

start_velocity          is the starting velocity, in meters per second

timestep                is the size of the timestep (t1 - t0) to use
in calculations, in seconds

method                  is either the word "euler", or the word "heun";
it indicates the algorithm to use for integration

output_file             is the name of a file in which you will
write values of time, height, and speed

```

Instead of assuming that the gravitational force on the projectile is constant, use the Law of Universal Gravitation to compute the force at each iteration. Assume for this assignment that the mass of the Earth is 5.98 x 10^(24) kg, and its radius is 6.37 x 10^6 m.

In this case, you should be able to compute the TRUE velocity and TRUE energy corresponding to any height; use the conservation of energy. That means that you know what the velocity should be when the cannonball returns to the Earth's surface.

Tabulate the following properties for each choice of timestep: t = 1.0, 0.5, 0.25 seconds, 0.125 seconds.

• peak height of the trajectory
• time when the projectile hits the ground
• velocity when the projectile hits the ground
• difference between your calculated final velocity and the true final velocity, in fractional terms

Your program should be able to employ either Euler's method or Heun's method to perform this simulation, depending on the value of the third argument. Which does a better job? How does the performance of Euler's method improve when you decrease the timestep by a factor of 2? How does the performance of Heun's method improve when you decrease the timestep by a factor of 2? Be quantitative.

Case B: air resistance only (no change in gravitational force)

(Assume that the gravitational force on the ball is constant for this case).

```  function cannonball_b(start_velocity, timestep, method, output_file)

where

start_velocity          is the starting velocity, in meters per second

timestep                is the size of the timestep (t1 - t0) to use
in calculations, in seconds

method                  is either the word "euler", or the word "heun";
it indicates the algorithm to use for integration

output_file             is the name of a file in which you will
write values of time, height, and speed

```

Now, however, there is air resistance pushing against the projectile as it moves through the air. Let us adopt the following rule for calculating the force of air resistance -- it's an approximation, but not a terrible one at moderate speeds.

```                                                    2
F(air)  =  C * (density) * (area) * (speed)

where

F(air)          is the force of air resistance (Newtons)

C               is the drag coeff, C = 1.0 for a smooth ball

density         is the density of the air at current altitude
(see below)

area            is the cross-section area of the projectile

speed           is the ball's current speed (meters/sec)

```

To be realistic, we must include the change of the density of air in the Earth's atmosphere with altitude. We don't ordinarily notice it, because the scale height of the change is much larger than the altitudes we typically cover in a short time. To a decent approximation, the density of air is

```			                                  (-height/8000 m)
rho(height)  =   (1.21 kg/m^3) * e
```

It might be very useful to write a very short Scilab function which calculates the density of air given a height above the ground. That might make your main program shorter and easier to understand, since some of the details would be "hidden" in the other function.

Tabulate the following properties for each choice of timestep: t = 0.01, 0.005, 0.0025, 0.00125 seconds.

• peak height of the trajectory
• time when the projectile hits the ground
• velocity when the projectile hits the ground

This time, we don't know (well, I don't know, anyway) the true velocity when the projectile hits the ground, so we can't easily figure out if our results are correct, or how incorrect they might be. Again, use both Euler's method and Heun's method. Does either method appear to converge as the timestep decreases? Be quantitative.

Bells and Whistles

1. Case C: changing gravitational forces AND air resistance. Once again, for each choice of timestep t = 0.01, 0.005, 0.0025, 0.00125 seconds, compute the peak height, total time in the air, and final velocity. Use Heun's method only.
2. In Jules Verne's book All Around the Moon, the sequel to From the Earth to the Moon, three travellers inside the projectile (!) discover that air resistance has reduced their original velocity by a certain round factor; that is, they have lost 1/N of their initial velocity by the time they leave the atmosphere. What is that factor? Perhaps it would help to scan chapters II or IV of the book ( which you can find on Project Gutenberg, if you don't have a copy). How does it compare to the effect of air resistance in your calculations? Make the rough approximation that all the air in the Earth's atmosphere lies less than 25 km from the surface.