Java simulation of solar system

[CAF123]

Homework Statement

I am creating an astrophysical simulation involving an arbritary number of particles and each particle is identified via a mass, label, position vector and velocity vector. I have made classes which can perform operations on these vectors and I also have a class which describes a particle in 3-D space (using this vector class).

With this, I have outlined an approach to calculate the force on some particle due to the others, the total energy of the system and methods to time-evolve the particle (LeapFrog algorithms). For the moment, we have to reduce the system to only 3 particles (Venus, Sun, Mercury) One of the questions was:

Count how many times each of the planets orbits the Sun during the simulation. First create an algorithm to test for a complete orbit. Do the two planets’ years have the correct ratio?

The Attempt at a Solution

So I need to find a way to test for a complete orbit. I understand this to mean that the particle passes the same point again. The number of orbits in the simulation will depend on the amount of time that the program is run for so I thought I could divide this by the period of the particle in question. I can't see how to do this because the path of the particle need not be circular and later we consider elliptical orbits anyway. A hint in the right direction would be great, thanks.

 

[Mark44]

Whether the orbit is circular (not realistic) or elliptic (realistic), if a planet passes through a given point at some time, it will pass through the same point after it completes one orbit, right? I'm assuming that the planetary orbit doesn't precess (I think that's the right term) or change shape.

 

[UltrafastPED]

Consider the coordinate system which is being used:
http://en.wikipedia.org/wiki/Orbital_mechanics

If you are using rectangular coordinates you could transform to spherical and test the angle against the starting angle; you would only need to do this when your orbital time elapsed has approached the orbital period - then test the angles.

 

[CAF123]

Hi Mark44,

Whether the orbit is circular (not realistic) or elliptic (realistic), if a planet passes through a given point at some time, it will pass through the same point after it completes one orbit, right?

Yes, that is my understanding. For a circular orbit, $r=\text{const} \Rightarrow T=2\pi r/v$. I can find r by the magnitude of the position vector and v I know from the LeapFrog algorithm. However, these eqns are true for circular motion only and not elliptical. I am also not sure if this is the right or best approach.

I'm assuming that the planetary orbit doesn't precess (I think that's the right term) or change shape.

I think we have to assume that the orbit follows one particular path, i.e doesn't move anticlockwise and then clockwise later on.

 

[CAF123]

Hi UltrafastPED,

If you are using rectangular coordinates you could transform to spherical and test the angle against the starting angle;

I am using rectangular coordinates, yes. If we place the sun at the origin initially (initial position at r=0 and with some velocity) then the planet is at a distance r away in sphericals. What angle are you referring to? Do you mean if the angles $\phi$ and $\theta$ coincide after one revolution?

you would only need to do this when your orbital time elapsed has approached the orbital period - then test the angles.

I don't think we know the orbital period.

 

[UltrafastPED]

You would need to impose a new coordinate system centered on the central star; your planetary orbitals will stay on their original plane (unless you have close encounters, or collisions). So for a given planet construct a circular coordinate system; then the angle coordinate from this system is your reference.

You can do this for each planet.

There are certainly other ways to do it, but it will involve some manipulation of the coordinates ... so use the most convenient system!

Note: orbital periods don't change - see Kepler's laws - so once you have completed an orbit, you can cut your computational load by tracking the "days of the year".

 

[CAF123]

You would need to impose a new coordinate system centered on the central star; your planetary orbitals will stay on their original plane (unless you have close encounters, or collisions). So for a given planet construct a circular coordinate system; then the angle coordinate from this system is your reference.

Ok, so what you mean is that essentially angle $\theta$ is fixed, ($\theta$angle from z axis) so we do not need to worry about it. The reference angle is the angle $\phi$ in the xy plane if I use SP's?

So I want to show that at a later time, the angle relative to the coordinate system and the distance of the planet from the central star is the same as a point before. This makes a complete orbit.

 

[D H]

Whether the orbit is circular (not realistic) or elliptic (realistic), if a planet passes through a given point at some time, it will pass through the same point after it completes one orbit, right? I'm assuming that the planetary orbit doesn't precess (I think that's the right term) or change shape.

Even elliptical orbits are not realistic in the N-body problem. Just look at the "year". There are sidereal years, tropical years, anomalistic years, heliacal years, Gaussian years, Besselian years, etc. One year and the next are not the same spans of time. It's quite a mess!CAF123, what I think you want is the sidereal year. Presumably your planets are orbiting in more or less the same orbital plane. (Planets that are highly inclined with respect to the invariant orbital plane soon become ex-planets.) The invariant orbital plane contains the center of mass of the solar system, aka solar system barycenter, and is normal to the total angular momentum of the solar system. Let that angular momentum vector define one of three principal axes; call this the z axis. You need two axes to define the invariant plane. Pick one axis that lies in this plane, call this the x axis. The third axis: It's defined by the two you already have: $\hat y = \hat z \times \hat x$.

What you want to look for is the time at which the y component of position changes signs when the x component is positive. You might want to protect against y quickly going through multiple sign changes in a short period of time.

 

[CAF123]

Hi UltrafastPED, I just released that since I did the vector and particle classes in rectangular coordinates, I am supposed to use rectangular coordinates for this question too. (Otherwise, I would have to go back and alter/redefine my vectors). To implement your method in the rectangular coordinates, the magnitude of the position vector from the central star is still $r=\sqrt{x^2+y^2+z^2}$ and I think I can express the angle as atan2(y,x) =.. (defined in http://en.wikipedia.org/wiki/Atan2). Then I need to find the time when this angle and the position from the central star are the same. I should also point out that I am just designing at the moment, coding happens later.

. Pick one axis that lies in this plane, call this the x axis. The third axis: It's defined by the two you already have: $\hat y = \hat z \times \hat x$.

What you want to look for is the time at which the y component of position changes signs when the x component is positive. You might want to protect against y quickly going through multiple sign changes in a short period of time.

Hi D H, I am not sure I understand why this would work. Could you elaborate?

 

[UltrafastPED]

You will never find "equals this value"; instead just watch for when it first becomes "greater than or equal to this value". Then you know that you have crossed that line.

 

[CAF123]

You will never find "equals this value"; instead just watch for when it first becomes "greater than or equal to this value". Then you know that you have crossed that line.

Ok, so look for the time when the position vector magnitude and angle are about the same as a point earlier in the motion. This arises because of the accuracy of the timestep right? If the timestep is too high, we might overjump the point at which they are equal.

 

[UltrafastPED]

Floating point numbers from calculations are only equal by accident ... thus if you need an actual equality you have to test a range of values centered on your "exact" value.

In your case the calculations are iterative, based on a time step for the approximation. This would require a larger interval if you were looking for "equality", approximately equal to one time step's worth of distance.

But in your case you have an angle sweeping around; all you need is to count the years, so just getting anything equal or greater will work ... depending, of course, on how your angle is generated ... is it a winding angle, or does it reset to 0 once per cycle? These are details that need to be taken into account in the actual calculations.

But from a design perspective, you just need "greater than or equal".

You will learn all about these issues if you take a course in numerical analysis ...

 

[CAF123]

Thanks. To obtain the aphelion and perihelion, in terms of design, is it sufficient to note that I want to find the point where the position vector is maximized and minimized respectively and then return the coordinates of that point in the defined coordinate system. Maybe I should say loop over all N particles, compute relative separation vector, find the magnitude, time evolve and finally return the coordinates of the point that had the max/min distance from the origin.

 

[UltrafastPED]

That sounds about right ...

 

[D H]

Hi D H, I am not sure I understand why this would work. Could you elaborate?

There are multiple meanings of the word "year". Our calendar year is approximately equal to a tropical year. What you want is something else, either the sidereal year or the anomalistic year (or both). The sidereal year is how long it takes to complete one orbit. The anomalistic year is how long it takes to go from aphelion to aphelion (or perihelion to perihelion). Using the Earth as an example, the mean sidereal year is about 20 minutes longer than the mean tropical year, and the mean anomalistic year is about 5 minutes longer than the mean anomalistic year.

Why these different measures, and why do they differ? The sidereal year measures how long it takes to complete an orbit, or revolve about the Sun by 360 degrees from the perspective of an inertial frame of reference. The anomalistic year differs from the sidereal year because of apsidal precession. This is caused mostly by gravitational influences of other planets. This apsidal precession is an effect you should be able to in your simulation if it is accurately modeling gravitation. The tropical year differs from the sidereal year because of the Earth's axial precession. I'm assuming that you aren't modeling axial precession (it's very hard to model this behavior), which is why I suggested that you want to measure the sidereal and/or anomalistic years.

What I suggested is a way to measure the sidereal year. What UltrafastPED wrote about is a way to measure the anomalistic year.

 

[Devils]

The Attempt at a Solution

So I need to find a way to test for a complete orbit.

If there is one heavy sun and 9 lighter planets, this is an n-body problem. You should be calculating the displacement vectors based on the gravitation pull of all the other objects.

What SHOULD happen is the orbits start looking elliptical, and the planets return 'near' where they started. Maybe each iteration is one Earth minute. I would probably use double precision floating point too.

Not being an astonomer (ha), I would find out the masses of the sun & planets, distance from the sun, and average year. That should be enough to get initial velocity approximations.

As suggested, test for the max and min distance from the sun to get orbit times.

Doing a 3D display would be cool too. Sounds like a nice hobby simulation. If that works add the moon!