Calculating all gravitational forces on the earth

[MillerGenuine]

Homework Statement

I would like to calculate the gravitational net force on the Earth by the sun, and all 7 planets (mercury-neptune).

Homework Equations

G* (m1*m2)/d^2

The Attempt at a Solution

I know this is a simple case of F_12 + F_21 + F_13 + F_31 ...etc

I am trying to make a model of our solar system, so i need all gravitational forces on earth, all gravitational forces on mercury, venus, mars,...etc. So I need to do action reaction for 8 planets, but I am having a bit of trouble visualizing the summation of forces for all 8 planets.

F_earth = G* (M_sun*M_earth)/d^2 + G*(M_mercury*M_earth)/d^2 + G* (M_venus*M_earth)/d^2...etc


would this extremely tedious summation be correct?
Maybe there is a more efficient equation do calculate net force on all planets in the solar system?

 

[Delphi51]

Looks okay, but of course you must do the x and y components separately. The z components out of the plane of the solar system probably don't matter very much, but then the pull of Mercury on the Earth doesn't matter either.

If it is a computer model, I suggest a 2D array for the masses, xy positions, xy velocities, xy accelerations. Then you should be able to loop through with a
FOR i = 1 to n do for j := 1 to n do
if i <> j then ...
loop to update each object's variables. Note that i's force on j is the same as j's force on i, so only half the calcs need be done.

 

[MillerGenuine]

but of course you must do the x and y components separately

How would I do this?

And yes it is a computer model, where I have my initial velocities in the y-direction. I have my initial velocities and mass for each planet, and i loop its position and velocity update from there. I am trying to loop its force but as you can see I am having some trouble with the net force on each planet. everything runs fine, good orbits, just need to change the net force on each planet.

If it is a computer model, I suggest a 2D array for the masses, xy positions, xy velocities, xy accelerations. Then you should be able to loop through with a
FOR i = 1 to n do for j := 1 to n do
if i <> j then ...
loop to update each object's variables. Note that i's force on j is the same as j's force on i, so only half the calcs need be done.

this sounds much more efficient, any possible way you could break down how to do this?
As of right now I am just using the method i showed above for the net force on each planet, which I believe will still work, (although more time consuming) is this correct?

 

[Stonebridge]

Don't forget the Moon, its force on the Earth is very significant.
If you are intending to calculate some sort of orbital path using these forces then the Moon is essential.

 

[Delphi51]

I think you would use the distances yi-yj and xi - xj to calculate the angle θ, then do F*sin(θ) and F*cos(θ) to get the components of the force.

Just a matter of making an array [1..9,1..6] called A so that
A[1,1] is Mercury's mass, A[1,2] is Mercury's x coordinate,
A[2,3] is Venus' y coordinate and so on. The calculation of the distance between planet i and planet j would be
d = sqrt((A[i,2] - A[j,2])² + (A[i,3] - A[j,3])²)

 

[MillerGenuine]

It seems that i got it so far. thanks for your help. I am trying to add the moon, but having trouble finding the moons velocity relative to the sun. I believe its about 1050 m/s relative to the earth, and the Earth's speed is about 29000 m/s relative to the sun. Any help?

 

[Delphi51]

If you are free to choose any initial positions, you could put the Earth on the x-axis and the moon on the x-axis so the 1050 will be simply added to the 29000. Once you get your model running, you could stop it when the moon is where you want it and have it display the moon's vx and vy.

This looks like it will be fun. Add a spacecraft ?

 

[MillerGenuine]

Looks like I will have to play around with it a bit.
Yea had I used my time more wisely I would definitely enjoy writing the program much more, but i was thinking on adding an asteroid or two on a collision course for earth. Would be fun to see how the other planets pull on the asteroids.

 

[Delphi51]

You could write it for just the sun and Earth starters. But write it with expansion in mind. That matrix will make the program easy to scale up.

There is a trick in doing the changes in distance that makes it much more accurate. Use the average velocity rather than the velocity at the beginning or end of the time increment.

 

[MillerGenuine]

Here is the code i wrote up, not sure how accurate (if at all) it is. Planet[0] is the sun and planet[1] is mercury. So I just used G(m1*m2)/d^2 nine different times for each planet, then did that for all other 7 planets including the sun. Kind of tedious, but all the planets seem to orbit the sun accurately.




planet[0].Force=G*planet[1].mass*planet[0].mass / mag2(planet[1].pos - planet[0].pos)*norm(planet[1].pos - planet[0].pos)
+ G*planet[2].mass*planet[0].mass / mag2(planet[2].pos - planet[0].pos)*norm(planet[2].pos - planet[0].pos)
+ G*planet[3].mass*planet[0].mass / mag2(planet[3].pos - planet[0].pos)*norm(planet[3].pos - planet[0].pos)
+ G*planet[4].mass*planet[0].mass / mag2(planet[4].pos - planet[0].pos)*norm(planet[4].pos - planet[0].pos)
+ G*planet[5].mass*planet[0].mass / mag2(planet[5].pos - planet[0].pos)*norm(planet[5].pos - planet[0].pos)
+ G*planet[6].mass*planet[0].mass / mag2(planet[6].pos - planet[0].pos)*norm(planet[6].pos - planet[0].pos)
+ G*planet[7].mass*planet[0].mass / mag2(planet[7].pos - planet[0].pos)*norm(planet[7].pos - planet[0].pos)
+ G*planet[8].mass*planet[0].mass / mag2(planet[8].pos - planet[0].pos)*norm(planet[8].pos - planet[0].pos)
+ G*planet[9].mass*planet[0].mass / mag2(planet[9].pos - planet[0].pos)*norm(planet[9].pos - planet[0].pos)
planet[0].velocity = planet[0].velocity + planet[0].Force/planet[0].mass * dt
planet[0].pos = planet[0].pos + planet[0].velocity * dt
planet[0].trail=curve(pos=[planet[3].pos],color=planet[3].color)
print planet[0].velocity

 

[Delphi51]

I can't fully understand it because I don't know what your mag2 and norm functions do, but I really like the record structure. You could make that force calculation more elegant by doing
planet[0].Force= 0;
For j = 1 to 9 do
planet[0].Force = planet[0].Force
+ G*planet[j].mass*planet[0].mass / mag2(planet[j].pos
- planet[0].pos)*norm(planet[1].pos - planet[0].pos)

Or better yet, do all the forces in a double loop:
For i = 0 to 9 do
[planet.Force = 0
For j= 0 to 9 do if i <> j then
planet.Force = planet.Force
+ G*planet[j].mass*planet.mass / mag2(planet[j].pos
- planet.pos)*norm(planet[1].pos - planet.pos)
]

There is no speed advantage to this approach, but it does greatly reduce the quantity of code you need to carefully write and debug. If there is a bug in it, it will affect everything and be seen while a bug affecting only the force due to Pluto might never be found.

planet[0].velocity = planet[0].velocity + planet[0].Force/planet[0].mass * dt
planet[0].pos = planet[0].pos + planet[0].velocity * dt

would work better if you used an average velocity:
v = planet[0].velocity
planet[0].velocity = planet[0].velocity + planet[0].Force/planet[0].mass * dt
planet[0].pos = planet[0].pos + (v +planet[0].velocity)/2 * dt

How in the devil are you doing this planet[0].pos update so that it does both the x and y positions?