Calculating The Speed Of A Moon?

[Meithan]

Verlet is not symplectic with 1/r potential. No known method is. Verlet is symplectic with r² potential, such as Hooke's law used in collisions, which is why it's a favorite in game dev. But no method is symplectic for all potential types.

Uhm, I may be confused then. I thought Strömer-Verlet integration was one of the prime examples of symplectic integration, e.g. see here. Do you have some references I can check out to better inform myself? I never studied this formally.

[K^2]

Uhm, I may be confused then. I thought Strömer-Verlet integration was one of the prime examples of symplectic integration, e.g. see here. Do you have some references I can check out to better inform myself? I never studied this formally.

He's kind of deriving fast and loose there. I haven't tracked it down to the exact assumptions he's making implicitly to make his claims, but I'm sure there are some.

Just take his Symplectic Euler example, and expand it for some arbitrary pair of coordinates (p, q) and (P, Q). To be Symplectic, it must approximately satisfy p_{n+1}Q_{n+1} - P_{n+1}q_{n+1} = p_n Q_n - P_n q_n. That works out for potential U ~ O(q²), which gives me error O(h²). But with U ~ O(1/q), I still have O(h) terms that I can't make disappear.

I do agree with his logic on Stormer-Verlet following from Euler. And, of course, for the special case of H(p,q) = T(p) + U(q), it reduces to Velocity Verlet. So as I've said above, for U = k|q_1 - q_2|², it is most definitely symplectic.

I'll definitely dig further into this to see if I can find where an assumption is made that's not consistent with gravity, but you can also just try it out with a two-body simulation with highly elliptic orbit to confirm that energy steadily grows regardless of h.

[Beale]

Yes, v is vx for x and vy for y. Same with acceleration and force.

Well, I've bungled it somewhere

My new add forces method, trying to implement the Velocity Verlet method.

(ax_old and ay_old are only set to 0 once, in the constructor, when you mention "maintaining" them?)

	
// Calculate force
	public void addForce(SimBody _bodyTwo)
	{
		SimBody _bodyOne = this;
		double dx = _bodyTwo._rx - _bodyOne._rx;
		double dy = _bodyTwo._ry - _bodyOne._ry;
		double dist = Math.sqrt(dx*dx + dy*dy);
		double F = (G * _bodyOne._mass * _bodyTwo._mass) / (dist*dist);
		_bodyOne._fx += F * dx / dist;
		_bodyOne._fy += F * dy / dist;

		double dt = 1000;

		double ax_new = _fx/_bodyTwo._mass;
		_vx = _vx + dt * (ax_old+ax_new)/2;
		_rx = _rx + dt * _vx + dt * dt * ax_new/2;
		ax_old = ax_new;

		double ay_new = _fy/_bodyTwo._mass;
		_vy = _vy + dt * (ay_old+ay_new)/2;
		_ry = _ry + dt * _vx + dt * dt * ay_new/2;
		ay_old = ay_new;
	}

New update method (does nothing now basically).

public void update(double _deltaTime) 
	{	
		updateTrail();
	}

20 bodies is definitely doable by direct force calculation; 300 I'd say so too, depending on what system it's running on, what programming language and how small you want the timesteps to be.

No, it's a Fortran code I wrote plotted with gnuplot.

Keep in mind the solar system is pretty stable: the orbits are nearly circular and they're pretty far apart. Close encounters are the bane of numerical N-body simulations: they'll quickly make the numerical error really bad.

If you're getting close encounters you'll probably need to make the timestep (how much each step of the calculation advances time) smaller, or even 'adaptive' (automatically adapt to be small when bodies are close and large when they are far away). But that's another story, I'd get the basic velocity verlet working first.

The "bottom" device is a 2-year old Nexus 7, with I think a ~1.2GHz processor(?), so I was not sure what number of bodies could be handled with that processing punch. Early on I ran into problem with a memory leak causing a very low framerate, which I attributed to a high-N, when in fact was not.

Fortran! I would not expect that. Although maybe lower level is best for this kind of thing.

Thankfully, as I'm trying to simulate the solar system with mostly only the planets with nearly circular orbits (The AU distance and eccentricity values pulled from wikipedia). So hopefully none of these close encounters to worry about!

Edited May 8, 2015 by Beale

[K^2]

You need to keep an ax_old and ay_old separate for each body, like you do rx, ry, vx, and vy.

