Parameters for Solar system simulation

[CAF123]

Homework Statement

I am constructing an N-body simulation in Java that will ultimately model the solar system. I am testing the reliability of it by slowly increasing the number of bodies in the simulation. So I am starting off with the Sun, Mercury and Venus. I am looking to find the initial position, velocity and masses of these three planets. Masses are easy enough to find, however I still need to find the initial position and velocity.

Homework Equations

I am told to find the initial velocity of the planets via, $$v = \sqrt{\frac{\mu}{r}},$$where ##\mu = m_1 + m_2## and ##G = 1## and I think ##r## is the distance of body 1 from body 2.

The Attempt at a Solution

If ##\mu = m_1 + m_2## then if one of the masses is the Sun, this will dominate the other planet. I think I can find the mean distance of each of Venus and Mercury from the Sun online. (Can I take this to be r?) So I can use the above formula to obtain the initial velocity. I am a bit confused though because this is a 3 body system, so how should I compute the resultant initial velocity using that formula?

Also, I can see that this eqn was derived redefining G=1, (that makes the equation dimensionally correct) but this does not alter the units of the quantities.

Many thanks

 

[gneill]

Spend a bit of time tinkering with the JPL Horizons Web Interface.

The Horizons system allows you to generate ephemerides for anybody in the solar system from just about any viewpoint you can name (use @sun to generate an ephemeris for an "observer" at the center of the Sun). Pick a particular date for your epoch, pick a target body (say Venus), and out pops an ephemeris with time-stamped coordinates as "seen" by the observer.
[STRIKE]
I think it's possible to change the outputs so that you can get rectangular coordinates, too.[/STRIKE] Nope. I checked.

Good luck!

 

[mfb]

I think I can find the mean distance of each of Venus and Mercury from the Sun online. (Can I take this to be r?)

That is good enough if you are not interested in details of the orbits (they are not circular, for example). If you want realistic orbits, use actual velocities and positions for the planets, there are databases that provide those values.

I am a bit confused though because this is a 3 body system, so how should I compute the resultant initial velocity using that formula?

The forces between the planets are tiny compared to the forces between planets and sun.

Also, I can see that this eqn was derived redefining G=1, (that makes the equation dimensionally correct) but this does not alter the units of the quantities.

It used a unit system where the units match. This can have weird consequences (like masses expressed as powers of length and time), but it works.

 

[CAF123]

That is good enough if you are not interested in details of the orbits (they are not circular, for example). If you want realistic orbits, use actual velocities and positions for the planets, there are databases that provide those values.

Do you have reference to a link? I tried searching for this, but I could only find the current velocities/position. For what I am doing, I think I have to use the expression given, so I guess the mean distance of the planets from the Sun will do.

The forces between the planets are tiny compared to the forces between planets and sun.

This makes sense, so then I could just apply the equation twice - for Mercury around the Sun and for Venus around the Sun.

It used a unit system where the units match. This can have weird consequences (like masses expressed as powers of length and time), but it works.

So I would need to find the masses of the planets in terms of length and time. So [v] = m/s and [G] = [m³kg^-1s^-2]. Hence I should express my masses like [mass] = [length³ time^-2] How would I obtain values like this so I can then use my formula?

 

[D H]

[STRIKE]
I think it's possible to change the outputs so that you can get rectangular coordinates, too.[/STRIKE] Nope. I checked.

You didn't check well enough. Change the "Ephemeris Type" (the first selection) from "Observer" to "Vector table".

 

[gneill]

You didn't check well enough. Change the "Ephemeris Type" (the first selection) from "Observer" to "Vector table".

Oh ho! Very nifty. Hadn't seen that before. Good info! Thanks.

 

[CAF123]

Spend a bit of time tinkering with the JPL Horizons Web Interface.

The Horizons system allows you to generate ephemerides for anybody in the solar system from just about any viewpoint you can name (use @sun to generate an ephemeris for an "observer" at the center of the Sun). Pick a particular date for your epoch, pick a target body (say Venus), and out pops an ephemeris with time-stamped coordinates as "seen" by the observer.
Good luck!

Thank you, I presume what they mean by 'initial velocity' and 'initial position' are those quantities when the planets begun their orbits around the Sun. So what epoch is this?

