Parameters for Solar system simulation

[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.

 

[CAF123]

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

That's me being silly late last night.

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.

That's interesting. I think my project is certainly under the category of small simulation, yet when I use velocities in km/s and position in km via the horizon interface, I cannot observe the true behavior of the bodies. I probably had the scaling wrong but I am going to implement the new unit system soon.
The horizon interface gives velocities in terms of AU/days. Since VU = DU/TU, VU = AU/58 (approx) days. So if I understand correctly, is it just a case of dividing all entries of velocity by 58?
Thanks.

 

[gneill]

That's me being silly late last night.


That's interesting. I think my project is certainly under the category of small simulation, yet when I use velocities in km/s and position in km via the horizon interface, I cannot observe the true behavior of the bodies. I probably had the scaling wrong but I am going to implement the new unit system soon.

The horizon interface gives velocities in terms of AU/days. Since VU = DU/TU, VU = AU/58 (approx) days. So if I understand correctly, is it just a case of dividing all entries of velocity by 58?
Thanks.

$$\frac{AU}{day} \times \frac{58.133 day}{TU} \rightarrow \frac{AU}{TU}$$

So it looks like you want to multiply by 58.133.

 

[CAF123]

$$\frac{AU}{day} \times \frac{58.133 day}{TU} \rightarrow \frac{AU}{TU}$$

So it looks like you want to multiply by 58.133.

I see. Last question for the moment: Why is using VU=DU/TU valid? I know d=vt, but this only holds for constant speed bodies.

 

[gneill]

I see. Last question for the moment: Why is using VU=DU/TU valid? I know d=vt, but this only holds for constant speed bodies.

Speed is defined to be the distance covered per unit time. It's the same idea as km/hr, where you're taking km as the distance unit and hr as the time unit. Regardless of whether a given motion is constant or not, the definition of velocity holds.

If you'd rather, define the VU as:
$$VU = \sqrt{\frac{\mu}{DU}}$$
You'll get the same value

 

[CAF123]

Yes, that makes sense. Thanks again gneill. I'll implement this new unit system and see if my simulation improves.

 

[CAF123]

Hi gneill,
I implemented the new system and it works for one timestep but not some others. Is there a reason why there exists a maximum time step? Is it the case that if you go above this max time step, then it is likely that a 'turning point' in the motion has been missed?

I also have to deduce the angle Halley's comet makes with the ecliptic. I googled for a known value and it appears it is about 18°. My method was to compute the arcsin of the z component of position over the total magnitude of position. This gave me a value of about 11°. Is this method okay and if so, I believe the difference would be due to the fact that by doing the trig like I did, I am assuming the comet approached in a straight line.

 

[gneill]

Hi gneill,
I implemented the new system and it works for one timestep but not some others. Is there a reason why there exists a maximum time step? Is it the case that if you go above this max time step, then it is likely that a 'turning point' in the motion has been missed?

That's possible. Integrations can go awry if timesteps "step over" significant changes in the trajectory, like "missing" perihelion. Try adjusting your timestep size downwards and see what happens.

I also have to deduce the angle Halley's comet makes with the ecliptic. I googled for a known value and it appears it is about 18°. My method was to compute the arcsin of the z component of position over the total magnitude of position. This gave me a value of about 11°. Is this method okay and if so, I believe the difference would be due to the fact that by doing the trig like I did, I am assuming the comet approached in a straight line.

That method would work if the position was perihelion or aphelion, then it would give you the orbital inclination. Can you see why this is so? Hint: imagine that the position happened to coincide with the ascending or descending node.

 

[CAF123]

That method would work if the position was perihelion or aphelion, then it would give you the orbital inclination. Can you see why this is so? Hint: imagine that the position happened to coincide with the ascending or descending node.

At the perihelion or aphelion, the total velocity of the particle is solely in the tangential direction. What do you mean by the ascending/descending node?

 

[gneill]

At the perihelion or aphelion, the total velocity of the particle is solely in the tangential direction. What do you mean by the ascending/descending node?

If the orbit is inclined to the ecliptic it must pass through the ecliptic at two points. Once heading "north" and once heading "south". Those are the ascending and descending nodes. The line of nodes is the line joining those two points.

 

[CAF123]

If the orbit is inclined to the ecliptic it must pass through the ecliptic at two points. Once heading "north" and once heading "south". Those are the ascending and descending nodes. The line of nodes is the line joining those two points.

So when it is instantaneously at the point where it crosses the ecliptic, position vector is parallel to this plane and so so the angle we would obtain would be 0. But what is stopping me from using any other point that does not coincide with the nodes? I.e a point just before the aphelion/perihelion or a point just after the node? Is it because at the aphelion/perihelion, the motion is at a turning point, we know exactly what the orbit looks like and hence computing it here gives the best estimate of the angle?
Thanks.

 

[gneill]

So when it is instantaneously at the point where it crosses the ecliptic, position vector is parallel to this plane and so so the angle we would obtain would be 0. But what is stopping me from using any other point that does not coincide with the nodes? I.e a point just before the aphelion/perihelion or a point just after the node? Is it because at the aphelion/perihelion, the motion is at a turning point, we know exactly what the orbit looks like and hence computing it here gives the best estimate of the angle?
Thanks.

It's purely a geometrical effect. You will get a continuum of angle values from zero at the the nodes to a maximum at peri/aphelion.

Orbits lie in a plane. The intersection of two planes is a line, and as it happens, it coincides with the line of nodes I mentioned before. You are looking for the angle of intersection of that orbital plane with the ecliptic plane; that's the inclination angle.

An orbital position vector at perihelion or aphelion will be at right angles to the line of nodes and lie in the plane of the orbit. The angle between that position vector and the ecliptic plane will thus be the inclination of the orbital plane.

If you have both the position and velocity vectors for a single instant you can determine the orbit and all its parameters. It is handy to know that what is called the eccentricity vector is a vector whose magnitude equals the eccentricity of the orbit and points in the direction of perihelion. It's an ideal candidate for determining the orbit inclination

 

[CAF123]

I attached a sketch of what I observed from my simulation for the angle of the comet. I also have the magnitude of the aphelion and if I search through my output files, I should be able to extract the z component of the total position. I think I can then find the angle as mentioned before.
Is this okay?
I realize that the position vector here may not lie exactly along the orbit of the comet, (despite what I have drawn), but would it be a reasonable approximation? Unfortunately, the work is due soon, so I probably wouldn't be able to implement the eccentricity, but it would be something to perhaps add in the report for further work.

 

[gneill]

I attached a sketch of what I observed from my simulation for the angle of the comet. I also have the magnitude of the aphelion and if I search through my output files, I should be able to extract the z component of the total position. I think I can then find the angle as mentioned before.
Is this okay?

That should work. Note that you could also use the angle between the velocity vector and the ecliptic plane when the body is at one of the nodes (detect the plane crossing (z crosses from + to - or vice versa at the nodes) and use the then current velocity vector).

I realize that the position vector here may not lie exactly along the orbit of the comet, (despite what I have drawn), but would it be a reasonable approximation? Unfortunately, the work is due soon, so I probably wouldn't be able to implement the eccentricity, but it would be something to perhaps add in the report for further work.

You mean the eccentricity vector? It just requires knowledge of one pair of r,v and a few vector operations. Your initial input values for r,v would suffice.

 

[CAF123]

Given the eccentricity vector, would I then follow the same procedure as before (calculate it's magnitude, divide by z component and take arcsin). In what way is this method better than using the position vector corresponding to the perihelion?

I made plots of energy difference between the total energy at some arbitrary point in the system and the initial energy vs time. (the energy difference is small) The form of these graphs is roughly sinusoidal and about symmetric with respect to the time axis. Why should these graphs be roughly sinusoidal? Is it just the fact that the max of the graphs represents the system at it's peak of kinetic energy and the min points, where the potential energy is greatest? (Something like a harmonic oscillator)

 

[gneill]

Given the eccentricity vector, would I then follow the same procedure as before (calculate it's magnitude, divide by z component and take arcsin). In what way is this method better than using the position vector corresponding to the perihelion?

The eccentricity vector can be determined starting with any r,v pair. As a bonus you also get the eccentricity of the orbit and the direction to perihelion.

I made plots of energy difference between the total energy at some arbitrary point in the system and the initial energy vs time. (the energy difference is small) The form of these graphs is roughly sinusoidal and about symmetric with respect to the time axis. Why should these graphs be roughly sinusoidal? Is it just the fact that the max of the graphs represents the system at it's peak of kinetic energy and the min points, where the potential energy is greatest? (Something like a harmonic oscillator)

The total energy of an unperturbed orbit should be a constant. Kinetic and potential energy should be exchanged perfectly and without losses or gains in the total. Deviations from a constant value represent artifacts of the integration algorithm. A perfect integrator would yield a straight line with no deviations. If you've got sinusoidal deviations that don't trend upwards or downwards but return to the starting value periodically, then the integrator isn't losing or "creating" energy over a full cycle (so your orbits aren't going to "run away").

Simple integrators like the Euler, leapfrog or Verlet algorithms are pretty good but not perfect. The leapfrog algorithm is remarkable for its simplicity and ability to conserve energy for orbital mechanics despite this simplicity, but as you can see, it is still not perfect.

 

[CAF123]

Thanks,

The eccentricity vector can be determined starting with any r,v pair. As a bonus you also get the eccentricity of the orbit and the direction to perihelion.

But what is special about this vector (i.e what makes it different from the position vector of the perihelion). Do I extract the ecliptic angle in the same way I did previously, but this time using the eccentricity vector?

 

[gneill]

Thanks,

But what is special about this vector (i.e what makes it different from the position vector of the perihelion). Do I extract the ecliptic angle in the same way I did previously, but this time using the eccentricity vector?

Nothing special other than it hands you the extra information, which can be useful if you're trying to determine the canonical orbital elements for an orbit from position and velocity.

Sure, that works.

 

[CAF123]

Nothing special other than it hands you the extra information, which can be useful if you're trying to determine the canonical orbital elements for an orbit from position and velocity.

Sure, that works.

I am using data from this page: http://ssd.jpl.nasa.gov/horizons.cgi#results
This gives position and velocity in km and km/s respectively so to compute the eccentricity, I converted everything to SI. My problem is when I do the computation of the eccentricity, via the two equations in Wikipedia, I end up with both times e = <-0.545, -0.746, 0.22> (which I confirmed with Wolfram alpha) and this means the arcsin does not exist.

 

[gneill]

I am using data from this page: http://ssd.jpl.nasa.gov/horizons.cgi#results
This gives position and velocity in km and km/s respectively so to compute the eccentricity, I converted everything to SI. My problem is when I do the computation of the eccentricity, via the two equations in Wikipedia, I end up with both times e = <-0.545, -0.746, 0.22> (which I confirmed with Wolfram alpha) and this means the arcsin does not exist.

Can you post the position and velocity values here? I won't see the "results" that you see on the link you gave; it'll show a results page from my last visit there.

 

[CAF123]

Can you post the position and velocity values here? I won't see the "results" that you see on the link you gave; it'll show a results page from my last visit there.

Sorry about that, here they are:
-3.056663785023214E+09 3.715913718019426E+09 -1.453046492475981E+09 (Position x,y,z km)
-2.107989578404891E-01 1.729355547234978E+00 -3.417437459119960E-01 (Velocity vx,vy,vz km/s)

 

[gneill]

Okay, using those values I find $\vec{e} = < 0.54757, -0.74966, 0.27407 >$. Its magnitude is 0.96796. arcsin(0.27407/0.96796) certainly has a value...

 

[CAF123]

Okay, using those values I find $\vec{e} = < 0.54757, -0.74966, 0.27407 >$. Its magnitude is 0.96796. arcsin(0.27407/0.96796) certainly has a value...

Sorry gneill for wasting your time, I was doing arcsin(0.96796/0.27401) despite the fact that I had written it down correctly. 16^o is quite a good result, I believe the actual result is about 18^o. Given that the method used to find this angle is not an approximation, the differences can be attributed to the integration algorithm I used?

 

[gneill]

Sorry gneill for wasting your time, I was doing arcsin(0.96796/0.27401) despite the fact that I had written it down correctly. 16^o is quite a good result, I believe the actual result is about 18^o. Given that the method used to find this angle is not an approximation, the differences can be attributed to the integration algorithm I used?

Well, I assume that these initial position and velocity values are the real deal from Horizons, and they never get near your integrator for this determination of the angle. So any discrepancy with a 'real' value would have to have another source. A couple of things come to mind. First, the inclination of a comet can change over time due to perturbations (particularly by Jupiter). Second, the coordinate system used may not correspond to that in which the original 18° figure was determined; in particular the Fundamental Plane (ecliptic plane) changes over time as the Earth's orbit evolves, and the Horizons values were probably for a barycentric coordinate system rather than a Sun-centered one. Any and all of these things could contribute to the difference. They'd have to be tracked down individually to stamp them out.