Also, a few typos. a = F/m1, not m2. And you have vx in ry update.

Edited May 8, 2015 by K^2

[Beale]

You need to keep an ax_old and ay_old separate for each body, like you do rx, ry, vx, and vy.

Also, a few typos. a = F/m1, not m2. And you have vx in ry update.

Fixing those typos seems to have done the trick, I feel silly for not spotting them.

Thanks greatly for the help K^2!

My only observation is now that the inner system seems very chaotic, perhaps things are too close, or my method still contains typos.

New method.

	
public void addForce(SimBody _bodyTwo)
	{


		SimBody _bodyOne = this;
		double dx = _bodyTwo._rx - _bodyOne._rx;
		double dy = _bodyTwo._ry - _bodyOne._ry;
		double dist = Math.sqrt(dx*dx + dy*dy);
		double F = (G * _bodyOne._mass * _bodyTwo._mass) / (dist*dist);
		_bodyOne._fx += F * dx / dist;
		_bodyOne._fy += F * dy / dist;

		double dt = 1000;

		double ax_new = _fx/_bodyOne._mass;
		_vx = _vx + dt * (ax_old+ax_new)/2;
		_rx = _rx + dt * _vx + dt * dt * ax_new/2;
		ax_old = ax_new;

		double ay_new = _fy/_bodyOne._mass;
		_vy = _vy + dt * (ay_old+ay_new)/2;
		_ry = _ry + dt * _vy + dt * dt * ay_new/2;
		ay_old = ay_new;
	}

Edit:

Chaotic indeed!

[K^2]

Oh, you should be adding all of the forces together, then running a single step of Velocity Verlet on the total force. Not doing an update per source.

[Beale]

Oh, you should be adding all of the forces together, then running a single step of Velocity Verlet on the total force. Not doing an update per source.

Ah! Okay, that can explain why the smaller bodies were having seemingly much stronger effects

So:

This section of the code is spun off into a new method (I.E. the update method), where Fx and Fy are the sum of all forces?

		
double dt = 1000;

		double ax_new = _fx/_bodyOne._mass;
		_vx = _vx + dt * (ax_old+ax_new)/2;
		_rx = _rx + dt * _vx + dt * dt * ax_new/2;
		ax_old = ax_new;

		double ay_new = _fy/_bodyOne._mass;
		_vy = _vy + dt * (ay_old+ay_new)/2;
		_ry = _ry + dt * _vy + dt * dt * ay_new/2;
		ay_old = ay_new;

Edit:

Yep! That's solved it. As far as I see...

All seems to be stable after a good few many orbits!

Of course, the moon still does as it does, but this would seem to be a distance / hill-sphere problem...

Still, is good to have a more accurate form of "integration" (? Trying to get the terminology right).

Thanks for the help!

Edited May 8, 2015 by Beale

[Brotoro]

Is the scale of that diagram accurate? The Moon looks awfully far from the Earth there.

[Beale]

Is the scale of that diagram accurate? The Moon looks awfully far from the Earth there.

