Trajectories to the Moon

[Not Myself]

Following up from another thread ...

I wrote a C program, 199 lines, but a lot of white space and declarations - very little in the way of actual calculation.  It uses the dumbest algorithm possible.  Calculate acceleration due to gravity (earth and moon), and next period's position and velocity are this period's position and velocity plus velocity and acceleration (times the length of the time step).

Started with a stationary earth, and a moon in a circular orbit, Newtonian gravity, no other objects in the universe, and uniform mass distribution within the two.  Start with an initial (x,y,z) position, and an initial velocity vector, and an initial time (needed to make sure you get to the moon, not where the moon was two weeks ago).  I started at an altitude (from the centre of the earth) of 6,400 km, within the plane of the lunar orbit.

One of my first trajectories got within 6,700 km of the centre of the moon, and then blew past it.  Changing the time step from 0.1 seconds to 0.01 seconds made a difference of about 100 km!  Moving from 0.01 to 0.001 made a difference of about 1 km, but also made it so slow that I didn't want to do it.  I need to improve my approximation algorithm, or use adaptive time steps (small time steps when close to the earth or the moon, big time steps when far away).

Toying with it, trying to find a good orbit that loops around the moon and then drops back, is challenging by trial-and-error (the method I am using).  The orbit is highly sensitive even to a 1 m/s initial velocity change!  I also need to develop better diagnostic/visualisation tools - right now, I'm just looking at a text output of a bunch of numbers, and try to adjust the initial starting time and velocity based on that.

So my conclusion so far is, even with the computers we have now, you still need to be at least a little intelligent in your choice of algorithm.

[Not Myself]

I might have misstated the time steps by a factor of ten, I'm not sure.  Didn't record all the early trajectories.

Also trying to think if there is an easy way to check for correctness.

ETA - never mind, was right the first time.

Also, adding an accel*(dt^2)/2 term to the update of position helps a lot, but not as much as I had hoped.

[Bob B.]

Following up from another thread ...

I wrote a C program, 199 lines, but a lot of white space and declarations - very little in the way of actual calculation.  It uses the dumbest algorithm possible.  Calculate acceleration due to gravity (earth and moon), and next period's position and velocity are this period's position and velocity plus velocity and acceleration (times the length of the time step).

Started with a stationary earth, and a moon in a circular orbit, Newtonian gravity, no other objects in the universe, and uniform mass distribution within the two.  Start with an initial (x,y,z) position, and an initial velocity vector, and an initial time (needed to make sure you get to the moon, not where the moon was two weeks ago).  I started at an altitude (from the centre of the earth) of 6,400 km, within the plane of the lunar orbit.

Your method sounds  much like what I did to create my trajectory animations.  However, I did it with an Excel spreadsheet, which I've discovered actually works pretty well for this type of thing.

Referring to the sentence bolded above, I have a suggestion that will improve your accuracy.  If I understand you correctly, you're doing the following for each time step:

v_f = v_o + a_o * dt

d_f = d_o + v_o * dt

