More precise orbital simulation program

[grav]

I thought this might be the best place for this thread. Actually, I don't know much about how simulation programs are made or how precise they are, but I make a few of my own, and I've noticed that adding a couple of simple lines to the calculations for the gravitation of a two-body system can make them much more precise, up to millions of times as much, so that they can be run faster as well. It works best for masses of similar quantity, though, such as binary star systems and such, while a small mass orbitting a much larger one will remain about the same.

The general program lines for gravity might normally look something like this:

d=((x(z)-x(3-z))^2+(y(z)-y(3-z))^2)^.5
vx(z)=vx(z)-G*M(3-z)*(x(z)-x(3-z))/(d^3)*ti
vy(z)=vy(z)-G*M(3-z)*(y(z)-y(3-z))/(d^3)*ti
x(z)=x(z)+vx(z)*ti
y(z)=y(z)+vy(z)*ti

which is repeated in a loop as the orbits proceed. Here, d is the total distance between the two bodies each time the loop is run through, x and y are the coordinates of each body, vx and vy the vectors of speed in the x and y directions, and z is the designated number for the bodies, 1 and 2, which alternate through these calculations to aquire their orbits. G is the gravitational constant and M is the mass of the body, which is given here as M(3-z) since we are considering the mass of the other body, so for z=1, 3-z=2 and for z=2, 3-z=1.

Okay, now here is what I have found to be a surprisingly precise program for at least a two body system of similar masses:

x(z)=x(z)+vx(z)*ti-G*M(3-z)*(x(z)-x(3-z))/(d^3)*(ti^2)/2*M(z)/(M(z)+M(3-z))
y(z)=y(z)+vy(z)*ti-G*M(3-z)*(y(z)-y(3-z))/(d^3)*(ti^2)/2*M(z)/(M(z)+M(3-z))
d=((x(z)-x(3-z))^2+(y(z)-y(3-z))^2)^.5
vx(z)=vx(z)-G*M(3-z)*(x(z)-x(3-z))/(d^3)*ti
vy(z)=vy(z)-G*M(3-z)*(y(z)-y(3-z))/(d^3)*ti

in exactly that order. I'm not sure how many programmers we have on this forum, but I think one will like the precision this program provides. If anyone else has any more suggestions for creating a faster or more precise simulation program, please let me know.

[tony873004]

Why are you dividing by d^3? Your acceleration formula should be GM/d^2.

Also, if you want to make it faster, you computer probably does d * d much faster than d^2. Try them both and compare.

[grav]

Why are you dividing by d^3? Your acceleration formula should be GM/d^2.

I am using vectors for the acceleration in the x and y directions, so a_x=(GM/d^2)*(x/d)=GMx/d^3 and a_y=(GM/d^2)*(y/d)=GMy/d^3. How do you do it?

Also, if you want to make it faster, you computer probably does d * d much faster than d^2. Try them both and compare.

Yes, that does make it slightly faster also. I've found that out while experimenting with different variations of writing the same program in different ways, along with a few other "tricks", such as only printing text periodically, since that seems to take an incredible amount of time. But while these can increase it somewhat as far as the program itself is concerned, the lines for the program in this thread can make it hundreds to thousands of times faster (for two-body systems of similar mass, anyway) by allowing one to increase the time increments with the same precision, or up to millions of times more precise for the same running time.

[tony873004]

I am using vectors for the acceleration in the x and y directions, so a_x=(GM/d^2)*(x/d)=GMx/d^3 and a_y=(GM/d^2)*(y/d)=GMy/d^3. How do you do it?

The same way you did it, at least in the Euler integrator. I just broke it into 2 steps, that's why I didn't recognize your d^3.

Here's a link to a thread on my forum where the Euler integrator source code appears. I compute all the accelerations and velocities first, then in another loop, I update the positions from the newly-computed velocities.
http://www.orbitsimulator.com/cgi-bi...1160874415/6#6
Also, sqrt(x) is faster than (x)^0.5. These little things make a big difference since they're inside a loop that will be executed thousands to millions of times.