Yup, the scale is accurate! (At least relative to what I've defined as 'AU').

The Moon starts off at the correct distance from Earth, that picture is after a few orbits around the sun, where it flies (or more accurately drifts) away

To calculate the Moon's speed, I'm using method of:

First calculate speed for circular orbit around sun,

then add onto that speed for circular orbit around earth.

Edited May 8, 2015 by Beale

[K^2]

Yup, the scale is accurate! (At least relative to what I've defined as 'AU').

So it's not a Hill Sphere problem, then. At least, not initially.

It'd be interesting to see this with a longer trail. In fact, could you make it almost exactly a year long, and do a screen capture exactly a year after the simulation started? It could be very informative in terms of what else might be going wrong that might be possible to improve. Even if not, the orbit of Earth-Moon system around the Sun is an interesting shape to look at. They really do co-orbit the Sun, unlike moons of other planets.

As I've indicated a number of times, I expect the Moon to drift away eventually in any case, but I feel like a "few years" of clean simulation isn't an unreasonable demand here.

Have you played with the time step at all, by the way? It's more than a quarter hour long right now. Maybe that's a bit much?

[Beale]

So it's not a Hill Sphere problem, then. At least, not initially.

It'd be interesting to see this with a longer trail. In fact, could you make it almost exactly a year long, and do a screen capture exactly a year after the simulation started? It could be very informative in terms of what else might be going wrong that might be possible to improve. Even if not, the orbit of Earth-Moon system around the Sun is an interesting shape to look at. They really do co-orbit the Sun, unlike moons of other planets.

As I've indicated a number of times, I expect the Moon to drift away eventually in any case, but I feel like a "few years" of clean simulation isn't an unreasonable demand here.

Have you played with the time step at all, by the way? It's more than a quarter hour long right now. Maybe that's a bit much?

That was a good idea, the following picture is exactly (or close) to the point where the Earth makes 1 orbit.

So, you can see where the moon begins to diverge.

Initially, I thought it was because the Moon had a higher orbital speed (addition of speed for "Earth orbit"), but if my KSP knowledge has taught me anything, the Moon's orbit should be eccentric with the highest point at the "north" of that image. If I make sense?

It does feel like the Moon is barely affected by the Earth's pull at all, or at least it is hard for me to observe.

Timestep: Yeah, I've tried a timestep of 10% and 1%, and still observe the same effects (but obviously is harder to make year long or longer observations).

Agains, thanks for the help with this troublesome thing.

[K^2]

Could you post initial mass, ry, and vx for Sun, Earth, and Moon that you are using?

[Beale]

Could you post initial mass, ry, and vx for Sun, Earth, and Moon that you are using?

Of course

mass

Sun = 330000

Earth = 1

Moon = 0.01

(Could make these more accurate, e.g Moon = 0.0123, if it makes difference).

_vx

Sun = 0

Earth = 2.098771068982989E-4

Moon = 2.090426021743164E-4

_ry

Sun = 0

Earth = 500

Moon = 504

I was a little confused as to why the velocities were such small numbers, but it works.

Error in my last post, the Moon had no alterations to its velocity, and was instead just initialized at "circular" orbital speed around the sun, so it is being "boosted" outwards by Earth?

Although, the behaviour is much the same if I do give the moon extra velocity.

The bodies are initialized like this:

		_bodies[3].setDiameter(1.0f);
		_bodies[3]._texture = "t_earth";
		_bodies[3]._name = "Earth";
		_bodies[3]._mass = 1;
		_bodies[3].setAxisTilt(23);
		_bodies[3].setTrailColour(Color.CYAN);
		_bodies[3]._ry = AU;
		_bodies[3]._vx = calculateOrbitalSpeed(_bodies[0],_bodies[3],solarmass);

		_bodies[4].setDiameter(0.15f);
		_bodies[4]._texture = "t_moon";
		_bodies[4]._name = "Moon";
		_bodies[4]._mass = 0.01;
		_bodies[4].setRotationSpeed(0.01f);
		_bodies[4].setAxisTilt(15);
		_bodies[4]._ry = _bodies[3]._ry + 4;
		_bodies[4]._vx = calculateOrbitalSpeed(_bodies[0],_bodies[4],solarmass);

For a little clarification, the "calculateOrbitalSpeed" method, though I don't think there is anything wrong here.:

	
public double calculateOrbitalSpeed(SimBody _bodyA, SimBody _bodyB, double _solarMass)
	{
		double G = 6.674E-11;

		double _distance = _bodyB.distanceTo(_bodyA);

		double _velocity = Math.sqrt(G * _solarMass / _distance );

		System.out.println("" + _bodyB._texture + " " + _velocity);

		return _velocity * -1;
	}

Velocity is returned at * -1, to make them go counter-clockwise around Sun.

And finally, how I was trying to calculate the correct speed for the Moon to orbit Earth.

		
_bodies[4]._ry = _bodies[3]._ry + 4;
_bodies[4]._vx = calculateOrbitalSpeed(_bodies[0],_bodies[4],solarmass);
_bodies[4]._vx += calculateOrbitalSpeed(_bodies[3],_bodies[4],_bodies[3]._mass);

Kind regards.

Edited May 9, 2015 by Beale

[K^2]

So that vx for the moon is without adding the orbital velocity around Earth, then? Is that what you have on that screenshot? Does it still happen if you add it in?

Also, why is the Moon so far away? I'm gathering you are using light-seconds as your distance units? Moon should be slightly over 1 light second from Earth, not 4.

[Beale]

So that vx for the moon is without adding the orbital velocity around Earth, then? Is that what you have on that screenshot? Does it still happen if you add it in?

Also, why is the Moon so far away? I'm gathering you are using light-seconds as your distance units? Moon should be slightly over 1 light second from Earth, not 4.

Yes, the last picture the moon only has velocity around Sun, sorry for the confusion.

With the addition of velocity around Earth, it looks like this (After 1 orbit).

So, a little more like what to expect if Moon was not affected by Earth at all,a slightly higher aphelion with the extra speed.

The distance unit is 500th of an AU (I know is odd), I've defined an AU as 500 "distance units" so the distance unit is 1 "unit".

I hope this helps.

I realise I've been a bit irresponsible not explicitly defining what distance units I am using, but not sure if this can be fixed (or if it matters?).

Edit: With some free time this evening, I can maybe tried code in rendering for lagrange points or SOI, to see if that maybe the problem.

Edit more: If I initialize the Moon very close to the Earth, and if the Moon gets initialized with no extra velocity, then it begins to drift away, before stopping and falling back towards Earth for a collision. Not sure how helpful this is, but at least puts my mind to rest that the Earth exerts some influence, despite appearances.

Edited May 9, 2015 by Beale

[K^2]

1/500th of AU is almost exactly 1 light second. The distance to the Moon is still much smaller. The semi-major axis is 0.00257AU, which is 1.285 of your units. But it still shouldn't have that effect.

You do iterate through all possible pairs for which forces are applied, right?

Sorry if I'm bugging you with too many questions, but your results for the Moon are a bit puzzling, given how well everything else seems to work.

[Beale]

1/500th of AU is almost exactly 1 light second. The distance to the Moon is still much smaller. The semi-major axis is 0.00257AU, which is 1.285 of your units. But it still shouldn't have that effect.

You do iterate through all possible pairs for which forces are applied, right?

Sorry if I'm bugging you with too many questions, but your results for the Moon are a bit puzzling, given how well everything else seems to work.

Bugging? No no! You are asking all the right questions

I've realised the rendered diameter of the planets are way off, taking 500th of an AU the current Earth is rendered at something close to 300,000 km in diameter.

Which explains why the position of the moon "looked right" because the two bodies appeared much bigger.

Every possible pair should be being calculated, but I'll admit it's hard to tell when my understanding of Barnes-Hut is quite weak.

I will try place the Moon at correct distance, switch to a direct method rather than barnes-hut and see what happens.

Kind regards.

[K^2]

Yeah. Given the size of the system, all you need to do is loop over all bodies, then loop over all other bodies again inside that loop to add forces, making sure to skip "pairs" where indices match. Then run separate loop to update all velocities and positions. You want that in separate loop to minimize errors.

[Beale]

Yeah. Given the size of the system, all you need to do is loop over all bodies, then loop over all other bodies again inside that loop to add forces, making sure to skip "pairs" where indices match. Then run separate loop to update all velocities and positions. You want that in separate loop to minimize errors.

Switched to direct method and...

It works! It bloody works!

Without altering the Moons distance at all, it comfortably entered (an eliptical) orbit around the Earth, which was quite stable for a few "months".

So, there seems to be the problem, my implementation of Barnes-hut is bunk.

Method

			
// Calculate forces for each body.
			for (int i = 0; i < N; i++) 
			{			
				_bodies[i].resetForce();
				for (int j = 0; j < N; j++) 
				{	
					if(i!=j)
					{
						_bodies[i].addForce(_bodies[j]);
					}
				}
			}	
			// Update Bodies
			for (int i = 0; i < N; i++) 
			{	
				_bodies[i].update(_deltaTime);
			}

Many thanks for your help!

Edit: Moving the Moon closer produces a circular, pretty damn stable orbit, fantastic!

Edited May 10, 2015 by Beale

[Meithan]

Switched to direct method and...

It works! It bloody works!

Oh, I'm glad you figured it out! So you'll have to double check your Barnes-Hut, although it's really not necessary for this few bodies (and in fact direct force calculation is more precise).

I was about to post a couple experiments I just did to see if the Moon remains in orbit in my old simulations. Seems it does:

If we zoom in a bit it's clear the Moon is actually orbiting the Earth, since it follows the Earth in its orbit around the Sun (note how it crosses the Earth's orbit back and forth):

Next, I tried to see how stable the Earth-Moon system is for longer periods of time, so I computed 1000 years of evolution. The result: rock stable.

There's no point showing the orbits because they just look almost exactly as the first picture above. It's more enlightening to plot the Earth-Moon distance as a function of time:

The vertical axis is the center-to-center separation between the two in thousands of kilometers (should be around 380,000 km). There are some spikes in the distance here and there, but it does not seem that the Moon is drifting closer or farther in the simulation.

Edited May 10, 2015 by Meithan

[Beale]

Oh, I'm glad you figured it out! So you'll have to double check your Barnes-Hut, although it's really not necessary for this few bodies (and in fact direct force calculation is more precise).

I was about to post a couple experiments I just did to see if the Moon remains in orbit in my old simulations. Seems it does:

http://i.imgur.com/KRmBb4Wm.png

If we zoom in a bit it's clear the Moon is actually orbiting the Earth, since it follows the Earth in its orbit around the Sun (note how it crosses the Earth's orbit back and forth):

