Partnership for Integration of Computation into Undergraduate Physics

**EXERCISE 1:** Write the function **gravAccel(h)** that takes the height above the ground as a parameter, and returns the resultant gravitational **acceleration**. Make sure to take the Earth's radius and the direction of the force into account. G = 6.674E-11 #Gravitational constant M = 5.972E24 #Mass of the Earth R_Earth = 6.371E6 #Radius of the Earth ___ **EXERCISE 2:** **a)** Write the function **drag(v, rho)** that takes velocity and air density as parameters and returns the **force due to drag**. Make sure that the force acts in the direction opposite of the velocity. ___ **EXERCISE 3:** **a)** Let's assume that the atmosphere has a uniform density of $1.225kg/m^3$ (the average density of air at sea level) all the way up. Using a gravitational acceleration of $9.81m/s^2$, a mass of $80kg$ and drag with an air density of $1.225kg/m^3$, set gravity and drag to be equal and find the resultant velocity. *(you can do this by hand with a calculator)* Your result from **a)** is the terminal velocity, the velocity where drag and gravity cancel each other out. When you reach this velocity, acceleration becomes zero and velocity becomes constant. **b)** Given your result in **a)**, would you be able to break the sound barrier by falling through an atmosphere with a constant air density of $1.225kg/m^3$? Justify your answer. ___ **EXERCISE 4:** There is a formatting error that prevents the symbol < from appearing in the code. I have added additional comments where needed. #Takes the height in meters as an argument and returns the air density at that height in kg/m^3 def airDensity(h): if h < 25000: #Less than T = -131.21 + 0.00299 * h p = 2.488 * ((T + 273.1) / 216.6)**(-11.388) elif 11000 < h <= 25000: #Less than, less than or equal to T = -56.46 p = 22.65 * np.exp(1.73 - 0.000157 * h) else: h = max(0, h) #Negative heights make no sense in this model, keep that in mind when analyzing results T = 15.04 - 0.00649 * h p = 101.29 * ((T + 273.1) / 288.08)**(5.256) #This formula was made to fit observed values of the air density. It is not a physical equation or law return(p / (0.2869 * (T + 273.1))) #The air density ___ **a)** Use the **airDensity(h)** function to find the air density at 39045 meters. Find the terminal velocity with this air density and the parameters from **exercise 3**. Would you be able to break the sound barrier by falling through air with this density? Explain your answer. **b)** What would the answer to **a)** be with the air density at 1000 meters? ___ **EXERCISE 5:** **a)** Write the function **acceleration(h, v)** that takes the height above the ground and velocity as parameters, and returns the acceleration from gravity and drag. Use the **gravAccel**, **drag** and **airDensity** functions you have defined previously. Make sure to use the acceleration from drag, and not the force. You can set the mass of the skydiver equal to $80kg$. ___ **EXERCISE 6:** #These values gives us a good look at a fall from 0 to 100km n = 5000 #The number of steps in the Euler-Cromer calculation dt = 0.05 #The size of each time-step in the Euler-Cromer calculation pos = np.zeros(n) #Height in meters. The array is length n vel = np.zeros(n) #Velocity in m/s acc = np.zeros(n) #Acceleration in m/s^2 time = np.linspace(0, n*dt, n) #This creates an evenly spaced array of time values, ranging from 0 to n*dt pos[0] = 41419 #The starting height of the fall in meters. Setting it equal to the fall height record ___ **a)** Do the Euler-Cromer calculation of a fall through a non-uniform atmosphere with the initial conditions from above, and an acceleration given by the **acceleration(h, v)** function. (The current position during the loop will in this case be pos[i]) Here is a template for how the Euler-Cromer method can look when storing values in arrays: for i in range(n-1): acc[i] = function(arguments) vel[i+1] = vel[i] + acc[i] * dt pos[i+1] = pos[i] + vel[i+1] * dt #The last value of this array needs to be calculated separately acc[n-1] = function(arguments) ___ **EXERCISE 7:** Plotting often involves a lot of functions to fill in all of the necessary information, though a simple "plot(time, pos)" can often be sufficient. Here is an example of how code to plot these results can look. Hopefully the previous calculations made these plots look reasonable plt.figure(figsize=(19,4)) #Makes room for three plots plt.subplot(131) #Placing the graphs side-by-side plt.plot(time, pos/1000) #Dividing by 1000 here to show height in kilometers plt.xlabel("Time [s]") plt.ylabel("[km]") plt.title("Height") plt.grid() plt.subplot(132) plt.plot(time, vel, color="orange") plt.xlabel("Time [s]") plt.ylabel("[m/s]") plt.title("Velocity") plt.grid() plt.subplot(133) plt.plot(time, acc, color="green") plt.xlabel("Time [s]") plt.ylabel("[m/s^2]") plt.title("Acceleration") plt.grid() plt.show() ___ **a)** If all went well, your acceleration graph should have a somewhat strange shape(it should not be just straight lines though). Justify this shape using the formula for air resistance. In particular, try to explain the value of the initial acceleration, the acceleration with the largest positive value and the acceleration near the end. **b)** What effect would increasing the weight of the skydiver have on your results? **c)** What would change about the shape of the acceleration graph if the air density was constant? ___ **EXERCISE 8:** **a)** Did the skydiver break the speed barrier (343m/s) during the fall? If so, do you think he or she could have done so with a lower initial height? ___ **EXERCISE 9:** In the previous section, you wrote code to calculate the motion of an object falling through the atmosphere. Your next task will be to turn this code into a function, **maxVel(h)**, which takes in a starting height as a parameter, runs the Euler-Cromer calculation and returns the top speed reached during the fall. You can then apply the bisection method to this function to find the input height necessary to break the sound barrier (for simplicity's sake, take the final **ceiling** as your answer). You are free to choose your own search space, but we would suggest starting with floor of 0 meters and a ceiling of 100 000 meters (known as the Karman Line, where space begins). ___ To implement the bisection method as described above, you will first need a function that helps you determine whether the sound barrier was broken during the fall or not. **a)** Write the function **maxVel(h)** that takes a starting height above the ground as an argument, runs the Euler-Cromer calculation like you did in **Exercise 5**, and returns the top speed reached during the fall. Numpy arrays have built in functions which make this easier, for example vel.min() and vel.max(). ___ **EXERCISE 10:** **a)** Implement the bisection method as described above. Start with a floor of 0m and a ceiling of 100 000m, and make it run 17 iterations (with a for loop for instance). Use the function **maxVel** from exercise 9. **b)** Print out the final value of the ceiling, which will be your answer. Also print out the final value of the floor. What is the biggest possible error your answer can have? How can you reduce it? ___ **Bonus Exercise:** **a)** Use the starting height you found in exercise 8 and calculate the position, velocity and acceleration during the fall. Make sure to not name your arrays and variables the same names as in the earlier code, as this might produce undesirable results. **b)** Copy and alter the plotting code to plot your results from **a)**