Instead try something like this (if C doesn't use For-Next loops, use the C equivalent):

a_f = a_o

For i = 1 to x

v_f = v_o + (a_o + a_f)/2 * dt

d_f = d_o + (v_o + v_f)/2 * dt

a_f = calculate updated value based on new position and time

Next i

Using x = 3 or 4 is probably adequate to significantly improve your results, though more is better.  The accuracy that required a very small time step using your method should be attainable with a larger time step using my method, and processing time should be faster (at least I think so).  I accomplish the same thing in Excel by enabling iteration.

[Bob B.]

The following questions were asked in the "too much education" thread, but I think it is a better place to discuss it.

Circumlunar Free Return Trajectory

How did you work that out?  The earth is fixed and the moon is easy under the assumptions made.  For the spacecraft, I take it it was numeric integration?

I think you pretty much answered this one for yourself.  I used a method very similar to the one you are using, except I calculate the velocity and position at the end of each time period using the average acceleration and velocity over the time period.  This is done by iteration, which I demonstrated in my last post.

Was the amount/direction of the initial burn determined by trial and error, or is there a better way to do it?

Trial and error, though it's possible to calculate approximate starting conditions.  Once you've done two or three runs, you can interpolate from the results to zero in on the desired trajectory.  It is then just a matter of tweaking until you get it just right. 

[Not Myself]

If I understand you correctly, you're doing the following for each time step:

v_f = v_o + a_o * dt

d_f = d_o + v_o * dt

That is what I did at first, but I later added a (1/2)*a₀*(dt²) term to the distance, which did lead to a substantial improvement.  The method you recommend would seem to eliminate the need for an acceleration term, as it explicitly estimates the average velocity over the entire time step.

Instead try something like this (if C doesn't use For-Next loops, use the C equivalent):

  C has, let's say, a rather "sparse" notation, that can appear cryptic and perhaps unintuitive if you're not used to it, but which is rather fast to type once you get the hang of it.  The notation would be

Code: [Select]

for(i=0;i<N;i++) {<code to be executed N times here>}

a_f = a_o

For i = 1 to x

v_f = v_o + (a_o + a_f)/2 * dt

d_f = d_o + (v_o + v_f)/2 * dt

a_f = calculate updated value based on new position and time

Next i

Using x = 3 or 4 is probably adequate to significantly improve your results, though more is better.  The accuracy that required a very small time step using your method should be attainable with a larger time step using my method, and processing time should be faster (at least I think so).  I accomplish the same thing in Excel by enabling iteration.

OK, I will have to try this.  I was slightly confused at first, because I thought you were dividing the time step into substeps - I didn't realise the iterations were successively better approximations over the entire time step at first.  But I will give it a whirl.

Another thing I am thinking about is having the size of the time step depend on how "curved" the gravitational field is.  The idea is, use small time steps when close to the earth or moon, and small changes in position could have potentially large changes in gravitational acceleration, but large time steps when out in the countryside, far away from both earth and moon.

I also need to look up the methods cjameshuff suggested in other thread.

Many thanks for help!

[Bob B.]

I didn't realise the iterations were successively better approximations over the entire time step at first.  But I will give it a whirl.

I found this to give much better results than anything else I've tried.

Another thing I am thinking about is having the size of the time step depend on how "curved" the gravitational field is.  The idea is, use small time steps when close to the earth or moon, and small changes in position could have potentially large changes in gravitational acceleration, but large time steps when out in the countryside, far away from both earth and moon.

I've done the same thing in the past.  I've generally linked the time step to the velocity, e.g. v<=x then dt=0.1, x<v<=y then dt=1.0, etc.  The theory being, the slower you're moving, the less change there is over the time period, and the less sensitive your calculations are to a bigger time step.

[cjameshuff]

Your method sounds  much like what I did to create my trajectory animations.  However, I did it with an Excel spreadsheet, which I've discovered actually works pretty well for this type of thing.

I'd strongly discourage using Excel for anything vaguely numeric (which, yes, covers most things you'd want to use Excel for). You simply can't rely on Excel doing anything like what you think you're telling it to do, without making invalid simplifications, losing precision, or just plain giving incorrect output.

For something better in just about every way (though oriented more towards statistics), look at R: http://robjhyndman.com/researchtips/rvsexcel/
http://www.burns-stat.com/pages/Tutor/spreadsheet_addiction.html

As for timesteps...keep in mind that floating point computer arithmetic is essentially binary scientific notation. You have a fixed number of bits multiplied by an overall scale factor represented as a power of 2. Gravitational accelerations in this sort of simulation tend to be quite small...0.006 m/s^2 at Earth's distance from the sun, while velocities tend to be relatively high...30000 m/s for Earth. Think of aligning the decimal points when adding numbers with a fixed number of digits, and rounding the result to the same number of digits. When adding the small number to the large one, you end up losing the least significant bits of the small one.

Reducing the timestep makes the difference in scales worse. 30 m/ms and 0.000000006 m/ms^2, a t^2 factor of 0.000001 s instead of 1 s, etc. You can easily get into situations where small numbers (like those based on acceleration) lose almost all their precision in your calculations, or even round to zero. And even when that doesn't happen, this roundoff error accumulates with each timestep, making things worse...smaller timesteps means more timesteps. So it can really pay off to use an algorithm that converges faster and allows you to use larger timesteps.

[Not Myself]

Agreed that Excel has issues, but I think Bob B. is (and I know I am) just having a laugh, not actually sending spacecraft to the moon

The big advantage I can see to Excel here is the ease of use.  For example, if I were to work it up to plot out a trajectory, I could define various metrics on the trajectory (how close to the moon does it come, how long does it take, etc.), and then use built-in tools to search for a good trajectory.  To do something like that with my C program would be a bit of work.  In other applications, I have definitely used Excel to come up with a first cut at the solution, then refined it on other software platforms.

I do know of one built-in worksheet function where Excel is just flat out wrong.  Not close but slightly inaccurate; just wrong.

I once taught an "application of statistics" kind of course, and my advice on software tools (including, but not limited to Excel) was to make sure the built-in functions were doing what you thought they were doing.  Quite a few people violated that principle

[Glom]

Off the starting mark.

Doing this in an object oriented programming language is tricky.

Curse Java's lack of operator overriding!  Yes, that's right, I'm doing it in Java because I haven't got the danged C++ to work on this machine.

[Not Myself]

OK, this is one of these things that looks easy, until you actually try it

Some comparisons.  First, the Coelacanth-A method, which is described by Bob B. as:

v_f = v_o + a_o * dt

d_f = d_o + v_o * dt

The table shows the position of a particular trajectory 144 hours after launch.

Code: [Select]

Time step    Distance to earth    Distance to moon
  1.0  sec     331955.4 km          88455.9 km
  0.1  sec     400272.7 km          95209.0 km
  0.01 sec     391201.0 km          95705.9 km
  0.001 sec    390374.0 km          95771.5 km
  0.0001 sec   390292.0 km          95778.2 km

As shown, the Coelacanth-A method really takes small time steps before it works well - even going from 1 ms to 0.1 ms makes a difference of almost 100 km in earth distance after 144 hours.  Presumably, if the time step were small enough, we would run into the problem cjameshuff describes, in which precision would be lost, and in the extreme case, the acceleration would round to zero when added to the velocity, and forecast trajectory would be a straight line.  It is not obvious that anything like this is happening, because the results do seem to be converging, although I don't know any simple way to check that they are converging to the truth.  Also, at least with my implementation and my computer, the 0.1 ms time step calculations are so slow they're painful.

The Coelacanth-B method is similar to the Coelacanth-A method, but includes an acceleration term in the distance update.

v_f = v₀ + a₀ * dt

d_f = d₀ + v₀ * dt+a₀*(dt*dt/2)

Same table for the Coelacanth-B method.

Code: [Select]

Time step    Distance to earth    Distance to moon
  1.0  sec     479512.6 km          99207.1 km
  0.1  sec     395052.5 km          95464.4 km
  0.01 sec     390741.1 km          95743.7 km
  0.001 sec    390328.5 km          95775.4 km
  0.0001 sec   390292.0 km          95778.2 km

The two methods seem to agree rather closely at 0.1 ms time steps, but the Coelacanth-B method seems to get there faster than the Coelacanth-A method (which would be expected).

Then there is the Bob B. method.

a_f = a_o

For i = 1 to x

v_f = v_o + (a_o + a_f)/2 * dt

d_f = d_o + (v_o + v_f)/2 * dt

a_f = calculate updated value based on new position and time

Next i

The Bob B. method will be indexed by number to indicate the number of iterations.  The Bob B.-1 method is equivalent to the Coelacanth-B method.  The Bob B.-2 method results appear below.  I didn't go to 0.1 ms, because it is so painfully slow, and the method appears to have converged sooner.

Code: [Select]

Time step    Distance to earth    Distance to moon
  1.0  sec     390315.1 km          95776.2 km
  0.1  sec     390283.2 km          95779.0 km
  0.01 sec     390282.9 km          95779.0 km
  0.001 sec    390282.9 km          95779.0 km

The Bob B.-3 method results appear below.

Code: [Select]

Time step    Distance to earth    Distance to moon
  1.0  sec     390315.1 km          95776.2 km
  0.1  sec     390283.2 km          95779.0 km
  0.01 sec     390282.9 km          95779.0 km
  0.001 sec    390282.9 km          95779.0 km

The Bob B.-4 method results appear below.

Code: [Select]

Time step    Distance to earth    Distance to moon
  1.0  sec     390315.1 km         95776.2 km
  0.1  sec     390283.2 km          95779.0 km
  0.01 sec     390282.9 km          95779.0 km
  0.001 sec    390282.9 km          95779.0 km

So, at least for this trajectory, the Bob B.-N method appears to strongly outperform the Coelacanth-B method for N>1, but there is little to no additional benefit for taking N>2.

Future possible improvements:

a) Make the time step depend on the spacecraft velocity and curvature of the gravitational field.  Based on my suspicion that most of the approximation error takes place while close to earth and close to the moon, but this is as yet unconfirmed.