[grav]

Well, I saw the program code you referenced. I like the way you saved time by running the loop for the coordinates separately so that it only needs to run through n times for the number of bodies rather than n^2. The thing is, though, I ran the configuration you are using, and it was way off compared to this one. Unless I put wrote yours incorrectly or something, I think you will be very pleased with the new configuration by itself, much less the improvement in the stability of the orbit. I had forgotten I originally experimented with the configuration (the order in which the calculations take place) to come up with the most stable orbit before I even tried to improve on that. The new lines will all have to be run together through the loop and so might appear that they would take more time, but you could probably run it up to hundreds of times faster and still gain a more stable orbit than that one seems to give, so I think you will like it.

I will rewrite the lines to be used and post it. Just to be sure I get it right, though, I have a couple of questions. How do you figure tM and the object mass. It appears that the increments of time are also figured in, but how exactly? And what is 398600440000000?

And just in case this turns out to be a million dollar deal or something (would be nice ), this thread constitutes a copyright.

[grav]

Just to be sure I get it right, though, I have a couple of questions. How do you figure tM and the object mass. It appears that the increments of time are also figured in, but how exactly? And what is 398600440000000?

Okay, I saw on your forum that 398600440000000 is calculated as G times Earth's mass. But where is the time? I guess it must be figured in separately for the coordinates when setting the graphics, right?

[grav]

Well, Tony, I have to hand it to you. It looks like you've got yourself a very good configuration after all. Perfect, even. I wrote it in wrong. You must have done some experimenting of your own or something. Mine has the same accuracy, but yours is apparently much more stable. Mine starts off just as precise, but then very slowly and steadily increases or decreases in size, so yours is actually better over the long run, and because of the way you've got your loops, faster even using the same increments for time. Apparently, by writing in the change in distance separate from the main loop, not only have you saved time in the program, but you've stopped any possible corruption of the orbits forthright. Guess I won't be making that million anytime soon . So the only thing I've managed to improve upon is my own earlier program, but you've done it one better, so much thanks for that. I'm going to keep trying to improve on it, though. Thanks for the tips.

As far as the timing in your program goes, then, about the only thing that might be improved upon to save more time is just the omission of a few lines, when figuring the acceleration. To run it just a little faster, you might try rewriting it as:

// CodeOptomizer.cpp : Defines the entry point for the DLL application.   
//   
   
#include "stdafx.h"   
   
BOOL APIENTRY DllMain( HANDLE hModule,    
   DWORD  ul_reason_for_call,    
   LPVOID lpReserved   
 )   
{   
    return TRUE;   
}   
   
int _stdcall RunLoop(int NumObjects, double * ObjMass, double * objx, double * objy, double * objz, double * objvx, double * objvy, double * objvz, double * objsize, long * CollisionObjectA, long * CollisionObjectB)    
{   
 if (NumObjects < 1) return 0; //Should be 1 or more   
   
 double tM;   
 double tM2;   
 double dx;   
 double dy;   
 double dz;   
 double D;   
 double f;   
 double fx;   
 double fy;   
 double fz;   
 int Collisions;   
 Collisions = 0;   
 for (int k=1;k<=NumObjects;k++) {   
  tM = ObjMass[k] * 398600440000000;   
  for (int j=k+1;j<=NumObjects;j++) {   
    dx = objx[j] - objx[k];   
    dy = objy[j] - objy[k];   
    dz = objz[j] - objz[k];   
    D = sqrt(dx*dx+dy*dy+dz*dz);   
    if  (((objsize[k] + objsize[j])/2) > D) {   
	Collisions = Collisions + 1;   
	CollisionObjectA[Collisions] = j;   
	CollisionObjectB[Collisions] = k;   
    }   
    if (tM > 0) { 
  f = tM/D/D/D;  
  objvx[j] = objvx[j] - dx*f;   
  objvy[j] = objvy[j] - dy*f;   
  objvz[j] = objvz[j] - dz*f;   
    }   
  tM2 = ObjMass[j] * 398600440000000;   
    if (tM2 > 0) {   
  f = tM2/D/D/D;
  objvx[k] = objvx[k] + dx*f;   
  objvy[k] = objvy[k] + dy*f;   
  objvz[k] = objvz[k] + dz*f;   
    }   
	
  }   
 }   
   
 for (int h=0;h<=NumObjects;h++) {   
  objx[h] = objx[h] + objvx[h];   
  objy[h] = objy[h] + objvy[h];   
  objz[h] = objz[h] + objvz[h];   
 }   
	
	
 //End of loop computations ------   
 return 1; //Finished successfully   
	
}

