
#19
May308, 07:48 PM

P: 84

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



#20
May308, 09:09 PM

P: 4,513

I hope I didn't PO lonedragon excessively... 



#21
May308, 09:26 PM

Mentor
P: 11,989

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 RungeKutta, but pretty straightforward in terms of understanding it. In 1D it works like this ... Do a Taylor series expansion to find [tex] x_{n+1} [/tex] given [tex] x_n [/tex] , where [tex]x_n = x(n \cdot \Delta t) [/tex] : [tex] 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 +[/tex] ... Similarly, [tex] x_{n1} = 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 + [/tex] ... where [tex] v_n, a_n, x'''_n [/tex] are the velocity, acceleration, and 3rd derivitive of x, respectively, at time [tex] t= n \cdot \Delta t [/tex] Add these two equations: [tex] x_{n+1} + x_{n1} = 2 \cdot x_n + a_n \cdot \Delta t^2 +[/tex] ... or [tex] x_{n+1} = 2 \cdot x_n  x_{n1} + a_n \cdot \Delta t^2 + [/tex] ... Note the cancellation of the oddorder terms, so that the error is actually of order [tex]\Delta t^4[/tex], not [tex]\Delta t^3[/tex]. All you need to calculate [tex]x_{n+1}[/tex] are the previous 2 values of x, as well as the acceleration at time [tex]n \cdot \Delta t[/tex], which you calculate from the Forces and masses in the system. In 2D you'd need to calculate the x and ycomponents of the accelerations, then apply the algorithm to find [tex]x_{n+1}[/tex] and [tex]y_{n+1}[/tex] . There might be a common name for this algorithm, if anybody knows it feel free to post the info. 



#22
May308, 09:55 PM

Mentor
P: 11,989

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 RungeKutta 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 3yearold handwritten notes that I have somewhere to figure that out again. 



#23
May408, 12:24 AM

P: 94

Anyone feel like sharing one of their simulators? Would be interesting.




#24
May408, 12:59 AM

P: 4,513

The guys at NASA, NSA, and whoever, should have devoted quit a few manyears 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.




#25
May408, 09:00 AM

Mentor
P: 14,467





#26
May408, 10:42 AM

Mentor
P: 14,467

There are many, many other orbital propagation techniques. One approach is to use Keplerian elements (e.g., [itex]a,e,i,\Omega,\omega,M[/itex]) 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 timevarying 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. 



#27
May408, 10:40 PM

P: 4,513

An interesting possibility to avoid closeapproach 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. 



#28
May508, 10:58 PM

P: 4,513

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




#29
May1008, 11:04 PM

P: 277

I'm also trying to make an nbody orbital simulator, its going quite well  i've recently upgraded from euler to fourth order rungakutta.
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! 



#30
May1508, 11:41 AM

P: 288

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 nearcollision "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 3D, 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. 



#31
May1708, 08:50 AM

P: 516

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?




#32
May1708, 09:43 AM

P: 288

I would suspect that you can't do much better than the fourth order adaptive stepsize Runge Kutta algorithm.




#33
May1708, 11:33 AM

P: 516

Does that mean repeating with a smaller time interval until all 3 steps of 4th order RungeKutta agree within an error margin?




#34
May1708, 12:51 PM

P: 288

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.0e10 both when the force is large and changing quickly AND when the force is next to nothing and staying put?)




#35
May2308, 01:20 PM

Sci Advisor
PF Gold
P: 1,542

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 higherorder 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 pointmass 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 pointmass sources: radiation pressure, relativity, nonspherical 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 


Register to reply 
Related Discussions  
Programming an Orbital Simulator  Programming & Computer Science  16  
Is it possible to create Orbital construction yards for Spaceships?  General Astronomy  6  
Looking for old DOS orbit simulator.  General Discussion  4  
Orbital simulator and non fixed orbit  Introductory Physics Homework  1  
good Flight Simulator  Aerospace Engineering  2 