b) Create a dynamic Bob B. method, in which the number of iterations is variable, depending on how much the spacecraft position and velocity are changing with each iteration.  However, if the trajectory tried here is typical, there would be little incremental benefit to more than two iterations, so it is not obvious that there is much benefit to be had from this strategy.

c) Use totally different methods, like those suggested by cjameshuff.

The idea that any old algorithm can be used, and the approximation error overcome by brute computing force, does not seem to have been confirmed here.

[Not Myself]

Curse Java's lack of operator overriding!

Yes, in C, I have functions called things like "vector_add", "vector_sub", "vector_scale", and "vector_len"   I need to relearn C++.

[Bob B.]

I'd strongly discourage using Excel for anything vaguely numeric (which, yes, covers most things you'd want to use Excel for). You simply can't rely on Excel doing anything like what you think you're telling it to do, without making invalid simplifications, losing precision, or just plain giving incorrect output.

I don't use Excel for anything requiring the kind of precision you're talking about.  Excel has always served me well and produced adequate results for what I'm using it for.  Also the ease with which I can produce graphs makes Excel ideal for my purposes.

[ka9q]

I've been programming a lot in C since 1977 or so but despite the popularity of C++, for a long time I never felt much urgency in learning it. I mostly did operating systems and network protocols where C was (and mostly remains) the language of choice.