http://i.imgur.com/PckYaJbm.png

Next, I tried to see how stable the Earth-Moon system is for longer periods of time, so I computed 1000 years of evolution. The result: rock stable.

There's no point showing the orbits because they just look almost exactly as the first picture above. It's more enlightening to plot the Earth-Moon distance as a function of time:

http://i.imgur.com/h2nH9YAm.png

The vertical axis is the center-to-center separation between the two in thousands of kilometers (should be around 380,000 km). There are some spikes in the distance here and there, but it does not seem that the Moon is drifting closer or farther in the simulation.

Thanks

I'll have to take another look yeah, but as you say there seems to be no measurable performance hit on the small system.

For now, the Barnes-Hut mode will be locked away in an options menu (Along with forward Euler).

In fact, here is a recording of the Moon orbit across a year, stable, no?

On your own simulation, colour me impressed! 1000 years seems really amazing to me

[Meithan]

In fact, here is a recording of the Moon orbit across a year, stable, no?

On your own simulation, colour me impressed! 1000 years seems really amazing to me

Colour me impressed, that's a very nice looking app! What language/platform is that written in?

[Beale]

Colour me impressed, that's a very nice looking app! What language/platform is that written in?

Ah, great appreciate it, many thanks

The App is written in Java, Eclipse, I use the libGDX framework and API for graphics rendering, loading.

