Hopefully you completed the Slingshot with Gravity! coding activity. (If you completed the Escape Velocity / Newtonian Black Holes activity that's great, but it's NOT a requirement for working on this activity.) The Slingshot with Gravity! activity explores the behavior of objects like planets or comets interacting with a massive object like the sun. An interesting part of that activity is showing that the gravitational force between two objects is inversely proportional to the distance squared. We found that this must be true: $$F_{grav} = \frac{G M m }{r^2}$$
and we have this formula instead of something proportional to $1/r^n$ because the planets are observed to have elliptical orbits and as the coding activity shows you only get elliptical orbits if $n = 2$.
In THIS activity we are going to keep using this formula but instead of keeping one object fixed (the sun) and letting the other move around it we are going to let BOTH objects move and see what happens.
Letting both objects move is important, because many of the extra-solar planets (a.k.a. exoplanets) that we know about we only discovered because of the motion of the star that they are orbiting! So we will talk about the little white circle as a planet orbiting the sun instead of being something much smaller, like a comet or a satellite.
Click here to open up a modified version of the Slingshot with Gravity code
Here is a link to just run the code if that's all you want to do
Our goal is to let the sun drift towards the top right. Look at the code:
function draw(){
// Update velocities
vx += deltaVx;
vy += deltaVy;
// Update location
x += vx*dt;
y += vy*dt;
// velocity is unchanged if there are no forces
deltaVx = 0;
deltaVy = 0;
r = sqrt((x - x_massive)*(x-x_massive) + (y - y_massive)*(y - y_massive));
Fgrav = G*M*mass/(r*r);
theta = atan2(y - y_massive,x-x_massive);
Fx = -Fgrav*cos(theta);
deltaVx = (Fx/mass)*dt; // = ax*dt;
Fy = -Fgrav*sin(theta);
deltaVy = (Fy/mass)*dt; // = ay*dt;
// This will clear the screen and re-draw it
display();
// Add more graphics here before the end of draw()
drawSun(x_massive,y_massive,0,0,0,0);
drawBlob(x,y,vx,vy,Fx,Fy);
} // end draw()
Notice that we are updating the position of the planet but we are not updating the position of the sun!
Step 1a. Add vx_massive and vy_massive to the beginning of the code
For the sun to move, it needs a velocity in the x direction (vx_massive) and a velocity in the y direction (vy_massive). Add vx_massive and vy_massive to the beginning of the code and set them both to 10.0. Here is some code you can copy if it helps:
vx_massive = 10.0; vy_massive = 10.0;
Note that we are calling this "massive" because it is the velocity of the massive object. Later the object will be a pulsar instead of the Sun. Pulsars and stars are both pretty massive compared to planets so let's just call it vx_massive instead of vx_sun and have to change it later.
Step 1b. Add a new section to the "Update location" part of the code
Having variables for the velocity of the massive object is great, but you will find that if you press play to run the code nothing different happens because we aren't using those variables to do anything. This is because the "Update location" section of the code is only updating the location of the planet. We need to update the location of the massive object too. You can use this code to do that:
x_massive += vx_massive*dt; y_massive += vy_massive*dt;
After copy this to your code, it should cause the sun to drift to the top right with the planet orbiting it as it goes. Pretty fun!
Step 1c. Modify drawSun so that it shows the velocity vector
In the code the drawSun function has a lot of zeros being passed to it. Here is what it looks like:
drawSun(x_massive,y_massive,0,0,0,0);
The first two zeros are for the velocity. Modify the drawSun function so it looks like this:
drawSun(x_massive,y_massive,vx_massive,vy_massive,0,0);
You should be able to see the velocity vector now if you run the code.
If you make these changes your code should behave like this
Ok, the sun is moving now, and that's good but it is not getting tugged this way or that by the planet. It is just moving in a straight line. We need to fix this! It will take a few steps to do.
The first step is figuring out the force in the x and y direction that the sun is experiencing. Once we have that we can at least show the force vector that should cause the sun to move away from the straight line. Then later we can configure the code to allow the sun to accelerate due to that force.
But it all starts with figuring out what the force is. There are two options for this: the hard way and the easy way. The hard way would be to go back to the slingshot with gravity activity and calculate the force from the planet on the sun ($F_{\rm grav} = G M m / r^2$) and then use trigonometry to calculate the components.
The easy way is to realize that for every force exerted there is an equal and opposite force (Newton's 3rd law). The code already computes the force of the sun on the planet. In the code this is Fx and Fy. The equal and opposite force to that is just -Fx and -Fy.
Let's see if this makes any sense by modifying the drawSun function. The last two zeros are for the force.
Change this:
drawSun(x_massive,y_massive,vx_massive,vy_massive,0,0);
to this:
drawSun(x_massive,y_massive,vx_massive,vy_massive,-Fx,-Fy);
Now run the code and see if the blue force arrow makes and sense. Does it at least point to the planet? Does it get larger if the planet is close to the sun?
If you make this change your code should behave like this
Earlier we added code because the location of the sun was not being updated. Now we are going to have to add code because the velocity of the sun needs to be updated.
Step 3a. Add variables for deltaVx_massive and deltaVy_massive
At the beginning of the code, add variables for deltaVx_massive and deltaVy_massive
Step 3b. Make sure the "Update velocities" section updates vx_massive and vy_massive
Remember that deltaVx_massive and deltaVy_massive don't do anything unless you use these variables to update vx_massive and vy_massive in the update velocities section
Step 3c. Calculate deltaVx_massive and deltaVy_massive
We need to set deltaVx_massive and deltaVy_massive to the appropriate values. Comparing to deltaVx and deltaVy is helpful:
deltaVx = (Fx/mass)*dt; deltaVy = (Fy/mass)*dt;
The code for deltaVx_massive and deltaVy_massive will be similar, except use -Fx and -Fy and don't use mass! That is the mass of the planet. The acceleration depends on the mass of the sun!
When you put all this in, you should notice that the sun gets a little tug in the direction of the planet and it does not move on a perfectly straight line towards the top right.
If you make these changes correctly your code should behave like this
See what happens if you set vx_massive and vy_massive equal to zero at the beginning of the code.
vx_massive = 0; vy_massive = 0;
Does the sun stay roughly in place? What happens and why?
If you make this change your code should behave like this
Try to get the sun to stay in roughly the same place by setting vy_massive to something negative.
If you choose the right value for vy_massive your code will behave like this
Optional: Use this code to measure the average velocity of the sun:
Near the beginning of the code:
vy_sum = 0; timesteps = 0;
Inside the draw function and after display()
timesteps += 1.0;
// check average velocity
vy_sum += vy_massive;
vy_ave = vy_sum/timesteps;
drawText('vy_ave = ',0.35*width,0.25*height);
drawText(vy_ave,0.5*width,0.25*height);
If you put this into your code correctly the program will behave like this
Advice: You won't be able to get vy_ave to be equal to zero. Just try to get vy_ave to vary from positive to negative around zero instead. It won't be perfect but it is more accurate than just guessing values for vy_massive and watching the sun to see if it will drift.
Why did we need to set vy_massive to something negative to prevent the system from drifting upwards? If you set the velocity of the sun to be zero initially, and then you set the velocity of the planet to be nonzero, then the total momentum is upwards (+y). What you're actually doing in Step 5 is setting vy_massive equal to some negative value until the total momentum is close to zero.
Use momentum conservation to come up with a value for vy_massive that keeps the total momentum in the +y direction equal to zero. How close does your estimate from momentum conservation compare to what you found in Step 5?
When astronomers go looking for extra-solar planets, very rarely are they able to see the planet orbiting the star. The planet is typically too dim or too close to the star for even our best telescopes to see it. Instead we have to infer the existence of the planet and the mass of the planet as best we can from observing the light from the star.
The situation is a little bit like taking our code and commenting out this line:
drawBlob(x,y,vx,vy,Fx,Fy);
change to this:
// drawBlob(x,y,vx,vy,Fx,Fy);
Now run the code and notice that you can only see the Sun moving. The planet isn't visible at all!
If you make this change your code should behave like this
If we are very far away from the star with our telescope we are not going to be able to notice the position of the star moving. The only thing we are going to be able to detect is if the star is moving towards or away from us. In a situation like this we will only be able to detect vy or vx.
Here is some code you can add at the end of the draw() function to plot vy_massive as a function of time:
graph1.addPoint(vy_massive); graph1.display();
If you put this into your code the program should behave like this
The other thing we should do is figure out the period of the motion. Let's borrow some code from the "Planetoids with a Spring!" activity to figure this out.
Add this code near the beginning where the variables are declared:
tcounter = 0; tlasttime = 0;
tcounter += dt;
drawText('counter time = ',0.7*width,0.8*height);
drawText(tcounter,0.75*width,0.75*height);
// if the x position of the planet is near the x position of the massive object
if (( abs(x - x_massive) < abs(vx)*dt) & (tcounter > dt)) {
tlasttime = tcounter;
tcounter = 0; //Reset the clock
}
drawText('half cycle time = ',0.7*width,0.65*height);
drawText(tlasttime,0.75*width,0.6*height);
With this code the program will write the "half-cycle time" to the screen. How does this relate to the period?
If you add this to your code correctly your program should behave like this
The mass units in the code are in terms of solar masses. So M = 1.0; at the beginning of the code means our star is the same mass as the sun. It is not hard to go looking in the code to figure out what the mass of the planet is. But imagine you are an astronomer and you can't just peak at nature's source code. Instead you have to figure out the mass of the planet from the motion of the star.
Here are all the ingredients you need to figure out the mass of the planet:
1. The max / min velocity of the motion of the star (Step 8)
2. The period of the motion (Step 9)
Note: this is useful because $v = \frac{2 \pi r} { T}$ where $T$ is the period.
3. Momentum is conserved (see step 6)
In Step 6. we talked about how the momentum in the y direction is $p_{\rm tot,y} = M v_{sy} + m_P v_{Py} = 0$. If we are just talking about the absolute value of the velocities then we can just say $M v_s = m_P v_P$.
4. The mass of the star ($M$)
In this case M = 1.0 which is just the mass of the sun. In a real exoplanet hunting scenario we can typically figure out the mass of the star even if its not moving because we can analyze the spectrum of the light.
5. The strength of gravity
In this case G = 100000.0 because we are using weird units and because the mass in the program is in units of solar masses instead of kilograms -- a big difference!
6. The force that the planet exerts on the star is the same as the force the star exerts on the planet
In other words: $$F_{Planet, star} = \frac{G M m_P}{r^2} = F_{star, Planet}$$
7. The acceleration of the sun is the force it experiences divided by the mass of the sun $$a_s = \frac{F_{Planet,star}}{M} = \frac{G m_p}{r^2}$$
8. The acceleration of the planet is the force it experiences divided by the mass of the planet $$a_p = \frac{F_{star,Planet}}{m_p} = \frac{G M}{r^2}$$
9. The acceleration of anything moving in a circle is $v^2 / r$
Goal: Calculate the mass of the planet in terms of what fraction of a solar mass that it is (ex. is it 0.001 or 0.01 or 0.5 ?)
Hint #1: the planet's mass won't be larger than the sun which is 1.0 in this case.
Hint #2. The mass of the planet is in the code. It's not a secret what the value is. The question is how you can figure out that number from the movement of the sun and knowing the mass of the sun and knowing how strong gravity is.
Replace the drawSun and drawPlanet commands with drawPulsar like this:
drawPulsar(x_massive,y_massive,vx_massive,vy_massive,-Fx,-Fy);
drawPulsar(x,y,vx,vy,Fx,Fy);
Now set the mass of both objects to 1.4 so that they are substantially more massive than the sun.
Finally, change the initial values of vy and vy_massive until you get the objects to have an interesting orbit like this
Fun fact #1: the orbit shown in the link is similar to a binary pulsar system that was part of the 1993 Nobel Prize in physics
Fun fact #2: All pulsars are neutron stars, but not all neutron stars are pulsars. Pulsars are neutron stars that emit a radio signal that we can detect with radio telescopes.
Most extrasolar planets have highly circular orbits like the one we've considered in this activity. But it's possible for the orbit to be more elliptical.
Go back to Steps 5 \& 6. but increase vy and make vy_massive more negative so you get something like this but don't make vy so large that the planet just flies away from the system (because the total energy would positive -- as discussed in the Escape Velocity exercise)
This will complicate the task of figuring out the mass of the planet substantially.