My code is implemented in rectangular coordinates (that is I have seven parameters for each particle - label, 3 position coordinates, 3 velocity coordinates).

I am still trying to navigate around the Horizon interface and when I have more time, I'll look at the user guide. My instructor actually provided the equation I mentioned in the OP, so I think I should use it. Thanks.

 

[gneill]

Thank you, I presume what they mean by 'initial velocity' and 'initial position' are those quantities when the planets begun their orbits around the Sun. So what epoch is this?

You can pick any instant in time as your epoch. Use the same instant for all the bodies and the positions and velocities will all be accurate to start your simulation from that instant.

My code is implemented in rectangular coordinates (that is I have seven parameters for each particle - label, 3 position coordinates, 3 velocity coordinates).

I am still trying to navigate around the Horizon interface and when I have more time, I'll look at the user guide. My instructor actually provided the equation I mentioned in the OP, so I think I should use it. Thanks.

If you assume circular orbits then that's fine. But your orbits will all be circular rather than elliptical, and you still have the problem of choosing suitable starting locations. Using actual ephemeris positions and velocities will give you the real thing.

Here's some sample output for an ephemeris for Mars:

 

[mfb]

So I would need to find the masses of the planets in terms of length and time. So [v] = m/s and [G] = [m³kg^-1s^-2]. Hence I should express my masses like [mass] = [length³ time^-2] How would I obtain values like this so I can then use my formula?

If you use F=Mm/r^2 as formula this is fine, if you use F=GMm/r^2 you can multiply all "mass" values with some power of G and whatever to get masses in kg.

 

[D H]

So I would need to find the masses of the planets in terms of length and time.

That's exactly right. That's what your μ was expressed in in the opening post. Look at the units. You had ##v=\sqrt{\mu/r}##. That means μ has to have dimensions of length³/time². Solar system modelers don't use G*M. They use μ. The reason is that a planet's "standard gravitational parameter", or mu, can be inferred from planetary observations. Planetary mass is computed from μ/G. Here's the problem: G is one of the least well known physical constants. This uncertainty in G means that it's best if G never enters into the equations.

 

[mfb]

Here's the problem: G is one of the least well known physical constants. This uncertainty in G means that it's best if G never enters into the equations.

This also means M is poorly known (well, relative to other measurements). Expressed in other words: the uncertainty on the mass of our sun is larger than the mass of earth.
We really need a space-based experiment to determine G...

 

[CAF123]

That's exactly right. That's what your μ was expressed in in the opening post. Look at the units. You had ##v=\sqrt{\mu/r}##. That means μ has to have dimensions of length³/time². Solar system modelers don't use G*M. They use μ. The reason is that a planet's "standard gravitational parameter", or mu, can be inferred from planetary observations.

##\mu = G(m_1 + m_2) = m_1 + m_2 \approx m_1 = m_s##, where ##G =1## and ##m_s## is the mass of the Sun and this dominates all the other masses under consideration. How do I obtain this value of the mass of the Sun in terms of units of length^3/time^2?
]

 

[mfb]

##\mu = G(m_1 + m_2) = m_1 + m_2 \approx m_1 = m_s##, where ##G =1## and ##m_s## is the mass of the Sun and this dominates all the other masses under consideration. How do I obtain this value of the mass of the Sun in terms of units of length^3/time^2?

You know M, you want to know GM. Just multiply the mass with G (expressed in the SI system).

 

[D H]

##\mu = G(m_1 + m_2) = m_1 + m_2 \approx m_1 = m_s##, where ##G =1## and ##m_s## is the mass of the Sun and this dominates all the other masses under consideration. How do I obtain this value of the mass of the Sun in terms of units of length^3/time^2?

Google the term standard gravitational parameter. That's the name for the product of G and mass.

 

[CAF123]

I see, I was overcomplicating things. Now that I have the magnitude of the initial velocity for Mercury around the Sun (computed using the expression in OP) in order to put this into the program I need to supply the x,y,z components of initial velocity in an input file.

To find these, I was thinking I could use gneill's suggestion of the Horizon interface since I can see a key on the ephemerides about half way down 'JDCT X Y Z <newline> VX VY VZ'. If I understand this correctly, the first row lists the components of position and the second row lists the components of velocity. I can use this to get my initial position. Calculating gives a 0.14% difference in the value of the initial position using Horizon and the value of the mean distance.