[tony873004]

You'll have to tell me where you made the change.

The reason you don't see time step in my code is that it is handled external to this code. When a user increases the time step by 2, the program increases the velocities of the objects by 2 and their mass by 2^2. This saves me from having to do the conversion in the loop, hence saving a lot of time.

[tony873004]

Check that you're using double precision variables. I find it odd that you say your program gives the same results for short simulations, but not long ones.

[grav]

You'll have to tell me where you made the change.

I just changed these

f = (1/D) * (1/D) * tM;   
  fx = (dx / D) * f;   
  fy = (dy / D) * f;   
  fz = (dz / D) * f;   
  objvx[j] = objvx[j] - fx;   
  objvy[j] = objvy[j] - fy;   
  objvz[j] = objvz[j] - fz;

and

 f = (1/D) * (1/D) * tM2;   
  fx = (-dx / D) * f;   
  fy = (-dy / D) * f;   
  fz = (-dz / D) * f;   
  objvx[k] = objvx[k] - fx;   
  objvy[k] = objvy[k] - fy;   
  objvz[k] = objvz[k] - fz;

to these.

 f = tM/D/D/D;
  objvx[j] = objvx[j] - dx*f;
  objvy[j] = objvy[j] - dy*f;
  objvz[j] = objvz[j] - dz*f;

and 

 f = tM2/D/D/D;
  objvx[k] = objvx[k] + dx*f;
  objvy[k] = objvy[k] + dy*f;
  objvz[k] = objvz[k] + dy*f;

Check that you're using double precision variables. I find it odd that you say your program gives the same results for short simulations, but not long ones.

I'm using UBasic, which gives me up to 2600 digits, but the more I use, the slower it runs. I generally run it to about thirty or fourty decimal places. I think the reason my program runs off like that is because I was finding for the coordinates within the main loop, instead of doing it separately like you are, as well as running a loop for each body individually, so the slight change in coordinates of one body with each run affects the next.

[HenrikOlsen]

Ah, you've just explained what I couldn't understand.

I got confused when one of you wrote that x*x is faster that x^2, when all the compilers I know would result in the same calculation for both, (unless you force 2 to be real or use an explicit power function)

Try to learn a compiled language, you'll be amazed that the speed increase, and possibly confounded by the lack of effect many of the tweaks such as the ones you're doing has, since the compiler does most of them already.

[grav]

Ah, you've just explained what I couldn't understand.

I got confused when one of you wrote that x*x is faster that x^2, when all the compilers I know would result in the same calculation for both, (unless you force 2 to be real or use an explicit power function)

Try to learn a compiled language, you'll be amazed that the speed increase, and possibly confounded by the lack of effect many of the tweaks such as the ones you're doing has, since the compiler does most of them already.

Thanks. I had been looking around for another program that might run faster and/or more accurately, but I was experimenting with different variations of formulas for gravity, and Basic is all I know, so I was hoping I could just change the few lines of calculations for gravity from another program to my needs, which wouldn't require much language skill, except for the layout I prefer.

[grav]

Well, it looks like changing the orientation of the calculations only exchanged one problem for another. Instead of the size of the orbits steadily increasing or decreasing with time, they now steadily precess counter to their lines of motion around each other, although the size of the orbits themselves don't change. It seems like there might be some simple way to help resolve one or the other of these.

[tony873004]

...Instead of the size of the orbits steadily increasing or decreasing with time, they now steadily precess counter to their lines of motion around each other, although the size of the orbits themselves don't change...

