View Full Version : Create an Orbital Simulator
I am toying with the idea of creating a general orbital simulator. Possibly lacking the proper terminology, here is a basic description of my goal, make that wild dream.
Step 1. Working in 2 dimensions, user enters the mass of 1 to N object, their initial position, their initial velocity, and a time calculation interval. For each time interval, the simulator calculates the forces that each body applies to each other body, the summed vector of that force, applies it to the initial velocity, and determine the position at the end of the time interval. It then shows the position of each body. Repeat.
Step 2. Add niceties like lines showing where each body has been for some number of iterations so we can better see the activity.
Some future step, make it three dimensional.
This seems a rather ambitious project so before I start:
Do you know of any simulators of this nature?
Do you have any advice for a novice in celestial mechanics?
Do you think this would be better to attempt on Microsoft Windows or under Linux, I have both.
Thank you.
Nabeshin
May2-08, 08:52 PM
I contemplated doing this but found it to be beyond my current physics (and probably computer programming) abilities. I'd be really interesting to know the method you use to accomplish all this!
Step 1. Working in 2 dimensions, user enters the mass of 1 to N object, their initial position, their initial velocity, and a time calculation interval. For each time interval, the simulator calculates the forces that each body applies to each other body, the summed vector of that force, applies it to the initial velocity, and determine the position at the end of the time interval. It then shows the position of each body. Repeat.
This approach will work fine with one object. The object will move in a straight line because there no forces act on it. This approach will not work very well with two or more objects. The approach you outlined is called the Euler method for solving an initial value problem (i.e., a system whose complete system state is known at some point and whose evolution is described by a set of ordinary differential equations). The Euler method is very simple to implement but unfortunately does an incredibly lousy job except for the simplest of systems.
Fortunately, a lot of numerical integration techniques exist that do yield good accuracy. They are all more complex than the Euler method, but many are only a bit more complex. A couple of nice balances between complexity and accuracy are the velocity verlet technique, which is widely favored in chemistry simulations, and the fourth-order Runge-Kutta (RK4) technique, which is widely favored in orbital simulations. A quick google on velocity verlet and Runge-Kutta will yield a lot of information.
Hello D H,
We quickly delve into the reasons for my post. I don’t understand why the use of the simple equations f = ma, and f = ( G * M1 * M2 ) / r squared will not work. Comparing two bodies, there is a force between them that will change their velocity. If the masses and velocities are anywhere close to “right” then the two masses should orbit. Assume I iterate fast enough to get sufficient resolution but slow enough to prevent thrashing with no results. I will probably start with a resolution of one hour, 8,760 iterations will simulate a year on the scale of our solar system.
I looked at Runge-Kutta on Wiki and do not understand it at all. I am not trying to predict or smooth our trajectories, just create an instance and calculate the values for the next instance.
My first hurdle will be how to express mass in an equation. I will be looking that up.
Thanks for taking the time to reply.
LoneDragon
May3-08, 09:31 AM
Both bkelly and DH are correct. If you use small timesteps (but very very small ones), the basic f=ma for N^2 interactions, will produce a proper simulation. But if two objects in that N body simulation *ever* approach closely enough so that the time-space resolution causes "thrashing" or a large error term in the step wise calculation, the simulation will only be characteristically correct, but will diverge greatly from reality. For the solar system, at a one hour resolution, you should have no problems for the solar system.
Runge-Kutta (RK) methods do not calculate smoothed out orbital calculations, or predictions, per-se. Rather, the multiple differnce terms you see in Runge-Kutta calculations, are for calculaing a quadratically accurate polynomial *approximation* of the orbit. Euler only calculates a first order polynomial *approximation* of the orbit.
Please read about Taylor series, rectangular integration rule, simpson's rule of integration, Euler method, and RK method, in detail and trials in code, to understand how Euler method, is like a first order rectangular integration, and RK on the fourth order method, is somewhat-like simpson's rule that calculates mini-chained parabolic integrals each step of the second order. It all has to do with the infinite Taylor series, and error terms of asymptotic convergent functions like gravitational fields. With RK, you could make the steps 3-5 days long, and be far more accurate with Mercury and close approaching comets, than you would using the Euler method. RK even comes in higher order methods that let you further increase the step size, as each step is even higher order accurate, with an exponentially shrinking residue in the Taylor series like infinite series. Trust me, if you learn RK, you will never regret it. Numerical Recipies in C, by Press, Teukolsky, et al, is a good book with such algoritms available in ready to use C code for either Linux or Microsoft C compiler environments.
Another method you may want to consider, if RK is too complex, is to use Burlish-Stoer (BS) methods. In short, first, you can take Euler method at 16 hours*1 step, 8*2, 4*4, 2*8, and 1 hours*16 step. Then do a least squares fitting of the 5 exponential refined ranges to extrapolate the asymptotic-convergent value using an appropriate asymptotic function, which you can vary it's parameters, until it fits the projected 16, 8, 4, 2, 1 series to find the convergent value at "0 hours". Then take that convergent value and take the 16 hour step at the infinitessimally projected 16 hour time step, using the asymptotic (0 hours*infinite step) projected step with approximated "perfect" precision. You will have to find the right equation for this application for the least squares fitting of the exponential approximations, and make sure the exponential ranged projections are not oscilatory, but asymptotic, so that the curve fitting is always accurate each step.
Incidentally, make sure any Euler, RK, or BS curve fit and planetary position and velocities are all done in double precision, and even then, that your numbers retain large precision in the LSB of the numbers. Study the IEEE (754) format for float and the one for doubles, paying special attention to the mantissa for precision checking. The Euler method, especially, if the time steps are too small, may become susceptible to quantization errors in extreme distance or closeness calculations, without checking precision in the mantissa. RK is the best, and BS can be used intermediately.
Step 1. Working in 2 dimensions, user enters the mass of 1 to N object, their initial position, their initial velocity, and a time calculation interval. For each time interval, the simulator calculates the forces that each body applies to each other body, the summed vector of that force, applies it to the initial velocity, and determine the position at the end of the time interval. It then shows the position of each body. Repeat.
I wrote an 3D, N-body orbital simulator, with displayed output, in basically the same manner you outline. Variables were updated, for the most part, in parallel. The model suffered orbits that would gradually increase in radius over iteration. I repaired it by breaking the parallel-set rule. Near collisions introduced a problem that I never did fix before loosing interest. A close approach would throw the two bodies out towards infinity. I attempted a conservation-of-energy kludge to fix the problem. It was a bust.
Wow, that was a lot to say. This will be more complicated that I imagined, even considering that I knew it would be more complicated than I knew. (Does anyone know of a happy face icon with a recursive grin?)
Simulating a solar system would be nice, but what I really want is the ability to start up a simulation with 3 or more stars in orbit around each other and see what happens. As noted by LoneDragon, that will involve some close approaches, high velocities, and obviously some time-space alterations (as compared to our perceptions of what is normal).
After reading Phrak’s post, I had best expect some anomalies in the behavior.
My thought is to use C++ or C# and create a class for celestial objects. For each body in a simulation, instantiate an object giving it the requisite mass, size, location, and velocity. (and other things I don’t know about yet.) The class would have a method that takes as its input parameter a pointer to another object. That method would then compare the parameters of the two bodies and calculate the resultant force at that point in time. In doing so, it would calculate the position and velocity for some specified time in the future. That time would be variable and should probably decrease based on distance, velocity, and I don’t know what else.
It seems like if I start it simple without close approaches, then when I get the mechanics of being able to calculate the resultant vectors of two bodies, and displaying an image of the bodies, then I can refine the equation within the class to any degree I want. Obvious performance will suffer, but that is the way it goes.
That reminds me. Back in the 1980s I created a program to display views of the Mandelbrot set. (This was when the standard speed was 4 MHz and a fast machine was 8 MHz.) I quickly had to use the 8087 for precision and speed. Using Turbo Pascal, I then dropped into assembly language to program the 8087 for speed and more accuracy. Finally when that was not enough resolution, I wrote some 128 bit fixed point arithmetic routines in assembly language and ran that. I never did find the threads between what appeared to be islands of points within the set.
But I digress.
I have seen images from simulations for galaxies, but I don’t know of any right now. It would be interesting to know how they ran those simulations, on what machines, what language, and how long it tool. A bit of searching is in my immediate future.
Thanks for your replies.
Both bkelly and DH are correct. If you use small timesteps (but very very small ones), the basic f=ma for N^2 interactions, will produce a proper simulation.
No, it won't. I work in this field and we have very explicit rules barring the use of the Euler method. The basic problem is that the Euler method conserves neither angular momentum nor energy. (RK4 has this problem too, but the error is vastly diminished compared to the Euler method.)
When writing any simulator, it is always a good idea to verify and validate the resultant simulation. One way to do this is to apply simulation to a case where the correct answer can be computed analytically. The case of 2 bodies has an analytical solution in Kepler's laws. The simplest case is two bodies in circular orbits about each other. The problem with the Euler method becomes apparent in this simple simulation: The bodies spiral away from one another. One can reduce the error by taking very small steps. Suppose you place two bodies in circular orbits about each other, simulate for what should be ten orbits, and examine the relative position error. To get a relative error of 1% or smaller, an Euler integrator would have to divide that 10 orbit time interval into well over 3 million subintervals. In comparison, an RK4 integrator would only need to divide that 10 orbit time interval into about 400 subintervals. Achieving a relative error of less than 0.1% requires about 900 subintervals for RK4 and is simply unachievable with Euler integration. The problem is that as one makes the time steps smaller and smaller, at some point the error starts to grow because (1+10-16)-1 = 0 using the ordinary double-precision arithmetic available on almost all computers now. Make the steps too small and the state freezes.
Edited to add:
RK4 is not the be-all and end-all of numerical integration techniques. It is the most widely-used numerical integration technique because it is easy to implement and offers fairly good accuracy. There are more complex Runge-Kutta techniques that do a much better job than RK4. For example, a Runge-Kutta-Feldberg 7/8 integrator achieves better than 10-12 relative accuracy for the 10 orbit scenario described above by dividing the 10 orbit interval into about 1000 intervals. In comparison, a RK4 integrator run with the same time step will achieve a relative accuracy of about 5*10-5. The relative error with the Euler method at this coarse step size is greater than one. There are other more advanced techniques that do much, much better than any Runge-Kutta integrator.
what LD doesn't tell you is that all n-body sumulations are only characteristically correct. how correct you want yours depends upon your intent in doing the simulation.
If you want some initial results you don't need polynomial approximations.
btw, a break-down at close approach is characterised by the displacements obtained between iterations per the radial displacement between objects (rather than the radius you assign the objects).
LoneDragon
May3-08, 12:54 PM
what LD doesn't tell you is that all n-body sumulations are only characteristically correct. how correct you want yours depends upon your intent in doing the simulation.
If you want some initial results you don't need polynomial approximations.
btw, a break-down at close approach is characterised by the displacements obtained between iterations per the radial displacement between objects (rather than the radius you assign the objects).
Actually, I did. All taylor series family like methods are asymptotic approximations, with an error residue. More than once I referred to approximation, error residue, and close approach gross errors.
LD: It all has to do with the infinite Taylor series, and error terms of asymptotic convergent functions " sums all of this the best, and even higher order terms like RK4 or RK6 or RK... predictor corrector methods are all approximations. My one "proper" simulation remark, only refers to a simulation like the solar system where things don't get too close to each other and are near circular orbits, excluding moons. A complex stellar environment, with many close approaches, agreed, does require the more precise methods, like RK4 or RK6 and so forth. All of them do not preserve everything, if you want to be technical, as nearly all natural three-body problems are analytically unsolveable, and N body too.
For bkelly,
For an approximation-near-analytic text, to test a simulation of yours, there are astronomical ephemerides, with heliocentric velocity and position information. JPL has online ephemerides which may also be used to test your method, before doing complex N body stellar simulations, with close approaches, which will require higher order terms, near-accuracy methods.
Wow. This is beginning to sound like an incredibly complex task. I don't know enough to really understand what you guys are saying, but I will continue looking to determine the odds of getting something running that remotely resembles what would happen in reality.
Going to www dot ask dot com (this thing seems to tell me I cannot post a URL so I "said" it instead) and asking what are ephemerides yeilds nothing useful. A google search tells me that they give the location of satellites for each day of the year. Haveing been a systems engineer tracking rocket launches and using tracking antennas, I know what a TLE is. The JPS site appears beyond me, but I will wander around there and see what I can understand.
However, to my knowledge, I don't want to predict anything. I want to run a simulation that takes steps based on the current conditions and calculates movements to the next position or iteration.
However, again, as I alluded earlier, I am vaguely aware of the severe limits of my knowledge in this field. The core question has become:
Concerning the task of creating a basic orbital or stellar simulator, providing the ability to designate a known starting point for each of N bodies such as stars: How difficult will it be to get something that reasonably approximates how two to a dozen stars or planets will orbit each other? This is just the mathematics of it all and ignores the problems of display on a computer screen.
LoneDragon
May3-08, 02:24 PM
Try Wikipedia. An astronomical ephemeris (blue books in the library for an american ephemeris I remember for each year) simply lists the positions and velocities of planets, and many other astronomical tables. For a test...
Euler method will *characteristically* work fine for you if all of the bodies don't have close elliptical approaches with large velocities that go though large angles in one step, which integrate increasingly improperly under those conditions in product. Then, if you have properly designed the input mass-position-velocity matrix, and next-step output mass-position-velocity matrix, and Euler step algorithm, and are happy you have worked out the bugs, you **may** replace the Euler test method step calculator, once you see initial Euler results, with a better algorithm. Modularity being key here.
Your code will work fine-enough, but inaccurately, so long as you are not trying to accurately predict, say, precise galactic arm characters, or precise galactic collision simulations.
Also, you can make a stop-gap measure that reduces the step size of close, high velocity, high angle around other body masses, by reducing the step size for *only* those bodies, leaving all of the other bodies at their normal step size. This is what is called an adaptive step size algorithm.
Incidentally, interestingly for the mathematicians only I guess, is that in the differential equations classifications, Euler's method is a first order RK, a second order modified Euler method is a second order RK (RK2), RK4 is obviously a fourth order RK and is closely related to the conceptually easy to understand Simpson's Rule for integration when done in one dimension, and there are higher order RK with more Taylor series expansion terms. But you can start with RK1 (Euler) and work your way up as you pick up the math.
I hope your simulation works well enough for you. I remember doing a solar system simulation on a commodore 64 in junior high which was alot of "fun" (please don't calculate my age) *grins*. Do make sure to calculate all of the bodies at once, from input array to output array in ping pong buffers, for if you advance through just one array and place the answers back in the array for the next step, it causes even more errors in a simple characteristic simulation as half the bodies are advanced, say, while the others are distorted in their interactions. To save time from memory swapping, just toggle referencing the ping pong buffers input-output back and forth with pointers, if your using C. C#, yer on yer own!
regards
Man you guys are smart. I too hope to become as smart as you guys one day. As of Right now I will use the excuse of being young. hehe
LoneDragon
May3-08, 03:28 PM
Man you guys are smart. I too hope to become as smart as you guys one day. As of Right now I will use the excuse of being young. hehe
LOL. Just study hard, read *everything* you can, and take whatever classes are possible *while paying attention*. Time will do the rest. But join the club, as I couldn't even multiply numbers until 3rd grade (I know, that sucks bad) ... hmmm, there must have been something in the water, causing the math spurt over the decades since *then*. LOL, then groans at self.
For bkelly:
To reiterate, as an engineer, I can say Euler will "git 'er done", for a quick peek at simulation programming, sufficient for a *basic* gravitation simulation, but with a strong emphasis on the *basic* simulation. So, "caveat emptor", or, "let the buyer beware", if this is a college math, and coding precision test,, versus a high school, personal interest, or introductory simulation level project. Still, it is good to study the other stuff, if you have the time, when regarding a related technical career.
Man you guys are smart. I too hope to become as smart as you guys one day. As of Right now I will use the excuse of being young. hehe
Hey! Don't include me in this catagory. I'm not BSing my way over someones head to give my gotta-be-smart-ego a plateform. But, hey, I could give it a whirl I 'spose. Let's all meet for tea on friday and we'll try it again.
Redbelly98
May3-08, 05:40 PM
I wrote an 3D, N-body orbital simulator, with displayed output, in basically the same manner you outline. Variables were updated, for the most part, in parallel. The model suffered orbits that would gradually increase in radius over iteration. I repaired it by breaking the parallel-set rule. Near collisions introduced a problem that I never did fix before loosing interest. A close approach would throw the two bodies out towards infinity. I attempted a conservation-of-energy kludge to fix the problem. It was a bust.
Too bad yours didn't work out. I did something similar (in Visual Basic, on a Windows PC) 3 or 4 years ago and was reasonably happy with the results. I remember modifying the force so that it didn't diverge when the bodies got too close to each other, for the reason you saw: all hell breaks loose otherwise. Of course that meant some inaccuracies in the trajectories, but for example I was able to give a cool demonstration of the gravity-assist technique for boosting satellites.
If you're interested I can try to look up exactly how I did the modification to the force.
FYI, I only worked with up to 10 or 20 objects at a time, the computations really bogged down beyond that (Remember, this was Windows and Visual Basic).
LoneDragon
May3-08, 06:04 PM
Hey! Don't include me in this catagory. I'm not BSing my way over someones head to give my gotta-be-smart-ego a plateform. But, hey, I could give it a whirl I 'spose. Let's all meet for tea on friday and we'll try it again.
Well, I am curious for some enlightenment, as you're speaking above *me*, now. Where are the inaccuracies, in the long list of posts, on this thread? Everyone seems to be on an even keel, as far as I can tell, in my humble engineer's math. I see you are an aerodynamics student, or engineer, from your other posts, so, sir, you can take a whirl if you know *that* level of simulation math, as I, at least, would like to know your own numerical opinions, and knowledge, to correct the errant information on this thread. You know, so that we may all be made higher.
Redbelly98
May3-08, 06:45 PM
Wow, that was a lot to say. This will be more complicated that I imagined, even considering that I knew it would be more complicated than I knew. (Does anyone know of a happy face icon with a recursive grin?)
Maybe the best thing right now is to use some basic algorithm, that you can easily understand, to get something working first. Then worry about making it more accurate or complex as a later project.
I found a simulator here: http :// kornelix . squarespace . com / galaxy /
(it seems I don't have URL priviledge, but google on "stellar simulation")
It runs under Linux so I fired up that computer and installed it. It runs well, but not much in the way of features. I haven't reviewed the code closely, but it seems to just use the F=MA type of approach. If I can get version 1.0 running, then I can look into better accuracy for close approaches.
There are two primary segments of this projects
Calculating the positions
Displaying the field
In this forum I want to discuss the calculations.
What type of coordinate system would be best, polar, Cartesian, or something else I don't know about? It seems to me that Cartesian would be better. If for no other reason, I have never messed with polar coordinates. But if polar would be much better, I would like to know.
If Cartesian, I suspect that it would be easier to deal with if the origin is in one of the corners rather than in the middle. My first thought would be the lower left corner. Thoughts about selecting the range: Astronomical Units, millions of miles, something else?
BTW: This site seems to respond much faster under Linux with Firefox than under XP with IE-7.
Too bad yours didn't work out. I did something similar (in Visual Basic, on a Windows PC) 3 or 4 years ago and was reasonably happy with the results. I remember modifying the force so that it didn't diverge when the bodies got too close to each other, for the reason you saw: all hell breaks loose otherwise. Of course that meant some inaccuracies in the trajectories, but for example I was able to give a cool demonstration of the gravity-assist technique for boosting satellites.
If you're interested I can try to look up exactly how I did the modification to the force.
FYI, I only worked with up to 10 or 20 objects at a time, the computations really bogged down beyond that (Remember, this was Windows and Visual Basic).
bkelley's the one you want to talk to about code. I wrote mine a few moons ago before they has invented rocks. Did you solve or notice the walk-away problem where circular orbits are really outward spirals? You can check if it's happening with a running total-energy value.
I hope I didn't PO lonedragon excessively...
Redbelly98
May3-08, 09:26 PM
bkelly,
In Cartesian coordinates it's easier to convert from the (x,y) calculated positions to displayed pixel locations, as only simple scaling factors are involved. Might also be easier to put the origin at pixel (0,0) in the display in the upper left, but if you don't it's still pretty straightforward to rescale from (x,y) to pixel location.
I remember an algorithm a former physics teacher gave me. Better than Euler, not as accurate as Runge-Kutta, but pretty straightforward in terms of understanding it.
In 1-D it works like this ...
Do a Taylor series expansion to find x_{n+1} given x_n , where x_n = x(n \cdot \Delta t)
:
x_{n+1} = x_n + v_n \cdot \Delta t + \frac{1}{2} a_n \cdot \Delta t^2 + \frac{1}{6}x'''_n \cdot \Delta t^3 + ...
Similarly,
x_{n-1} = x_n - v_n \cdot \Delta t + \frac{1}{2} a_n \cdot \Delta t^2 - \frac{1}{6}x'''_n \cdot \Delta t^3 + ...
where v_n, a_n, x'''_n are the velocity, acceleration, and 3rd derivitive of x, respectively, at time t= n \cdot \Delta t
Add these two equations:
x_{n+1} + x_{n-1} = 2 \cdot x_n + a_n \cdot \Delta t^2 + ...
or
x_{n+1} = 2 \cdot x_n - x_{n-1} + a_n \cdot \Delta t^2 + ...
Note the cancellation of the odd-order terms, so that the error is actually of order \Delta t^4, not \Delta t^3. All you need to calculate x_{n+1} are the previous 2 values of x, as well as the acceleration at time n \cdot \Delta t, which you calculate from the Forces and masses in the system.
In 2-D you'd need to calculate the x- and y-components of the accelerations, then apply the algorithm to find x_{n+1} and y_{n+1} .
There might be a common name for this algorithm, if anybody knows it feel free to post the info.
Redbelly98
May3-08, 09:55 PM
bkelley's the one you want to talk to about code. I wrote mine a few moons ago before they has invented rocks. Did you solve or notice the walk-away problem where circular orbits are really outward spirals? You can check if it's happening with a running total-energy value.
I didn't check circular orbits, that would have been a good test. What I did notice was that if two objects got very close, they could suddenly go flying off the screen at incredibly high speed.
p.s. I just looked over my code, and found that I used the algorithm I described in message #21. Somehow I thought I later changed this to Runge-Kutta but apparently not.
I can't figure out, from a quick glance at the code, just what I'm doing to the force to make it "behave" at close range. I would have to find some 3-year-old handwritten notes that I have somewhere to figure that out again.
Anyone feel like sharing one of their simulators? Would be interesting.
The guys at NASA, NSA, and whoever, should have devoted quit a few man-years at this. So they may have found a different method than interating t, and updating V_i and X_i. There's the Hamiltonian, if you've heard of it. It's an alternative, equivalent to Newtonian mechanics. What I'm trying to say, is that it may be possible to iterate values in some other mathematical space, rather than velocity and position. Angular momentum and Kinetic energy perhaps; I have no idea really. Just offinging possibilities.
I wrote an 3D, N-body orbital simulator, with displayed output, in basically the same manner you outline. Variables were updated, for the most part, in parallel. The model suffered orbits that would gradually increase in radius over iteration.This is a direct consequence of using the Euler method. The Euler method displays this problem even with two bodies and even without any eccentricity.
I remember an algorithm a former physics teacher gave me. Better than Euler, not as accurate as Runge-Kutta, but pretty straightforward in terms of understanding it. ... There might be a common name for this algorithm, if anybody knows it feel free to post the info. The method you described is the basic verlet method. As you noted, it is better than Euler. An even better approach and still very simple to understand is Heun's method.
I found a simulator here: http:// kornelix.squarespace.com/galaxy/. It runs well, but not much in the way of features. I haven't reviewed the code closely, but it seems to just use the F=MA type of approach. The cited simulator updates the velocity and then updates the position using the updated velocity. This is not a very good integrator. A simple improvement is to update the position using the average of the velocity before and velocity after the velocity update. This simple improvement is Heun's method, a second-order Runge-Kutta integrator. Heun's method is considerably better than the midpoint method, another second-order Runge-Kutta integrator.
I haven't reviewed the code closely, but it seems to just use the F=MA type of approach.*All* of the Runge-Kutta techniques as well as the verlet method "just use the F=Ma" type of approach. Some just do it a lot better than others.
I didn't check circular orbits, that would have been a good test.That is always a good test. A simulator that can't do a very good job reproducing a circular orbit isn't worth very much at all. Highly elliptical orbits and hyperbolic orbits are also good tests because the close approaches will stress even good integrators.
The guys at NASA, NSA, and whoever, should have devoted quit a few man-years at this. So they may have found a different method than interating t, and updating V_i and X_i. There's the Hamiltonian, if you've heard of it. It's an alternative, equivalent to Newtonian mechanics.
NSA, as far as I know, doesn't address this issue. OTOH, NASA certainly does. Hamiltonian-based techniques are called symplectic integrators; google this phrase for more info. The advantage is that these techniques inherently do a better job at conserving energy than do non-symplectic integrators (e.g., the Euler method and the standard RK4 method). The basic verlet method is a symplectic integrator, for example, and it does a much better job than does the Euler method even though both evaluate state derivatives once for each propagation interval. I believe JPL uses symplectic integration techniques for developing ephemerides.
There are many, many other orbital propagation techniques. One approach is to use Keplerian elements (e.g., a,e,i,\Omega,\omega,M) rather than position and velocity. In the case of two bodies, all Keplerian elements but the mean anomaly are constant and the mean anomaly is a linear function of time. Propagation is easy and exact. Making sense out of the elements (e.g., plot positions on a computer screen) takes a little work, but not too much. Now consider a stellar system which comprises a single star and several planets. The planets' orbits can still be described in terms of Keplerian elements, but these now exhibit time-varying behavior. The variations are small, so perturbation techniques can be put into play. NASA used such an approach during the Apollo era, and still uses it in some arenas today.
... There are many, many other orbital propagation techniques. One approach is to use Keplerian elements (e.g., a,e,i,\Omega,\omega,M) rather than position and velocity. In the case of two bodies, all Keplerian elements but the mean anomaly are constant and the mean anomaly is a linear function of time. Propagation is easy and exact. Making sense out of the elements (e.g., plot positions on a computer screen) takes a little work, but not too much. Now consider a stellar system which comprises a single star and several planets. The planets' orbits can still be described in terms of Keplerian elements, but these now exhibit time-varying behavior. The variations are small, so perturbation techniques can be put into play. NASA used such an approach during the Apollo era, and still uses it in some arenas today.
I hadn't considered perturbations from elliptical orbits, but then again, I haven't considered orbital simulations in a long while.
An interesting possibility to avoid close-approach divergence for any pair of bodies, might be to switch to elliptical orbits + perturbations (for the pair only) based upon some nearness condition, then switch back. It may even be possible to survive a direct hit where the mass centers pass through one another.
In the above I failed to mention the possibility of close encounders of hyperbolic and parabolic trajectories where these might also be degenerate (head on collisions).
lzkelley
May10-08, 11:04 PM
I'm also trying to make an n-body orbital simulator, its going quite well - i've recently upgraded from euler to fourth order runga-kutta.
I've been using Cartesian coordinates for all of my work so far, and for the actual calculations i plan to continue - but i need to be able to use keplerian orbital elements instead of state vectors for initial conditions. I was wondering if anyone out there had experience or direction on how to convert between orbital elements and state vectors besides drowning myself in geometry. Thanks!
csprof2000
May15-08, 11:41 AM
Regardless of what everyone is telling you, your method will work fine. Just be prepared to use very small time steps if you are using Euler's method.
Besides inaccuracies that can occur due to the approximate nature of the solution, the near-collision "explosion" can be solved by simply defining the force piecewise. Calculate the distance between each pair of bodies (this must be done anyway) and see if it's less than some minimum value (this can be as simple or as complicated as you like). Then decide whether or not you want to "do" the force there, and act accordingly. Otherwise, you could use a similar method to use conservation of momentum to have them collide (elastically or inelastically) like billiard balls. Or have them lump together and run until all masses have lumped together.
There's a lot of potential for fun and good programming in what you're talking about.Good graphics can be done in Java, C, C++, C#, Visual Basic, etc. For 3-D, you may even look into DirectX or OpenGL.
The programming won't be hard, but there are some things you can do to make it easy on yourself. Look for libraries with functions / objects you can use. I recommend taking the opportunity to explore OO approaches, if you're not already familiar with these. You can make objects for vectors and massy particles and this will greatly ease the notation. Have fun with it.
Ulysees
May17-08, 08:50 AM
May I ask, for a given accuracy of 10^-5 say, and for a given orbit lasting 1 year, which integration method is the fastest?
csprof2000
May17-08, 09:43 AM
I would suspect that you can't do much better than the fourth order adaptive stepsize Runge Kutta algorithm.
Ulysees
May17-08, 11:33 AM
Does that mean repeating with a smaller time interval until all 3 steps of 4th order Runge-Kutta agree within an error margin?
csprof2000
May17-08, 12:51 PM
I'm not sure precisely what happens, but I'm pretty sure that a timestep is taken with one RK method and then another is taken with one which is of one higher order, and if the difference is within the error bars, it goes with it. Otherwise, it tells you what the timestep should be to get your accuracy. It has the benefit of remaining as accurate as you tell it (within reason) and not being *too* when the forces are negligibly small (for instance, why use a timestep of h = 1.0e-10 both when the force is large and changing quickly AND when the force is next to nothing and staying put?)
tony873004
May23-08, 01:20 PM
You should start with Euler's method because it is simple and you understand it. Then if you like, you can upgrade your method later.
I wrote Gravity Simulator: www.gravitysimulator.com , which sounds like what you want to do. It uses Euler's method. Using Euler's method I can:
*Predict moon phases at least 50 years into the future
*Propogate Apophis 22 years into the future, and miss by only 200 km its close approach distance to Earth in 2029 (which btw is also the amount the higher-order Picard's method misses by, so the inaccuracy seems to have more to do with forces not modeled, rather than with integrator inaccuracies).
*Correctly predict the Newtonian predicted value of the precession of Mercury's perihelion
*accurately recreate virtually any point-mass Newtonian simulation done by professional astronomers except for the ones that last billions of years.
*and the list goes on. My web site has more examples
I'm not an expert on the higher order routines, but most of the bad things I've heard about Euler's method seem to be overrated and seem to disappear with a slow enough time step.
In the real universe, there are lots of other forces besides Newtonian gravity from point-mass sources: radiation pressure, relativity, non-spherical objects, etc. I'm guessing that errors caused by ignoring these overwhelm the errors caused by Euler's method. For example, the toughest thing to compute is the orbits of tight moons such as Phobos and Deimos. After 2 years of simulating, even at a slow time step, they're not where they are supposed to be. The size, shape and orientation of their orbits are still fine, but their position on that orbit is off, by lots for Phobos, and by a little for the more distant Deimos. But this is more likely due to treating Mars as spherical, than with numerical errors. So even a higher order routine will miss this one too.
I recently wrote a 2nd order routine. It was more accurate than Eulers at time steps less than 64 seconds, but less accurate than Euler at time steps higher than 64 seconds. In the 2nd order routine, the planets would spiral out. It was explained to me like this: In Euler (1st order) every time you cut your time step in half, you cut your error in half. In 2nd order routines, every time you cut your time step in half, you cut your error 4x. But the inverse of this is true too. Every time you double your time step, your error grows by 4x for the 2nd order routine and by only 2x for Euler. Since most simulations I like to perform require a time step much high than 64 seconds, I won't find much use for the 2nd order routine. Next I'd like to try RK4. It will be interesting to compare to the current Euler routine. Will I like it better, given that it seems like Euler is very underrated?
Don't be discouraged from using Euler if you don't understand the higher order methods. You'll have lots of fun with it. And don't be scared away from 3d. It's not that much different from 2d. Pothaorean theorem works in 3d too: d=sqrt(x^2+y^2+z^2), a=GM/d^2, r=r+a
Redbelly98
May23-08, 07:48 PM
You should start with Euler's method because it is simple and you understand it. Then if you like, you can upgrade your method later.
Definitely good advice here. Get something working so you can start having fun with it. Then worry about what is the "best" method.
tony873004
May24-08, 12:27 PM
...I've been using Cartesian coordinates for all of my work so far, and for the actual calculations i plan to continue - but i need to be able to use keplerian orbital elements instead of state vectors for initial conditions. I was wondering if anyone out there had experience or direction on how to convert between orbital elements and state vectors besides drowning myself in geometry. Thanks!
Why do you need to use Keplerian orbital elements as your initial conditions? From JPL Horizons, you can get state vectors.
I've got code for converting both directions from state vectors to orbital elements and back again. Send me a PM if you want a copy.
tony873004
May24-08, 01:35 PM
One thing that needs to be asked is why does accuracy matter? Are you trying to predict solar eclipses far into the future? Are you trying to pilot a spacecraft to Mars where missing by just a few kilometers will doom your mission? Or are you simply trying to demonstrate a gravity assist, Kozai mechanism, precession of nodes, etc.? For most of these accuracy isn't really that important. I'll show a few examples of tests I've performed:
This is an animation of the Moon's orbit around the Earth. It was made with Gravity Simulator's 1st order Euler's method. The time step was 16000 seconds. At this time step, the simulation is very innacurate. This would be evident if you look at the Moon's position at a known future date. For example, on August 12, 2045 there will be a total eclipse of the Sun across America. So if I paused the simulation there, I should see the Moon directly between the Earth and Sun. But with the errors introduced by the numerical method, the Moon will likely be elsewhere in its orbit at this time. However, despite the error, many features of the Moon's orbit remain accurate: Its semi-major axis has held constant. Its eccentricity and inclination remained constant as well. And its nodes still precess at their 18 year cycle known as the saros cycle, as evident in this animation. In fact, the saros cycle remains consistent no matter how large of a time step I take as long as it is slow enough that the Moon remains in its orbit.
http://www.orbitsimulator.com/gravity/saros.GIF
This next image compares the the precession of Mercury's perihelion between the Newtonian-predicted value of 531 arcseconds per century, Euler's method, and a 2nd-order method:
http://www.orbitsimulator.com/astr699/mercuryprecession.GIF
Euler's method actually out-performs the 2nd order method. It requires time steps of less than 64 seconds before the 2nd order method begins to catch up. But notice that every time the time step is halved, the 2nd order method's error is cut by 4x, and Euler's method is only cut 2x. So there will reach a point where 2nd order is more accurate, but it will be when the time step is very slow.
These next two images show a massless object on a chaotic trajectory in the Earth/Moon system. The simulation is run for a few months, then the color of the object is changed from magenta to green, and time is reversed. If everything is perfect, it will retrace its same path. The time step used is 16 seconds.
http://www.orbitsimulator.com/astr699/euler.GIF
http://www.orbitsimulator.com/astr699/secondorder.GIF
You'll notice that the 2nd order method is nearly perfect, while Euler is not. Is this a victory for the 2nd order method? I'm not sure. Some routines have the feature that they're perfectly time-reversable. So even if there are errors, the errors will be mirrored exactly when time is reversed. I'm not sure if my 2nd order routine has this feature. I know that Euler does not.
Just for fun, here's my failed attempt at the Leapfrog method. I obviously did something wrong
http://www.orbitsimulator.com/astr699/leapfrog.GIF
Here's a graph created with Gravity Simulator (Euler) showing Earth's changing eccentricity over the course of 1 million years. It clearly shows the small 100,000 year cycles and the larger 400,000 year cycle predicted by analytical methods. Earth's orbit is chaotic enough that I would imagine that the graph would be quickly out of phase regardless of the method used, although the 100,000 year and 400,000 cycles would remain evident despite the phase error.
http://www.orbitsimulator.com/astr699/EarthEcc.GIF
In the real universe, there are lots of other forces besides Newtonian gravity from point-mass sources: radiation pressure, relativity, non-spherical objects, etc. I'm guessing that errors caused by ignoring these overwhelm the errors caused by Euler's method.
No, they don't. Moreover, if your orbit simulator maintains the Moon's semimajor axis, then you are not using the basic Euler method. You are most likely using the Euler-Cromer method instead. The two methods:
Euler method
\begin{aligned}
\mathbf x(t+\Delta t) &= \mathbf x(t) + \mathbf v(t)\Delta t \\
\mathbf v(t+\Delta t) &= \mathbf v(t) + \mathbf a(\mathbf x(t))\Delta t
\end{aligned}
Euler-Cromer method
\begin{aligned}
\mathbf v(t+\Delta t) &= \mathbf v(t) + \mathbf a(\mathbf x(t))\Delta t \\
\mathbf x(t+\Delta t) &= \mathbf x(t) + \mathbf v(t+\Delta t)\Delta t
\end{aligned}
The two methods differ only in the calculation of the updated position vector. In the basic Euler method, the pre-update velocity is used while in Euler-Cromer, the position is updated using the post-update velocity. This slight change makes a big difference. Euler-Cromer is a symplectic integrator, meaning that it more-or-less conserves energy in situations with a time-independent Hamiltonian (e.g., gravity). The degree to which it fails to conserve energy is a function of step size and numerical errors inherent in using finite-precision arithmetic. The basic Euler method is not symplectic; orbits spiral out.
tony873004
May24-08, 02:47 PM
Thanks for that, DH. My method updates velocity and then position. I never knew that was called Euler-Cromer. I actually made it up myself. (can we call it Euler-Cromer-Tony? :) ).
I had always heard that Euler should cause the orbits to spiral out, so I wondered why mine never did.
Ulysees
May24-08, 03:28 PM
I should add the more important problem with Euler variants is not accuracy. It is stability. For example a mass attached to a spring without any damping should oscillate forever (that's the analytical solution), but with Euler it oscillates with an increasing amplitude. No matter how small you make the integration step, the amplitude of oscillation is increasing.
4th order Runge Kutta is not so unstable. So its 3 steps are better than 3 normal steps of Euler.
I should add the more important problem with Euler variants is not accuracy. It is stability. For example a mass attached to a spring without any damping should oscillate forever (that's the analytical solution), but with Euler it oscillates with an increasing amplitude. No matter how small you make the integration step, the amplitude of oscillation is increasing.
That is always a problem with the basic Euler method. It is not so much a problem with the Euler-Cromer method, which like RK4, exhibits conditionally stable (but the stability region is considerably smaller with Euler-Cromer than RK4). The simple act of updating the position with the already-updated velocity makes Euler-Cromer conditionally stable and symplectic. RK4 is considerably more stable and considerably more accurate than Euler-Cromer, but it is not symplectic.
Ulysees
May24-08, 08:01 PM
And what's the benefit of symplectic-ness?
And what's the benefit of symplectic-ness?
Ignoring errors resulting from the use of finite-precision arithmetic, a symplectic integrator comes very close to conserving energy. How close depends on the order of the integrator and the product of the integrator step size and the largest critical frequency of the system. When applied to a two body system, orbits will remain at more-or-less the size and shape with a symplectic integrator.
This is not the case for non-symplectic integrators such as basic Euler or the RK methods. With Euler the deviation is immediate. Orbits spiral out quickly. With RK4 the degradation is much, much slower, but still occurs. An RK4-based integrator will be closer to truth than the simplest symplectic integrator, Euler-Cromer, when used over just a few orbits. On the other hand, using an RK4 integrator to propagate over hundreds of orbits is a bad idea. At some point the energy will have changed enough to make all results garbage. Euler-Cromer suffers from reduced accuracy but benefits from increased global stability.
Just as RK4 is much better than RK1 (i.e., Euler's method), higher-order symplectic integrators are a lot better than Euler-Cromer. Higher-order methods include Verlet methods, Forest-Ruth, and Chin. Chin's Algorithm C is the new gold standard.
tony873004
May24-08, 10:13 PM
Your replies to my posts and Ulysees posts are very facinating. Do you have any links describing Forest-Ruth, Verlet or Chin?
Something I wanted to work on this summer was a higher order integrator for Gravity Simulator. I thought I wanted to do RK4, but you may have discouraged me from that.
Why would anybody want to use Euler's method for anything, being that if you just solve for v before you solve for x, it becomes much nicer? Didn't Euler realize this?
Your replies to my posts and Ulysees posts are very facinating. Do you have any links describing Forest-Ruth, Verlet or Chin?
Something I wanted to work on this summer was a higher order integrator for Gravity Simulator. I thought I wanted to do RK4, but you may have discouraged me from that.
Why would anybody want to use Euler's method for anything, being that if you just solve for v before you solve for x, it becomes much nicer? Didn't Euler realize this?
RK4 is a very, very good single-step integrator. (All of these techniques are single-step integrators which means that they don't keep a history. The three internal steps in RK4 are all part of the same step.) RK4 is typically more accurate than all the techniques mentioned but Chin's Algorithm C. That RK4 is not symplectic only becomes a significant factor when the integration is performed over a very long time span.
Ideally your simulator should be constructed so that the propagation technique is completely independent of the rest of the simulation. You should be able to switch propagators without changing any code. This is a completely achievable task and how easily it can be done is a good test of the quality of your code. (Chin's Algorithm C is an exception as it requires the gradient of the acceleration for one internal step.)
Euler died 225 years ago. Give him a break. He predates the concept of a Hamiltonian by 100 years and predates the concept of a vector by even more than that. The simple act of using the post-update velocity to update the position looks like a dumb mistake from a mathematical perspective -- unless one looks at things from the perspective of a Hamiltonian.
As far as references,
Verlet. Use Google. The first page of hits are all true hits, including the two Wiki articles.
RK4. The wikipedia article on RK4 (http://en.wikipedia.org/wiki/Runge-Kutta#The_classical_fourth-order_Runge.E2.80.93Kutta_method) is compact and accurate. Other articles on Runge-Kutta include http://www.nrbook.com/a/bookcpdf/c16-1.pdf and http://numericalmethods.eng.usf.edu/ebooks/runge4th_08ode_ebook.pdf.
Chin and Forest-Ruth. You can find Chin's papers at arXiv.org (http://arxiv.org/find/grp_physics/1/au:+chin_siu/0/1/0/all/0/1). His papers also discuss FR since his method is an offshoot of that method.
Ulysees
May25-08, 08:07 AM
Have you done any long-term simulations like the formation of a solar system? I was looking for a simulator suitable for testing Titius-Bode law.
Also, there are some who say that if a planet is hit by a large enough object or otherwise explodes, the entire solar system will collapse eventually. Any thoughts on this?
tony873004
May27-08, 05:19 PM
Euler died 225 years ago. Give him a break.
Fair enough. I'll cut his some slack! Thanks for the links.
Have you done any long-term simulations like the formation of a solar system? I was looking for a simulator suitable for testing Titius-Bode law.
The formation of the solar system is too complicated to simulate with point-mass n-body code. When the primordial particles are small, gravity does not play a roll. Once they're large enough for gravity to play a role, certain collisions cause them to stick together, while others break them apart. And there's millions of them. My simulator does not do very well with even hundreds of particles.
However, if you're willing to make some assumptions, there's a few simulations you can run. Check this one out. It's called "moonbuilder": http://www.orbitsimulator.com/gravity/articles/moonbuilder.html . It begins with 100 particles each with 1/50 the mass of the real moon. All collisions cause them to merge together. It only takes a few weeks for the hundred to merge together into 2 moons, one that remains in orbit, and one that gets ejected. And different starting conditions will yield different results, so if you try it again and again, you may get systems with 2 or 3 moons. The more you try it, the more you should see some trends develop.
I've tried the same thing, except placing 100 particles at 1 AU from the Sun to form a planet. Sometimes I get trojans or horseshoe objects. Again, every run is different, but long-term trends can be spotted. So if you're clever, and willing to make some simplifications and assumptions, there are definately some things you can try.
I don't think it would do to well trying to confirm or reject Titus-Bode. You'd be trying to build multiple planets at the same time, which would require lots of particles. You could play around with "Moonbuilder", perhaps spacing the initial particles further apart, and trying to find the ideal conditions that lead to 2 or 3 stable moons, and see if there's any kind of trend to the spacing. But I'd be careful about the conclusions you draw with such a simplified model.
Also, there are some who say that if a planet is hit by a large enough object or otherwise explodes, the entire solar system will collapse eventually. Any thoughts on this?
I've heard things like Earth's presence supresses a long-term resonance between Venus and Jupiter, and without Earth, Venus would not be stable. Maybe that's what you're referring to? I'm not sure how accurate that statement is. When I heard it, the first thing I tried to do was simulate the solar system without Earth. Venus stayed put, although I only simulated for a few million years. It's possible that it takes longer than this for Venus to get ejected, and its also possible that my time step was high enough that numerical errors dominated, and I wont get an ejection no matter how long I simulate. It takes a pretty big time step to get millions of years into the future.
Ulysees
May27-08, 06:33 PM
Ιsn't it reasonable that Titius-Bode would be close to being satisfied with just 1000 particles? Why do you believe the statistics would only converge to an analytical pattern such as Titius Bode only with zillions of particles? We don't need to poll millions of voters in order to predict the outcome of an election, just a 1000 voters converge to a pattern too.
tony873004
May28-08, 12:33 AM
That was the point I was trying to make: set up what you feel would be a representative system using fewer particles. But you need to understand how your model differs from reality. The newly-formed planets interact with the intermediate belts of trillions of lefover particles. This causes their orbits to circularize and to migrate to their final positions. That would have to be modeled too.
Ulysees
May28-08, 06:06 AM
That would have to be modeled too.
And what's stopping you from imagining a handful of planets made up from 100 particles each, rotating in a cloud of 500 other particles, and this occuring automatically starting from an initial cloud. I only want to test Titius Bode as a mathematical property that should originate from the inverse square law in any scale, atomic or galactic. Not confine it to the actual masses, distances, etc of a solar system.
tony873004
May28-08, 12:43 PM
And what's stopping you...
A lack of interest in Titus-Bode prevents me from trying to verify it through simulations. But if it interests you, that doesn't need to stop you. Anyone with a Windows computer can run my simulator.
...I only want to test Titius Bode as a mathematical property that should originate from the inverse square law in any scale, atomic or galactic...
There are 6 opportunities in our solar system for Titus-Bode to be revealed: the planetary system, and 5 moon systems (>2 moons: Jupiter, Saturn, Uranus, Neptune, Pluto), none of which demonstrate Titus-Bode. The planetary system fails because there's no major planet between Mars and Jupiter, and Neptune is out of place.
Modeling the formation of planetary system is a vastly more complex task than a simple orbital simulator. Planetary growth is a complex and not yet fully understood process. I assume Tony's simulator treats the planets as point masses. This is a good assumption for a stabilized system with only a small number of orbiting objects. It is not a good assumption for modeling planet growth. Tony's simulator runs on a single off-the-shelf computer. Modeling protoplanetary systems requires a lot of coupled, high powered computers because one needs to model millions of tiny particles. These little particles collide (point masses don't collide), have angular momentum due to rotation (point masses don't have angular momentum due to rotation) and exchange angular momentum with the Sun and each other. Once again, point masses don't do this.
tony873004
May28-08, 02:32 PM
...I assume Tony's simulator treats the planets as point masses....These little particles collide (point masses don't collide),..
They're treated as point masses for the sake of computing acceleration, but they do have size, so they can collide and stick. But every collision simply combines their mass and momentum into a new object if their distance is less than their combined radii. But the points you make are what I was getting at a few posts ago (#48). The real system has lots of stuff going on besides point mass interactions. There are more types of collisions than the simple one my code models. And the real system has gas that behaves differently than solid matter. So I can only agree with you that this type of a simulator is too simple for task. It doesn't mean than you can't make some assumptions and generate simplified models to explore various aspects. But to build a solar system from scratch, which you can draw confident conclusions about, is well beyond its capabilities.
Ulysees
May29-08, 06:20 PM
I assume Tony's simulator treats the planets as point masses. This is a good assumption for a stabilized system with only a small number of orbiting objects. It is not a good assumption for modeling planet growth.
Well don't forget Newton's law of gravity only makes sense for point masses. Everything else is approximated by point masses, for the purposes of this law. Conglomerates of point masses can have angular momentum but there is no need to introduce angular momentum in the model. Just point masses connected with breakable and rejoinable spring-damper links, should be enough to model just about anything we might be interested in gravity-wise, and I think Tony is not far from this. The difficulty is with the estimation of suitable parameters for something that matches the real thing, but then we are not always interested in matching the real thing but just want to explore possibilities.
Regarding Titius-Bode that I am interested in testing, there is no reason to assume it applies exclusively to the scale of our solar system. I would expect it to also apply to much bigger and much smaller systems, because it seems to be a statistical result, like Gaussian distributions that appear in nature.
By the way, Titius-Bode -like equations predict a planet between Mars and Jupiter, and indeed there is the asteroid belt there, ie a planet that failed to form or a planet that formed and was later broken apart somehow, perhaps in a collision with a large object like the one that hit the earth when the moon was formed.
Modeling protoplanetary systems requires a lot of coupled, high powered computers because one needs to model millions of tiny particles.
There's a simulation of the collision of the Mars sized object with the earth when the moon was formed. Judging from the number of particles, it seems able to run on todays pc's. It ends up with the moon at the right distance from earth, and the right masses for the earth and moon.
Well don't forget Newton's law of gravity only makes sense for point masses.
Wrong. It is an trivial matter to compute the gravitational influence of a non-point mass body that has spherical shape and a uniform density on some object. An ellipsoidal body is not too much harder. We use complex mathematical models such spherical harmonics to describe real body such as the Earth. Some spherical harmonics models of the Earth's gravity field include EGM96 (http://cddis.nasa.gov/926/egm96/egm96.html) and the GRACE GGMO2 gravity models (http://www.csr.utexas.edu/grace/gravity/). For the Moon, there are the GLGM-2 (http://cat.inist.fr/?aModele=afficheN&cpsidt=2847930) and more recent ones based on Lunar Prospector data. JPL has developed spherical harmonics models for Venus, Mars, and even some asteroids.
Regarding Titius-Bode that I am interested in testing, there is no reason to assume it applies exclusively to the scale of our solar system.
One of the cornerstones of the scientific method is the concept of falsification. Neptune falsifies the Titus-Bode law. Most astronomers view the Titus-Bode law as mere numerology. I suggest you post your concepts on this law in the scepticism and debunking forum at PF.
Ulysees
May29-08, 07:37 PM
Newtons law of gravity does not apply to ellipsoid shapes last time I checked, only spherical shells are equivalent to point masses, and only on the outside (shell theorem).
But if you have proof that ellipsoid shapes are equivalent to point masses too, as required by Newton's law, I'd like to see it.
Or if you have proof that ALL spherical harmonics are equivalent to point masses for Newton's law to apply, I'd like to see it.
Cause what I have seen from my numerical calculations (in thread
Numerical integration of gravity at surface of an object (http://www.physicsforums.com/showthread.php?t=212711)), Newton's law of gravity only matches spherical shells. And therefore it can only be used exactly for spherically symmetrical distributions of matter like an onion, everything else is by appoximation.
I never said ellipsoidal shapes are equivalent to point masses. I said one can compute the gravitational attraction toward an ellipsoidal body. The gravitational attraction at some point \mathbf x outside some body is volume integral
\mathbf a(\mathbf x) = G\int_V \frac{\rho(\mathbf {\xi})}{||\mathbf {\xi}-\mathbf x||^3}(\mathbf {\xi}-\mathbf x)\,d{\mathbf {\xi}
where the integration is taken over the massive body in question.
tony873004
May29-08, 07:59 PM
Newton's Laws: (from memory)
1. An object at rest remains at rest, or an object in motion stays in linear constant motion unless acted upon by a force.
2. F=ma
3. For every action, there is an equal but opposite reaction.
There's nothing here that excludes ellipsoids. Newton spent a lot of time developing Calculus simply so he could mathamatically justify his assumption that a sphere simplifies into a point mass. He did this by integrating the volume of a sphere. Now we can save a lot of time by using this nice shortcut. But there's no reason you can't integrate an ellipsoid, or any other shape, without going beyond Newton's laws. And if you wanted to do it numerically, if I created an ellipsoid from a collection of 1 million point masses, and then simulated a particle orbiting the ellipsoid, using Newton's laws, I imagine it would work quite nicely, showing the proper orbital precession that gives us "walking" orbits and "sun synchronous" orbits.
Ulysees
May29-08, 07:59 PM
I never said ellipsoidal shapes are equivalent to point masses.
You said "wrong" below this:
Well don't forget Newton's law of gravity only makes sense for point masses.
Wrong.
Ulysees
May29-08, 08:03 PM
Newton's Laws: (from memory)
1. An object at rest remains at rest, or an object in motion stays in linear constant motion unless acted upon by a force.
2. F=ma
3. For every action, there is an equal but opposite reaction.
We're talking about Newton's law of gravity, not Newton's laws of motion.
But there's no reason you can't integrate an ellipsoid, or any other shape, without going beyond Newton's laws.
I've done this for a living. It's another thing we're talking about here (whether Newton's law can be used exactly in anything but point masses and equivalent spherical shells)
You said "wrong" below this:Well don't forget Newton's law of gravity only makes sense for point masses.Wrong.
That is because what you wrote is wrong.
I've done this for a living. It's another thing we're talking about here (whether Newton's law can be used exactly in anything but point masses and equivalent spherical shells)
I gave the exact answer in post #58. Whether that integral is computable in closed form or has to be approximated numerically is a different question. The integral is the exact answer to the acceleration due to gravity. BTW, this is exactly what I do for a living (one of the things I do for a living). Moreover, NASA has many people at JPL and Goddard working on developing high-precision gravity models. Until recently they used ephemeris data from satellites that had some other primary objective as the basis for these gravity models. NASA and ESA together are now flying an expensive pair of satellites whose sole objective is developing even higher-precision gravity models of the Earth. Do you think they would this just for fun?
tony873004
May29-08, 08:32 PM
If you're referring to F=GMm/d^2, then yes, this works only for point masses. But throw a double integral symbol in front of it, with the limits of integration describing your desired shape, and you can do any shape you want. You're still applying the point mass formula GMm / d^2 (or one M for acceleration) to each differential element. And if you want something that is not uniform in density, and you have your density function, then write it the way DH did.
Ulysees
May29-08, 08:33 PM
NASA and ESA together are now flying an expensive pair of satellites whose sole objective is developing even higher-precision gravity models of the Earth. Do you think they would this just for fun?
Certainly not, but we were talking about the formation of a solar system from dust, and you are introducing ellipsoids. Don't you think they're from another application?
Ulysees
May29-08, 08:36 PM
throw a double integral symbol in front of it, with the limits of integration describing your desired shape, and you can do any shape you want.
Not in the context of dust forming proto-planets in order to test Titius bode.
Even if an analytical formula came out from integration, you still can't use it in this context.
tony873004
May29-08, 08:43 PM
Not in the context of dust forming proto-planets in order to test Titius bode.
I agree with you. You're not going to be using Newton's laws of gravity for dust. Objects need to be a few km wide before their gravity is significant enough to consider. I don't think it's fully understood yet, what happens prior to this point.
And that's one of many reasons why my program will have a difficult time helping you verify TB. My program is GM/d^2 inserted into an endless loop.
Ulysees
May29-08, 08:52 PM
Why don't you post a copy of your program in case anyone wants to play with it. Or have you already?
Ulysees
May29-08, 08:56 PM
Objects need to be a few km wide before their gravity is significant enough to consider.
Well you'd start with thick "dust" then, wouldn't you. Like that simulation of the formation of the moon from a collision of a large object with Earth: both objects turn into thick "dust" at places.
Ulysses, please stop with the Titus-Bode law. It is numerology. No causal mechanism, and debunked by Neptune.
Ulysees
May29-08, 10:21 PM
Red cars cost more to insure in the UK. Insurance companies do not seek a cause why red cars have higher statistics for accidents. They just use the result and charge more.
So if you do not know a cause for an observed statistical trend, that does not mean no cause exists. It just means you don't know a cause yet. In trials of new drugs they give a license if the drug gives results that are "statistically significant", and they have equations for deciding if a result is
a. statistically significant (and therefore caused by the drug) or
b. not significant (and therefore just chance), and the drug is rejected.
So refusing to even simulate the formation of planets roughly in search of ANY statistical pattern (if it exists it should appear after many repeats), refusing to even test is not a scientific approach, it is dogmatic thinking that has no place in science.
We have now observed other solar systems. In several of these systems, the large gas giants orbit very close to their parent star. The Titus-Bode Law does not explain this phenomenon. In our solar system, the Titus-Bode law predicts a planet between Mars and Jupiter. No planet exists there; just a bunch of asteroids that collectively mass far less than any planet and that could not have formed into a planet because of the disruptive influence of Jupiter. The Titus-Bode Law does not explain this phenomenon. Neptune is also a member of our solar system, and it completely defies the Titus-Bode law.
Falsifiability and falsification are cornerstones of the scientific method. Scientific theories can never be proven true, but they can be proven false by one well-confirmed piece of falsifying evidence. Several scientific hypotheses to explain some phenomena have fallen by the wayside because of evidence against them. The Titus-Bode law is one such rejected hypothesis. It is now viewed as little more than numerology. The analogy to the insurance industry is a false analogy because science operates with the scientific method while the insurance industry does not. The medical industry operates on the edge of the scientific method because of ethics, the complicated nature of drug interactions, and a host of other reasons. A single bad reaction or failed cure is not an adequate reason to reject a drug that does work in some people.
The existence of Neptune is more than enough to reject the Titus-Bode law in physics and astronomy. Please desist with this stuff on the Titus-Bode law.
Aethelwulffe
Jun2-08, 01:11 AM
Back to the subject of creating an orbital simulator, are any of you guys familiar with Orbiter? It is a free space-flight simulator. I myself am a pretty prolific add-on developer for it, having done some mars mission stuff as well as some (well about 700 mb actually) SCI-FI add-ons for it. There is quite a great support community, and a lot of fun things to do with existing packages and add-ons, not to mention the SDK.
This software is being used in lots of space flight simulations for education, as well as a rendering engine for video presentations of space missions. Best of all, it is FREE, has a HUGE in-depth API and SDK, as well as tons of support from lots and lots of top-notch developers that all contribute for free and make their works free to the public, as well as help set up simulations in classrooms all over the world. If nothing else, the global/local coordinate system and other elements would be a good thing to check out for someone starting their own project...though you would be a few million man-years behind that particular project which has been developed since 1999 or so.
Just google Orbiter Space Flight Simulator and check it out if you care to. It's no joke, and it's no cheezy bit of work.
vBulletin® v3.8.7, Copyright ©2000-2012, vBulletin Solutions, Inc.