Similarly, for the velocity I obtain a 6.9% difference. I think the only reason for this difference is the propagation of the error when I take the sqrt of ##\mu/r## since the differences in r are barely noticeable. Is this right? And thanks gneill for providing that link - it really helped.

Thanks everyone.

 

[gneill]

Another (small) source of difference would be that the Horizons data has the solar system barycenter as its origin. This is close to, but not identical with, the center of the Sun (occasionally it lies outside of the Sun). The barycenter is determined by the instantaneous location of all the planets, not just the one you're looking at at the moment!

If you're modelling the Sun as free body, too, you'll want to give it initial location and velocity parameters.

 

[CAF123]

Another (small) source of difference would be that the Horizons data has the solar system barycenter as its origin. This is close to, but not identical with, the center of the Sun (occasionally it lies outside of the Sun). The barycenter is determined by the instantaneous location of all the planets, not just the one you're looking at at the moment!

In the coordinate origin data entry, I changed this from the default 'Solar system barycenter' to the Sun. But see comment below; I may change this.

If you're modelling the Sun as free body, too, you'll want to give it initial location and velocity parameters.

Since I set the Sun at the origin, I think I can let the initial position of the Sun to be at <0,0,0>. How would I obtain the velocity of the Sun? This will change depending on the point of observation. Since I set the Sun as my origin, the velocity of the Sun in this 'frame' is trivially zero. Do you think it is better, therefore, to assign values to everything with the Solar system barycenter as my origin?

Thank you.

 

[gneill]

Since I set the Sun at the origin, I think I can let the initial position of the Sun to be at <0,0,0>. How would I obtain the velocity of the Sun? This will change depending on the point of observation. Since I set the Sun as my origin, the velocity of the Sun in this 'frame' is trivially zero. Do you think it is better, therefore, to assign values to everything with the Solar system barycenter as my origin?

Thank you.

I would say yes. Otherwise you may find your entire simulation "creeping off the page" as time goes on. Why? Because your coordinate system origin won't be located at the center of momentum of the system, wherever it turns out to be once the 'physics' get going in the simulation. Unless the Sun happens to be at the true center of momentum for the ensemble, it's going to wander.

Two choices:
1) Set initial positions and velocities for all bodies, including the Sun, with the barycenter as the origin.
2) Set initial positions and velocities for all bodies and assume the Sun is motionless at the origin. THEN run through the data to determine the location and velocity of the center of mass of the ensemble. THEN tweak all the velocity values to make the center of mass stationary.

 

[CAF123]

So do you mean to say about the barycenter the total mechanical momentum of the system is conserved? Choice 1 seems simpler given I can use the Horizon interface.

 

[gneill]

So do you mean to say about the barycenter the total mechanical momentum of the system is conserved? Choice 1 seems simpler given I can use the Horizon interface.

It's a closed system, so yes, momentum is conserved and the motion of the center of mass is inertial (constant velocity). You can choose your coordinate system to 'ride along with' this center of momentum, thus keeping your simulation "pinned" to the center of the coordinate system. Otherwise the center of mass of the system would drift linearly with time. It can be annoying to watch an interesting simulation gradually drift off the screen!

Yup, choice 1 is a good way to go. Even so, I would still recommend, as a matter of course, tweaking the velocities after assembling the system to zero any residual motion of the center of mass. This is because perhaps not all the existing solar system bodies will be included for every run, and the initial positions/velocities given by Horizons take all bodies into account. Imagine if Jupiter and the other gas giants are left out but the Horizon's location and velocity for the barycenter includes their influence...

Just find velocity vector for the center of mass (add all the momentum vectors, divide by total mass to yield velocity vector of center of mass). Subtract this vector from all the initial velocity vectors. Done.

 

[CAF123]

I extracted the masses, initial positions and initial velocities of the Sun, Mercury and Venus about the Solar system barycenter from the Horizon interface. I am wondering though, I implemented G=1 in my code and so I think this would vary the magnitudes of these masses, positions and velocities. Is this right? If so, I am struggling to come up with suitable scale factors.