My code does that if you use too large of a time step. Take Mercury for example. If you use a time step like 64K seconds, Mercury's orbit noticabely precesses. And it's precessing backwards from the direction it precesses in real life. If you slow the time step down, the precession gets less and less. If you slow it below 2k seconds, you've slowed it past 0, and it is now free to precess in the natural direction. As you slow it further, the rate of precession asymptotically approaches the correct Newtonian-predicted value of about 530 arcseconds per century. Here's a graph I made of time step vs precession in arcseconds per century.


Try running it to 15 decimal places. I doubt you need much more than that. 15 is what my program uses (Double precision). With a slow enough time step, if you propogate Apophis' orbit from today's vectors, to its 2029 approach of Earth, it only misses JPL's predicted close approach distance by about 200 km. That's not too bad considering both Apophis and Earth travelled over 20 billion kilometers each, on completely separate trajectories during this 22 year period to get to their encounter. For comparison, if someone could hit a golf ball from San Francisco to Tijuana, Mexico, this would be like predicting within 1/10th of a millimeter where it would land. This shows you just how powerful a Newtonian-only model, using a simple Euler algorithim with double precision variables and a slow time step can be.

I suspect that the 200 km difference has little to do with my use of only 15 digits. Their simulator takes GR into account, and perhaps more gravitational perturbers from the asteroid belt, as well as non-spherical gravity.

Since you say that UBasic lets you choose how many digits you use, if you get your code working, I'd be interested in seeing a plot of how different number of digits effects the results. How many are required for good results? How many are overkill?

[publius]

Tony,

Is the source code for your simulator available? If it's "trade secrets", then no problem.

The reason I ask is I'd like to run it through Intel's C++ compiler and compare it. That sucker does vectorizing and parallelizing for multi-CPUs (and multi-core and hyperthreading) and I always like to see if it's as good as they like to brag about. It has impressed me with other stuff.

-Richard

[tony873004]

Tony,

Is the source code for your simulator available? If it's "trade secrets", then no problem.

The reason I ask is I'd like to run it through Intel's C++ compiler and compare it. That sucker does vectorizing and parallelizing for multi-CPUs (and multi-core and hyperthreading) and I always like to see if it's as good as they like to brag about. It has impressed me with other stuff.

-Richard

Gravity Simulator is written in Visual Basic 6.0, with the exception of the number-crunching routine. It's written in C++ since C++ is much faster than VB 6.0, and this is the only part of the program where speed is advantagous.

The C++ code is compiled into a .dll. The VB code calls this routine once per iteration. I provided a link to the C++ code above. No trade secrets here. It's just a basic Euler's method propogating formulas from your High School Physics book. I guess it would be as easy as copying it, pasting it, and compiling it into a .dll using Intel's compiler. I used Visual C++.Net.

Assuming that works, and assuming you have Gravity Simulator installed on your computer, you simply need to replace the file "codeoptomizer.dll" with your replacement .dll. Just rename the old one incase you want to restore it. Then to test speeds, you want to run Gravity Simulator with the graphics turned off. This forces the program to only crunch the numbers.

I'd be very interested to see if a different compiler would make a big difference. So if you attempt this, let me know.

BTW, if you want to see a comparison between the speed of VB 6.0 and .NET compiled C++, simply rename codeoptomizer.dll, so VB can't find it. This forces it to use VB to do the math. I believe the difference in speed is a factor of 6.

[publius]

Tony,

I see. I doubt the Intel compiler would do much better on that simple function in the .dll that MS's now. What it would shine at if the actual time-stepping, the big loop was done as well. I don't know see enough if the above function to do much.

If the whole thing were there, Intel's compiler has vectorizing and parallelizing, plus profile guided optimization. The latter is something else. You do an "instrumented build", where the compiler puts metric code in that keeps track of the data flow and stores the results in a data file. You turn that instrumented build loose on typical input data conditions. It's slow as heck, but it is carefully keeping track. Now, the compiler then takes that data set to better optimize the final build.

