Partnership for Integration of Computation into Undergraduate Physics

##Exercise 1--Set up and some math## In the simplified scenario, the origin of the *xy* axes is where the asteroid would hit the Earth if there were no intervention. The Earth is on the +*x* axis with initial position $x_{E}$, moving in the negative *x* direction with speed $v_{E}$. The asteroid is on the +*y* axis with initial position $y_{a}$, moving in the negative *y* direction with speed $v_{a}$. The initial separation between the Earth and the asteroid is *d*. Draw a diagram to depict the initial situation, including the axes as well as arrows that show the initial velocities. To ensure that a collision occurs absent intervention, the values of the initial positions, the initial separation, and the initial speeds cannot all be independent. We will assume the particular initial separation *d* (which later will be how far from the Earth the intervention occurs) and the particular initial speeds $v_{E}$ and $v_{a}$. From those quantities, derive equations for $x_{E}$ and $y_{a}$ . Hint: The Earth and the asteroid must both arrive at the origin at the same time *t*. ##Exercise 2--Preliminary computations## The code below shows how to set up "variables" for the speed **vE** of the Earth (orbital speed) and the velocity vector **velEarth** of the Earth. (On the left of **=** is the name of a variable (a space in computer memory) and on the right is the value to put there.) In vpython, **vec(x,y,z)** creates a vector object with x-component **x**, y-component **y**, and z-component **z**. Also, since vpython does not understand units of physical quantities, you should write all numbers assuming that they have SI units, so that the calculations come out with consistent SI units. ``` vE = 2.978e4 velEarth = vec(-vE,0,0) ``` This choice of **velEarth** fulfills the requirements given in Exercise 1, that the Earth moves in the negative *x* direction with speed $v_E$ . Type the code into vpython. Then add code to create the speed **vA** and the velocity vector **velAst** for the asteroid--you can look up typical numbers for the asteroid speed, and remember that in this scenario the asteroid's motion is in the -*y* direction. Finally, set up a variable **dmoon** for the distance from the Earth to the Moon (to help for a sense of scale) and a variable **d** for the initial distance between the Earth and the asteroid (later to be the time of intervention). For instance, you can set *d* to equal thirty times the distance from the Earth to the Moon. In python, that calculation is **30\*dMoon**. The code given below shows how to do the calculation of **x**, the initial *x* coordinate of the Earth, denoted $x_{E}$ in Exercise 1. (In python **sqrt(number)** gives the square root, and **number\*\*2** gives the square of the number.) The code then sets up the vector **posEarth** for the position of the Earth. Finally, the code prints the velocity and position vectors of the Earth. The **print()** "function" will print the quantities inside the parentheses (the "arguments" of the function, where commas separate the different arguments of a function.) An argument consisting of text within quotation marks will print the text as it appears, while an argument consisting of a variable name not in quotation marks means that the value of the variable gets printed. Add this code to the previous code. Then add lines to calculate **y**, the initial *y* coordinate of the asteroid, denoted by $y_{a}$ in Exercise 1, and the position vector **posAsteroid** of the asteroid. Finally, add lines to print the position and velocity of the asteroid, and run the program to check that the numbers are correct. ``` x = d / sqrt(1+(vA/vE)**2) posEarth = vec(x,0,0) print("Earth velocity", velEarth, "m/s") print("Earth position", posEarth, "m") ``` ##Exercise 3--Changing positions## The code below first sets up the variable **closest** that will track the closest distance between the asteroid and the Earth as time passes, by making **closest** equal to the initial distance. The code then sets up the initial time **t** and the time step **dt**. The code then does a loop of iterations of updating the position of the Earth, assuming that the velocity remains constant. The line **while...:** begins the loop and specifies that it will run as long as the condition after **while** is true. In this case, the calculations continue until the asteroid crosses the orbit of the Moon on the other side of the Earth. The indented lines (type four spaces or a tab at the beginning of each indented line) are the lines that repeat until the loop ends. After the loop ends, the program continues with the unindented lines that come after the loop. If you are not familiar with computer coding, lines such as **t = t + dt** may seem strange and mathematically ridiculous. The program will calculate the expression on the right side using the existing values of variables, then set the variable on the left side to be the result of the calculation. In this case, the program will take the original value of **t**, add **dt** to it, then place the result of that calculation back into the memory space labelled **t**, thus replacing the value of **t** with a new value to update time. In the same way, **posEarth** gets updated by adding **velEarth** multiplied by **dt**. This code is equivalent to the mathematical expression $$\vec{R_{Ef}} \ =\ \vec{R_{Ei}} \ +\ \vec{v_E}*\Delta t . $$ The difference is that one variable, **posEarth**, is used to store each value of the position vector as it changes instead of giving the position vector a new name after each update calculation. For constant velocity, the computation could be done all at once, without the loop. Using the loop structure here is to give some familiarity with the later structure needed when velocity is changing. Type in the code. Add lines (immediately after the line beginning with **\#**, which denotes a comment not code) to update the asteroid position (assuming constant velocity). Within the loop the variable **dist** is recalculated as the magnitude (**mag**) of the relative position vector between the Earth and the asteroid (**relPos**). Then the variable **closest** is updated using the **min()** function, which "returns" the minimum of the two variables within the parentheses. Run the program to confirm that, absent any intervention, the asteroid will hit the Earth. ``` closest = d t = 0 dt = 0.1 while posAsteroid.y > -dMoon: t = t + dt posEarth = posEarth + velEarth*dt # add indented code here to update posAsteroid relPos = posAsteroid - posEarth dist = mag(relPos) closest = min(closest, dist) print("closest approach", closest, "m") ``` ##Exercise 4--Attaching a rocket## In this scenario, a rocket is attached to the asteroid to apply a constant force. Before the loop in the code from Exercise 3 (preferably near the beginning of the program), set up one variable for the mass of the asteroid (you can perhaps look up the mass of an asteroid whose impact would cause destruction on a large scale), a variable for the magnitude of the force (the rocket thrust, look up a typical number to use), and another variable for the force vector. You may want to review vector components and find out how to do the needed trig calculations in vpython. You will now need to add a line inside the loop to use the applied force to update the velocity of the asteroid. This line should be right before the line that calculates the position of the asteroid, similar to the second line in this code snippet: ``` while t <= tfinal: vel = vel + (force/mass)*dt pos = pos + vel*dt t = t + dt ``` The method of first calculating velocity then calculating position is known in physics as the Euler-Cromer method (other names are also used). If the position is calculated before the velocity, the calculations are mathematically unstable (small rounding errors get exponentially bigger). Add the needed line of code inside the loop from Exercise 3. Run the program. Is the Earth saved? If not, experiment with changing the thrust of the rocket and/or its direction, as well as changing how far from the Earth the asteroid is when the rocket starts pushing on it (labelled **d** in the previous exercise). Give conceptual discussion of the results. ##Exercise 5--Theory of hitting the asteroid## Saving the Earth in this scenario will happen by hitting the asteroid, for instance with a slug or spaceship that will stick to the asteroid. Assuming that the velocity and mass of the asteroid and the velocity and mass of the slug before the impact are known, derive the equation for the velocity vector of the asteroid (with slug imbedded) after the slug hits it. ##Exercise 6--Computations for hitting the asteroid## For simplicity, we will assume that the slug hits the asteroid at the initial position (when the asteroid is distance **d** from the Earth). The velocity in the update loop, therefore, is the final velocity from conservation of momentum. You will need to add lines before the loop in the code from Exercise 3 to set up the mass of the asteroid as well as the velocity vector and the mass of the slug before the hit. You will then need to add a line to calculate the velocity vector of the asteroid after the hit. You may want to review vector components and find out how to do the needed trig calculations in vpython. Inside the loop, no change is needed. The velocity of the asteroid after being hit, calculated before the loop, is assumed constant for the remainder of this simulation. Run the program. Is the Earth saved? If not, experiment with the mass of the slug, its velocity (try different speeds and directions), and the distance of the asteroid from the Earth when the asteroid is hit to find a way to save the Earth. For the same slug speed, which direction of its velocity is more effective? Can you give a conceptual explanation for the result? ##Exercise 7--Effect of gravity without intervention## The code below sets out a structure for including gravity. There will be two loops--one to run time backwards starting from the collision of the asteroid with the Earth and ending earlier, to establish initial conditions that would lead to the collision (doing this mathematically would be harder than in the earlier exercises). The second loop will then run time forward from the initial conditions to confirm that with those initial conditions, absent intervention the asteroid would collide with the Earth (within a margin of error). In this exercise, you will fill in missing lines of code. In the Exercise 8, you will add an intervention (asteroid hit by slug or asteroid pushed by rocket) to save the Earth. The first two lines of code show how to set up a graph and include a plot of the *x* coordinate of Earth vs time. You will add code to create other plots to help visualize that the asteroid will in fact hit the Earth without intervention and, in Exercise 8, how an intervention will change the situation. In this code, gravity acts on the asteroid. To simplify the simulation, the motions of the Earth and the Moon are simply calculated as functions of circular orbits instead of including gravity acting on them. The parts of this exercise are imbedded within the code, which contains several comments **# exercise** that are points where you need to add code, using what you've seen in the previous exercises and what is in the sample code below. Add the missing code. Run the program to confirm that, given the initial conditions found in the first loop, the asteroid will hit the Earth in the second loop within a margin of error. Comment on the accuracy of the calculation, which may be important for Exercise 8. ``` # set up the graph and add individual plots g1 = graph(ymax=3e7,ymin=-3e7,xmin=-0.0005e6,xmax=0.0005e6, title="euler-cromer dt=1 s, tfinal=2e6 s") # setting up a plot for x-component of Earth position gEarthx = gcurve(color=color.green,label="Earth x") # exercise # set up plots for Earth y, asteroid x and y, and the # distance between the Earth and the asteroid # exercise # create, and initialize or calculate, variables for # Earth-Sun distance, position of Sun (Earth at origin at collision), # angular frequency of Earth orbit, radius of Earth, position # of asteroid at collision (surface of Earth), Earth-Moon # distance, angular frequency of Moon orbit around Earth, # G, and masses of asteroid, Earth, Moon, and Sun # set time 0 at collision, run backwards t = 0 dt = -10 # choose amount of time to run backwards tf = 2e6 # set up the loop to run time backwards # use Euler-Cromer method while t > -(tf): # Motion of Earth # use function, assume circular orbit about sun # set up to be at origin at time zero posEarth = posSun + dSun*vec(-sin(omegaE*t), cos(omegaE*t),0) # exercise # do the same for the Moon # set up for "half-moon" at time zero posMoon = # function of time, Moon orbits Earth position # motion of asteroid, using gravity from Sun, Moon, Earth # calculate individual fields, then add up for acceleration # effect of Sun Sun2Ast = posAst - posSun gSunAtAst = -G*mSun*norm(Sun2Ast)/(mag(Sun2Ast))**2 # exercise # do the same for the Earth and the Moon Earth2Ast = posAst - posEarth gEarthAtAst = -G*mEarth*norm(Earth2Ast)/(mag(Earth2Ast))**2 Moon2Ast = posAst - posMoon gMoonAtAst = -G*mMoon*norm(Moon2Ast)/(mag(Moon2Ast))**2 # total acceleration of the asteroid aAst = gSunAtAst + gEarthAtAst + gMoonAtAst # exercise # update velocity and position of asteroid, and update time velAst = velAst + aAst * dt posAst = posAst + velAst*dt t = t + dt # print the result, which will be the initial # conditions for asteroid approaching Earth print("At time", t, "s") print("Asteroid position", posAst, "m \n and velocity", velAst, "m/s") print("Earth position", posEarth, "m") print("Moon position", posMoon, "m") # restart time, to run forward # t = 0 dt = -dt # create and initialize variable to track minimum # Earth-asteroid distance dist = mag(posAst-posEarth) # set up loop to run for time a bit bigger than tf # update velocities and positions # update minimum Earth-asteroid distance while t <= 0.01*tf: # exercise # write the needed code for Earth, Moon, asteroid motions t = t + dt # calculate distance and track minimum distance dist = min(dist, mag(Earth2Ast)) # this is within the loop # if time is in a small range of expected time of # collision, track x and y coordinates # sample to graph x-coordinate of Earth vs time gEarthx.plot(t,posEarth.x) # exercise # do the same for the other coordinates and for the # distance between the asteroid and the Earth print("min distance on forward run", dist, "m") ``` ##Exercise 8--Effect of gravity with interventions## Adapt the program from Exercise 7 by adding either something hitting the asteroid (see Exercise 6) or a rocket attached to the asteroid (Exercise 4). Adjust the parameters of the intervention to save the Earth.