I am meant to assume a circular orbit to begin with so ##m_1v^2/r = Gm_1m_2/r^2## holds. This means ##v = \sqrt{Gm_2/r} = \sqrt{\mu/r}##. Since G=1, to obtain units of m/s for velocity, should I simply scale the mass of the Sun provided by the Horizon interface by 6.67*10^-11?

Is this the best way? I am not sure how the other values from the interface would be affected though.
Thanks.

 

[CAF123]

I am a just bit unsure of how to get the scaling correct really. I have supplied simple numbers to the code to test that it is working and it does, but now I need to get the scaling correct when I use the data from the horizon interface.
Thanks.

 

[CAF123]

gneill, can you help at all or is my question unclear? I am now thinking of perhaps using astronomical units to make the quantities bigger so that I can see the simulation. For random numbers of position, velocity = √(μ/r), and mass=100, the simulation was good for 3 bodies. When I implement real data via Horizon, I get a what looks like repulsive force between the bodies. So, I think I should scale the data in such a way (which is my problem) to observe the true behaviour of the 3 bodies.
Thanks.

 

[gneill]

gneill, can you help at all or is my question unclear? I am now thinking of perhaps using astronomical units to make the quantities bigger so that I can see the simulation. For random numbers of position, velocity = √(μ/r), and mass=100, the simulation was good for 3 bodies. When I implement real data via Horizon, I get a what looks like repulsive force between the bodies. So, I think I should scale the data in such a way (which is my problem) to observe the true behaviour of the 3 bodies.
Thanks.

Here's a suggestion for a system of units for a simulation. Suppose we call them MU, DU, TU, VU. Right away you can pick MU to be the mass of the Sun and DU to be an astronomical unit. Then use the known value of ##\mu_{sun}## to set TU:
$$TU = \sqrt{\frac{DU^3}{\mu_{sun}}}$$
VU is simply DU/TU.

If you set suitable constants to these values then you can do your conversions back and forth.

Note that in these units the gravitational parameter and gravitational constant become:
$$\mu = 1\;\frac{DU^3}{TU^2} \\
G = 1 \frac{DU^3}{MU\cdot TU^2}$$
also note that "in real life", TU is about 58.133 days and multiplied by ##2\pi## is one year. The velocity unit VU is equal to the average orbital speed of Earth. So while the mass unit is Sun-based, the other units scale by values related to the Earth orbit. Typical of Man's geocentricity I suppose

 

[CAF123]

Thanks gneill, what does the units MU, DU, TU and VU represent? I figured MU is the 'mass unit' which you defined by 1 solar mass. 'VU' is the velocity unit, TU = 'time unit' and 'DU' = distance unit?
Thanks.

 

[D H]

Right.

Another way to look at it: 1 time unit is 365.25/(2*pi) days (astronomers use the julian year), 1 length unit is the astronomical unit, and 1 mass unit is the solar mass.