That can do wonders. But to do that here, the Intel compiler would need to see the whole picture, with the main time step loop in there as well. And it would need to know the better data picture of where the inputs to the loop are stored, etc, etc.

-Richard

[tony873004]

...To run it just a little faster, you might try rewriting it as...

Thanks. That make it 8% faster.

[grav]

Since you say that UBasic lets you choose how many digits you use, if you get your code working, I'd be interested in seeing a plot of how different number of digits effects the results. How many are required for good results? How many are overkill?

I've run the program for 9, 14, 19, and 24 digits in the decimal place, which are point 2, 3, 4, and 5 in the system. The integer part can be as great as necessary up to the maximum 2600 allowed, though, so the total number of digits altogether varies. I ran them in intervals of 10000 to 100000 seconds per loop. Apparently, 9 decimal places doesn't even put the orbit on the map. The program won't even begin to work, so I can see why you would need to double that. With fourteen decimal places, it works fine. But it also works precisely the same at 19 and 24. So the minimum number of decimal places required would be at least 10, but more than 14 makes no difference, so the 15 you use is plenty, unless 15 is the total number of digits including the integer part, in which case I couldn't be sure. I use measurements in units of seconds and meters, so the integer part would be the greatest portion for these. It seems I've actually been using 24 for the fraction part myself (point 5), so I can probably double my speed by using 14 with the same performance.

If you want to compare my optimum results to yours, I used the orbit of Mercury around the sun only, as a two body system. For intervals of 10000 seconds, I get a negative precession of -10000 arc-seconds per century. For 25000 seconds per loop, I get -64000 arc-seconds. For 50000, I get -260000, and for 100000, I get -1000000.

[grav]

If you want to compare my optimum results to yours, I used the orbit of Mercury around the sun only, as a two body system. For intervals of 10000 seconds, I get a negative precession of -10000 arc-seconds per century. For 25000 seconds per loop, I get -64000 arc-seconds. For 50000, I get -260000, and for 100000, I get -1000000.

I've noticed that pattern runs fairly steadily at prec=-(ti^2/10000) arc-seconds per century, where 'ti' is the time interval of each loop. I think I'll see if I can't pick this apart to see where it's coming from and try to come up with something just as steady to counterbalance it.

[grav]

Okay. I've tried running it for a perfectly circular orbit to see what I could come up with. But instead of circular as it should be, though, the orbit becomes slightly elliptical and the difference in distance between that of the closest approach and furthest separation between the bodies is precisely equal to vo*ti, where 'vo' is the sum of the original velocities of travel of the bodies, parallel and in opposite directions from each other. Now we're getting somewhere. All I have to do now is to figure out some way to add that distance back into the formulas to provide a more precise and stable orbit.

[tony873004]

...All I have to do now is to figure out some way to add that distance back into the formulas to provide a more precise and stable orbit.

You'll be beating yourself up with methods like this. Your next logical step is to code a Runge-Kutta 4 routine.

[grav]

You'll be beating yourself up with methods like this. Your next logical step is to code a Runge-Kutta 4 routine.

You got it. I've been trying to see where the coordinates according to the program and that of a perfectly circular orbit differ, so that I can find a formula to better orders of approximation. But I'll have to do it in such a way that it works out for elliptical orbits, straight line paths, etc, as well. I'll look into the Runge-Kutta 4 routine.

For a circular orbit, the first run through the loop seems to fall off by a factor of (1-3*GM/d^3*ti^2/2) for x and (1-GM/d^3*ti^2/2) for y. In other words, in factors of a*ti^2/(2d). I'll have to run it at other original angles as well, such as 45 degrees, where vx and vy are the same and see how they differ at each angle around the orbit.

For a straight line, the correct formulas are v'=v+a*ti for speed and x'=x+v*t+a*t^2/2 for distance with constant acceleration. For each loop, the program would run it at

1) v'=v+a*ti
x'=x+v'*ti=x+(v+a*ti)*ti=x+v*ti+a*ti^2