libGDX very nice, allows really easy cross-platform development, enable me to test on Android and desktop (Windows/osx/unix) simultaneously!

It does need explicit asset management though, which is a little alien when I was used to Java's garbage collector. But, the control has fostered a bit of an interest in languages without a GC!

Small extra question if you know: how comparable we're the perturbations (?) In your simulated Moon's orbit to the perturbations of the real Moon's orbit, does it not deviate a distance like that IRL? I should brush up on my basic solar system knowledge eheheh.

Edited May 10, 2015 by Beale

[Meithan]

Ah, great appreciate it, many thanks

The App is written in Java, Eclipse, I use the libGDX framework and API for graphics rendering, loading.

libGDX very nice, allows really easy cross-platform development, enable me to test on Android and desktop (Windows/osx/unix) simultaneously!

It does need explicit asset management though, which is a little alien when I was used to Java's garbage collector. But, the control has fostered a bit of an interest in languages without a GC!

And this is meant to run on a mobile? Very nice. I've done my share of programming in my life but I've never wrote code for mobile, and frankly I've been wanting to learn. Do you think you could share your code? Just for personal, educational use (for me)?

Small extra question if you know: how comparable we're the perturbations (?) In your simulated Moon's orbit to the perturbations of the real Moon's orbit, does it not deviate a distance like that IRL? I should brush up on my basic solar system knowledge eheheh.

Nah, that was just a toy simulation, it's not meant to be a faithful representation of reality. The Moon's orbit is slightly eccentric: 363,000 x 405,000 km. So the main explanation for the variations you see is just the Moon's eccentricity. The rest is numerical error. In reality, there are also several perturbations that make things a little more complicated: tidal forces due to the Earth and the Sun's gravitational attraction are probably the largest perturbations.

Fun fact: the Sun exerts twice as much gravitational force on the Moon as the Earth does. Go ahead and compute it if you want to check. It's the only natural satellite in the solar system for which that's true.

[K^2]

Hence, some people pushing to re-classify Earth-Moon system as a double planet.