Note: The value of G is not quite 1 with these units. It's very close to one, and in the end it doesn't matter because G should never enter into the picture. You should be using the product of G and mass, the standard gravitational parameter or mu (each body has it's own mu); rather than G*mass as two separate quantities.

 

[CAF123]

Right.

Another way to look at it: 1 time unit is 365.25/(2*pi) days (astronomers use the julian year), 1 length unit is the astronomical unit, and 1 mass unit is the solar mass.

How did you deduce that TU = √DU³/μ_s is equivalent to 365.25/(2*pi)? I was wondering also, in the Horizon interface, there is an in built time unit of AU/days for the velocity of the bodies. Is this the same unit here?

Note: The value of G is not quite 1 with these units.

In SI, G = 6.67*10^-11kg^-1m³s^-2. If we initially define G = 1 DU³/(MU*TU²) then the expressions for the rest of the units can be derived exactly, so why would G not be 1?

 

[gneill]

How did you deduce that TU = √DU³/μ_s is equivalent to 365.25/(2*pi)? I was wondering also, in the Horizon interface, there is an in built time unit of AU/days for the velocity of the bodies. Is this the same unit here?

The formula is Kepler's third law with μ as the proportionality constant. I believe that the ##2\pi## has to do with the fact that we're defining the units based upon Earth's "circular" orbit. The DU, VU and TU are all related by that circle.

 

[CAF123]

Is the time unit in built into the Horizon interface the same one you proposed? It gives an option to express length in AU and velocities in AU/days, which given the discussion about the equivalence of TU in terms of 365.25/2*pi days, makes me wonder.

 

[gneill]

Is the time unit in built into the Horizon interface the same one you proposed? It gives an option to express length in AU and velocities in AU/days, which given the discussion about the equivalence of TU in terms of 365.25/2*pi days, makes me wonder.

No, not the same units. AU/day is a convenience for astronomers since distances in the solar system are generally expressed in AU and the day as a unit of time suits many of the formulae used in the industry (not to mention daily data acquisition schedules!).

 

[CAF123]

Ok here is my analysis:
In SI, G = 6.67*10^-11kg^-1m³s^-2. Define 1 MU = 10³⁰kg => kg^-1 = 10³⁰(MU)^-1.

Define 1 AU = 1DU = 149597871km = (p*10³)m = αm => m³ = (1 DU)³/α³.

Subbing into the required form;

1 DU³MU^-1TU^-2 = 6.67*10^-11 (10³⁰)MU^-1 (DU/α)³s^-2. Rearrange for TU in terms of s gives TU = s/21115.4. This means my unit of time is no longer the second, but this fraction of a second. Is this right given the definitions of DU and MU?

So then VU = DU/TU = 1 AU/s γ, where γ = 21115.4. Does this mean when I take data from Horizon, I need to put all positions into AU (which is easy - it does it for me) and convert all velocities given in AU/day into this unit for VU. How would I do that if that is correct?

Thanks.

 

[D H]

Just don't use G. You do not need it. Here are the values for GM in units of (AU)³/day² from JPL's DE405 (the DE421 values are slightly different):

• Mercury - 4.91254745145081187e-11
• Venus - 7.24345248616270166e-10
• Earth+Moon - 8.99701134671249882e-10
• Mars - 9.54953510577925806e-11
• Jupiter - 2.82534590952422643e-07
• Saturn - 8.45971518568065874e-08
• Uranus - 1.29202491678196939e-08
• Neptune - 1.52435890078427628e-08
• Pluto - 2.18869976542596968e-12
• Sun - 2.95912208285591095e-04

 

[gneill]

The MU, DU, TU, VU units have corresponding values in metric units. So it's an exercise in unit conversion. if you're given a velocity in AU/day, and you want to convert it to DU/TU, then you know that DU = AU (they were defined to be the same) and TU = 58.133 days. So:
$$ <given\; \frac{AU}{day}> \times \frac{1\;DU}{1\;AU} \times \frac{58.133\;day}{TU} = ~~<result> \frac{DU}{TU}$$

and DU/TU = VU so you're done.

 

[CAF123]

The MU, DU, TU, VU units have corresponding values in metric units. So it's an exercise in unit conversion. if you're given a velocity in AU/day, and you want to convert it to DU/TU, then you know that DU = AU (they were defined to be the same) and TU = 58.133 days. So:
$$ <given\; \frac{AU}{day}> \times \frac{1\;DU}{1\;AU} \times \frac{58.133\;day}{TU} = ~~<result> \frac{DU}{TU}$$

and DU/TU = VU so you're done.

I derived in my last post that TU = s/21115.4. Can I show that this is equivalent to the 58.133 days?

 

[gneill]

I derived in my last post that TU = s/21115.4. Can I show that this is equivalent to the 58.133 days?

I don't see how a fraction of a second can be the same as 58 days...

Sorry, I didn't follow what you were doing in that post.

Is the problem that you're having difficulty doing unit conversions to the internal units? If so, why not go through an intermediate step; convert all given values to base metric units (kg, m, s) to begin with. Then its just a matter of multiplying or dividing by the appropriate conversion constant to convert to MU, DU, TU,...

For example, the 58.133 days of the TU is, in seconds, 5022732 s. If you divide a 'real' time in seconds by 5022732 then you have that time in TU. Similarly, if you divide a 'real' velocity given in meters per second by 29784.5 then you'll have that velocity in VU.

The idea of having your own set of units for the simulation is strictly to make some of the calculations faster by making often used constants have numerical values of 1. These days with faster computers with 64 or 128-bit arithmetic it's not as big a deal for small simulations, and one can work directly in kg,m,s values if one wishes.