But I'd been intrigued by operator overloading ever since C++ came out, so I decided I'd give it a try as soon as I found an application that could really use it.

That application was Reed-Solomon error correction. Like many (most) block-type error correction codes, the Reed-Solomon code is an "algebraic" code based on the properties of Galois (finite field) algebra. Basically, you construct codewords that are overdetermined systems of Galois polynomials, and the redundancy lets you detect and solve for a certain number of missing or erroneous coefficients.

I had already implemented RS and other error correcting codes in C, but the resulting program tends to be hard to follow because the algorithms are lost in the details of implementing basic Galois operations like multiplying, dividing and evalulating polynomials.

So it seemed natural to first write a set of classes to implement the primitive Galois field operations and then implement a RS codec by writing the formulas just as they appear in the textbooks. People should be able to read my program source and see the underlying algorithms uncluttered by their implementation details in ordinary C.

I succeeded in that my C++ RS codec produced correct results and the actual codec source was compact and easily understood by anyone with the textbooks. But I was never able to get it to execute faster than about 25% of my existing RS codec in C. Even with all the usual optimization tricks I still spent a lot of time creating, destroying and copying C++ objects.

It's been a very long time since speed was the overriding consideration in almost everything having to do with computers, including the choice of  programming language. But I don't think we'll ever reach the day that speed becomes completely unimportant. So I continue to have very mixed feelings about C++.

Way back in the early 1980s I wrote a general purpose orbit modeling package in C. More recently I have been slowly writing a new one in C++ that uses operator overloading to implement basic vector operations and the like. I suspect it will also be much slower than my existing C version. But I figure that writing at a higher level of abstraction, without the continual clutter of implementing these basic operations in C, will let me achieve correct operation more readily when I first implement some new and complex algorithm. If I then need more speed, I can use my correct but slow C++ version as a verification reference for an optimized C version.

[Glom]

Before I go any further with the "science" I'm taking a break to learn some 3D as I need to visualise what it's doing so far.  In particular, I want to make sure I got my signs the right way round.  I was about to start putting in my own projection algorithms, before thinking that's not the smartest way to do things.

My programme is a bit more of an n-body simulator, or at least most of the classes are designed for n-body simulation.  Obviously the main class will be specific to this case where there are 3 bodies, Earth, Moon and satellite.

Huh, whadya know, Java already has a 3D vector in the javax.vecmath package.

[ka9q]

Hi, I know this is an old thread but it's still a good one. I found it when searching for previous discussions about lunar trajectories.

A couple of suggestions about computing trajectories. While I never got around to finishing it, some years ago I looked at what it would take to do a good N-body trajectory predictor. I learned just how hairy numerical integration can be, how many different ways there are to do it, and how each tends to be optimal for different kinds of problems.

One of the most popular numerical integration methods is the 4th order Runge-Kutta, and I see a lot of people using it with pretty good results. The main complaint in the textbooks seems to be that it's slow, but these are also old textbooks.

It's my understanding that in the real world, people rarely just integrate all the forces. Most real trajectory problems are close approximations to 2-body problems with perturbing forces that, while very significant, are still small compared to the force from the central body. This led to something called Encke's method, whereby you use ordinary, simple, Keplerian equations to solve the 2-body problem and then numerically integrate all the other forces to keep track of their effect. When the effect of the perturbations gets large enough, you update the Keplerian elements to bring the error back to zero, and continue.

This apparently gets very good results with minimal computer time. It should work well on lunar trajectories because there's only a short time when earth and lunar gravity are comparable. Most of the time you're orbiting either the earth or the moon.