2)v''=v'+a*ti=(v+a*ti)+a*ti=v+2a*ti
x''=x'+v''*ti=(x+v*ti+a*ti^2)+(v+2a*ti)*ti=x+2v*ti +3a*ti^2

3)v'''=v''+a*ti=(v+2a*ti)+a*ti=v+3*a*ti
x'''=x''+v'''*ti=(x+2v*ti+3a*ti^2)+(v+3*a*ti)*ti=x +3v*ti+6a*ti^2

So for 'n' number of runs through the loop, we get

v(n)=v+n*a*ti
x(n)=x+n*v*ti+[n*(n+1)/2]*a*ti^2

Since the total time 't' is t=n*ti, this becomes

v(n)=v+a*t
x(n)=x+v*t+a*t^2/2+a*t*ti/2

v(n) is correct for a smooth transition, but x(n) is not. It is off by an additional a*t*ti/2. If we figure that over 'n' number of loops, then each loop runs off by a*ti^2/2, same as that for the circular orbit. Now, that is how I made mine more precise to begin with, by running it to this next order of approximation, subtracting a*ti^2/2 from the distances with each run through the loop, along with a ratio of the masses, making it millions of times more precise. But your configuration already seems to account for this for the most part, by figuring the distances separately from the main loop, causing only the precession, but for the same reason, so I'm still trying to edge the difference of a*ti^2/2 in there somewhere. I'll keep working on it.

[Flaney]

Go here for adaptive stepsize control for RK

http://www.tat.physik.uni-tuebingen....mrec/f16-2.pdf

Joy

[publius]

Tony,

You didn't tell me the main program would simply execute the calculation itself if the "codeoptimizer.dll" failed to load. I compiled the above snippet, and thought it was working, but slower......... Just remembered something about MS linker/compiler behavior. First I compiled it as C++ code, and so the exported name was a decorated C++ name. Fixed that, but then rediscovered the odd MS behavior that if you a "declspec" a C stdcall name, it will still be decorated, "ie _runloop@12". But if you export it with a .def file, you get the undecorated name, 'runloop'. Sigh...........

I'll tell you what the results are when I get it working. Intel and MS base code look the same, except if I go to SSE2 instructions. Both compilers will use that for double precision. I don't know if the Intel compiler really has much to vectorize.

-Richard

[tony873004]

Tony,

You didn't tell me the main program would simply execute the calculation itself if the "codeoptimizer.dll" failed to load...

Yes I did

...BTW, if you want to see a comparison between the speed of VB 6.0 and .NET compiled C++, simply rename codeoptomizer.dll, so VB can't find it. This forces it to use VB to do the math...

But I've just discovered in most circumstances that the C++ dll doesn't speed things up at all, and on some computers it actually slows it down. So I retract my claim of 6x faster that I made above. (One of the biggest hassels of programming is assuming that something will behave on other computers the same way it behaves in my computer.) The only way to see for yourself is to temporarily re-name codeoptomizer.dll and compare the speed of a simulation to its speed using codepotomizer.dll. If you re-compile it using a different compiler, you can probably guage the relative difference between your compiled version and the .NET compiled version.

If it helps you, here's the VB code that calls the .dll.
This declares it. This snipet of code is only executed once:

Private Declare Function RunLoop Lib "CodeOptomizer.dll" (ByVal NumObjects As Long, _
ByRef ObjMass As Double, ByRef objx As Double, ByRef objy As Double, _
ByRef objz As Double, ByRef objvx As Double, ByRef objvy As Double, _
ByRef objvz As Double, ByRef objsize As Double, _
ByRef CollisionObjectA As Long, ByRef CollisionObjectB As Long) As Integer

And this calls it. This code is executed once per iteration.

xxx = RunLoop(objects, ObjMass(0), objx(0), objy(0), objz(0), objvx(0), objvy(0), objvz(0), objsize(0), CollisionObjectA(0), CollisionObjectB(0))

Also, if you make a change to codeoptomizer.dll, either by replacing it with a different version, or by renaming it as to hide it, you must exit Gravity Simulator and start again as it only reads it in once at the beginning of the program, and won't recognize your change.

[publius]

Tony,

Indeed you did tell me. Doh! Read the whole thing, Richard, read the whole thing carefully.............. More later when I get around to it.....


-Richard

[HenrikOlsen]

Is it actually truly called CodeOptomizer?

Reminds me on a time I had to support code originally written by a guy who is severely dyslexic, all the function names where not quite written in English.

[grav]

...To run it just a little faster, you might try rewriting it as...

Thanks. That make it 8% faster.

If you liked that little bit, you should really like this. It has nothing to do with what I've been working on, but it should more than double your speed, depending on the time consumption of the rest of the program (the more time the rest of the program takes up, the more this will help), with the same accuracy, or you can run it at half of the time intervals for the same speed, but with twice the accuracy. If you're lucky enough that it quadruples it, you can have twice the accuracy and twice the speed. Just type in the three lines in blue. You can experiment with the value of 'qw' to find the point where the speed tops off. But with qw=10, it'll only plot the points one out of every ten times, so you don't want to go over qw=1000 or so to keep a smooth progression in the graphics. Unfortunately, however, the speed tops off for me at about qw=5 at about twice as fast, so I don't have the luxury of worrying about it . I have to go as high as qw=100 just to the speed up from there to about 2.5 . Maybe yours will do better.

// CodeOptomizer.cpp : Defines the entry point for the DLL application.   
//   
   
#include "stdafx.h"   
   
BOOL APIENTRY DllMain( HANDLE hModule,    
   DWORD  ul_reason_for_call,    
   LPVOID lpReserved   
 )   
{   
    return TRUE;   
}   
   
int _stdcall RunLoop(int NumObjects, double * ObjMass, double * objx, double * objy, double * objz, double * objvx, double * objvy, double * objvz, double * objsize, long * CollisionObjectA, long * CollisionObjectB)    
{   
 if (NumObjects < 1) return 0; //Should be 1 or more   
   
 double tM;   
 double tM2;   
 double dx;   
 double dy;   
 double dz;   
 double D;   
 double f;   
 double fx;   
 double fy;   
 double fz;   
 int Collisions;   
 Collisions = 0;
 qw=10;
 for (int qq=1;qq<=qw;qq++) {   
 for (int k=1;k<=NumObjects;k++) {   
  tM = ObjMass[k] * 398600440000000;   
  for (int j=k+1;j<=NumObjects;j++) {   
    dx = objx[j] - objx[k];   
    dy = objy[j] - objy[k];   
    dz = objz[j] - objz[k];   
    D = sqrt(dx*dx+dy*dy+dz*dz);   
    if  (((objsize[k] + objsize[j])/2) > D) {   
	Collisions = Collisions + 1;   
	CollisionObjectA[Collisions] = j;   
	CollisionObjectB[Collisions] = k;   
    }   
    if (tM > 0) { 
  f = tM/D/D/D;  
  objvx[j] = objvx[j] - dx*f;   
  objvy[j] = objvy[j] - dy*f;   
  objvz[j] = objvz[j] - dz*f;   
    }   
  tM2 = ObjMass[j] * 398600440000000;   
    if (tM2 > 0) {   
  f = tM2/D/D/D;
  objvx[k] = objvx[k] + dx*f;   
  objvy[k] = objvy[k] + dy*f;   
  objvz[k] = objvz[k] + dz*f;   
    }   
	
  }   
 }   
   
 for (int h=0;h<=NumObjects;h++) {   
  objx[h] = objx[h] + objvx[h];   
  objy[h] = objy[h] + objvy[h];   
  objz[h] = objz[h] + objvz[h];   
 }   
   }	
	
 //End of loop computations ------   
 return 1; //Finished successfully   
	
}

[tony873004]

Thanks, I'll try this later today when I get home. I already added your previous suggestion (8%, which sometimes is as much as 14%) to the lastet beta version. Together with other changes, it's now 60% faster. The link is in my web forum. I'm anxious to try your new idea.