# Thread: More precise orbital simulation program

1. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Grav,

Sorry, but your 'qq' correction makes Mercury's orbit precess like a mad hoola-hoop around the sun at a 64K timestep. Worse precession yet. It may be the time-step needs to be adjusted in there.

Tony is simply carrying the timestep with the masses and velocity. Figure out what your 'qq' looks like with the time-step multiplied through.

-Richard

2. Originally Posted by publius
BTW, what exactly is the "qq" correction above? Is that the "better ellipse" thing? Or something else? I was going to sit down, if I felt like it, and play with some higher derivatives -- that may be the essence of what 'qq' is, at least assuming the curve is elliptical.

-Richard
It's not the better ellipse, but it is based upon the same concept. Since UBasic will run up to 2600 digits, I've been running them out to about 400 at very small time intervals, about 10^-40 seconds, and seeing what the difference is for running the loop through once at 10^-36 seconds and a million times at 10^-40. Since the smaller time interval is more precise, I have been trying to adjust the calculations accordingly, to get the same result while running through the loop just once instead of a million times. With this many digits, one can see a pattern forming and can kind of pick out the calculations to be used. I can tell they're on the order of [G*(m1+m2)*ti*ti/d/d/d]^n for consecutive orders of approximation, but I'm still having trouble getting the coefficients right. Here's a sample of a typical run.

Code:
```?vx(1)
-4999.99999999999999999999999999999999999999999999999999999999999999999999999999
50000049999999999999999999999999999999999999999999999999999999999999999999999933
70836591673370829999999999999999999999999999999999999999999999999999999999990554
13543207131309950554229849422999999999999999999999999999999999999999999998608841
25644813184545704824757430760808465673931249999999999999999999999999999791188147
00475333911252154742256531145147902781918850597994216562499999999999968268978937
7657936
?vy(1)
0.00000000000000000000000000000000010000000000000000000000000000000000000000000
00000000000000000000000000000000000661666674165997500000000000000000000000000000
00000000000000000000000000000000072237084573792672915416760099375000000000000000
00000000000000000000000000000009062702637379254891319449962062293172657832812500
00000000000000000000000000001214342881698955396469123581177089028406309980018862
21493667968750000000000000169082550020736959321544005215654776002472460542345525
1743
?x(1)
-0.00000000000000000000000000000049999999999999999999999999999999999999999999999
99999999999999999999999999999998333333333334999999999999999999999999999999999999
99999999999999999999999999999867416416666832083583333300499999999999999999999999
99999999999999999999999999986505861361126723892430562115037874991321874999999999
99999999999999999999999998454261166286967128689731285668434446388181258148015871
56249999999999999999999810170003741383589568329061228804145899783425748677744756
7473
?y(1)
-4999999999.99999999999999999999999999999999999999999999999999999999999999999999
49999949999999999999999999999999999999999999999999999999999999999999999999999983
45830000001654169999999999999999999999999999999999999999999999999999999999998796
04497444551405342222235656826999999999999999999999999999999999999999999999886715
76352785941501302496790882734311929372943749999999999999999999999999999987856510
42751612969450966388100950437700780765073235712076486562499999999999998590970291
3516005103101```

3. Originally Posted by publius
Grav,

Sorry, but your 'qq' correction makes Mercury's orbit precess like a mad hoola-hoop around the sun at a 64K timestep. Worse precession yet. It may be the time-step needs to be adjusted in there.

Tony is simply carrying the timestep with the masses and velocity. Figure out what your 'qq' looks like with the time-step multiplied through.

-Richard
Is there a variable originally used to store the time-step interval somewhere in the program?

4. Established Member
Join Date
Mar 2005
Posts
1,249
Originally Posted by grav
Is there a variable originally used to store the time-step interval somewhere in the program?
Yes, its called "timestep". See post #108. But you don't have access to it through the .dll.

Sorry I haven't been following this thread too closely lately. I've got lots of less interesting stuff keeping me busy.

5. I guess I'll try to work out the units for mass and velocity in my program the way Tony has his where the intervals of time are incorporated and set the program to "point 9" which gives me 14 digits after the decimal place. If I get it to work for that, it should work fine for the program you are using.

6. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
You know, I've been looking at Mercury's (truly anomalous ) precession here in Gravity Simulator. Heck, there is not that much difference between the original routine, my original "average velocity) correction to the position update, and Grav's latest "qq" thing.

Start the thing out with the "onlyplanets" dataset, and zoom in on the sun to about 1 - 2AU scale so Mercury's orbit looms large. Now, let it plot. At a 64K timestep, that sucker "smears" out quite a donut shape due to the precession.

As near as I can figure, at 64K, all three sweep out that full donut, meaning a full circle of precession by about the year 2300. Drop that down by half just to 32K, and that settles that mess down considerably. It's still precessing, and in the wrong direction, but it settles it down greatly.

-Richard

7. Originally Posted by publius
You know, I've been looking at Mercury's (truly anomalous ) precession here in Gravity Simulator. Heck, there is not that much difference between the original routine, my original "average velocity) correction to the position update, and Grav's latest "qq" thing.

Start the thing out with the "onlyplanets" dataset, and zoom in on the sun to about 1 - 2AU scale so Mercury's orbit looms large. Now, let it plot. At a 64K timestep, that sucker "smears" out quite a donut shape due to the precession.

As near as I can figure, at 64K, all three sweep out that full donut, meaning a full circle of precession by about the year 2300. Drop that down by half just to 32K, and that settles that mess down considerably. It's still precessing, and in the wrong direction, but it settles it down greatly.

-Richard
Oh, you'll notice a big difference with the qq. Have you tried f3=(f1+f2)*ti in the code, where 'ti' is whatever variable for the time-step interval being used?

8. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
I confirmed something. The precession seems to be going roughly
as h^2. At 64K it takes 300 years to sweep out the donut. At 32K, it takes it to the year ~3150 to make the full donut, or 1200 years. That's 4x the time.

-Richard

9. Well, I haven't finished reprogramming the "Tony" way yet, but apparently the units are correct for 'qq' they way are now in the code I posted. It must be the number of digits that's the problem, then. 'qq' must not be getting up much past 1 or something. I might have to multiply and then divide back out, as much as I know you can't stand that. How many digits can you have in your numbers? Is there a limit to how many actual digits can be included in the integer part of it? What about the integer and the fraction part together?

10. Established Member
Join Date
Mar 2005
Posts
1,249
Originally Posted by publius
I confirmed something. The precession seems to be going roughly
as h^2. At 64K it takes 300 years to sweep out the donut. At 32K, it takes it to the year ~3150 to make the full donut, or 1200 years. That's 4x the time.

-Richard
Yes, it's definately not linear. But at some point your formula will break down. At 2K the precession halts and changes direction to the correct direction, and by around timestep= 32 seconds Mercury will precess very close to its Newtonian-predicted rate of 531 arcseconds per century.

BTW... if you want custom time steps, so you don't have to choose powers of 2 (1, 2, 4, 8...) you can use the autopilot to set it.

Also, if you change the graphics interval in the Preferences menu, it plots less often and drastically speeds things up (beta version only, but I think you're running that)

11. Established Member
Join Date
Mar 2005
Posts
1,249
Originally Posted by grav
Well, I haven't finished reprogramming the "Tony" way yet...
I'm famous

12. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Originally Posted by grav
How many digits can you have in your numbers? Is there a limit to how many actual digits can be included in the integer part of it? What about the integer and the fraction part together?
Grav,

We're running on your typical Intel x86 hardware, which uses the
IEEE-whateverNumber floating point standard. Floating point is done in hardware. Double precision is a 64-bit floating point format. These are *binary* floating point representations (ie 1.mantissa * 2^exponent), so the number of decimal digits doesn't work out exactly. But roughly, double precision has 15 decimal digits, with an exponent range of 10^(+/- 308).

Now, internally, the x87 standard uses 80 bits, ("extended double"). However, the SSE instructions do the arithmetic at the storage precision.

My Intel FORTRAN compiler implements a quadruple precision, 128 bits in software (very slow -- well, what "software" means is you, typically using the integer registers, "manually" do the bit operations required to do the arthimetic), giving effectively 30 digits of decimal precision.

-Richard

13. Well, I don't get it. I just programmed it in the way Tony does it (which is very efficient, by the way) and it worked just fine, great in fact. Even at 64K there is no precession. I don't start to see a change until I run it at about 200K or so, and barely then. The only extra thing I'm doing is multiplying by 'qq', and that is always just above 1, with only 14 digits past the decimal place, for 15 digits in all, so that multiplied by any value, regardless of the number of digits used for the other number, should increase it by that much all in all. I don't see any way this could not be working for you. Could you copy and paste the code you used for this?

14. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Grav,

Just what is a typical value for your 'qq'? I used the exact formula you gave. By "14 digits" do you mean qq is typically this:

1.0000....1, that is 13 or so zeros before anything? If so, that's the problem. qq, and probably the terms building are *so close to 1* for 15 digits of precision, that the result is just 1.

And really I'm surprised a correction that small would make any difference at all, even using all the digits you have.

-Richard

15. Originally Posted by publius
Grav,

Just what is a typical value for your 'qq'? I used the exact formula you gave. By "14 digits" do you mean qq is typically this:

1.0000....1, that is 13 or so zeros before anything? If so, that's the problem. qq, and probably the terms building are *so close to 1* for 15 digits of precision, that the result is just 1.

And really I'm surprised a correction that small would make any difference at all, even using all the digits you have.

-Richard
Okay, I added a couple of lines to find the highest and lowest values for qq as it ran through the loops. At 64K, the lowest value is 1.00026535443448 and the highest is 1.00092731832708.

16. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Originally Posted by grav
Okay, I added a couple of lines to find the highest and lowest values for qq as it ran through the loops. At 64K, the lowest value is 1.00026535443448 and the highest is 1.00092731832708.
Well, that ought to do something, then. It's getting late, and I'm going to pack in. What I need to do is run the thing in the debugger where I can watch the variables. The trouble is doing it from a .dll. Tony's VB code loads and calls this as a .dll, and getting the debugger "attached" to that will require me to read some docs.................:sigh: If it's possible. It ought to be, but, well, MS tools have evolved a lot since I was on top of them, and this old dog can't learn new tricks as fast as he once could.

-Richard

17. Thanks, publius. Just so you know, it's also important that the code is written in exactly the order I wrote it. I had to rearrange a couple of things the way it needed to be for this and it must be run in precisely that order.

18. The object masses used in the program you guys are using aren't negative to account for G acting in a negative direction or anything, are they? If so, then the bodies originally move in the opposite direction than what I have, and that would actually give twice the precession. You would have to use qq=1 - f3/6 + f3*f3/15 - f3*f3*f3/28 in that case to account for the negative values of the masses.

19. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Well, I got the debugger running (after reading a bunch of crap).

The value of 'qq' is indeed "1.000000000000000" most of the time, save for a few big masses, such as Mr. Sun and Jove.

Also, with the debugger going, I noticed the main program, gravitysimulator.exe was throwing "Inexact Floating Point Result" exceptions every so often. That is in the VB code, not the .dll. Whether that has anything to do with anything, or is just superfluous I don't know. The exception is just silently handled internally by the run-time code.

That may be nothing -- you can set the FPU to complain about everything, even when it rounds, which is normal, and it may be the VB run-time is doing that. There weren't too many of those during the run (100 or so, from the looks), but that does trigger time-consuming exception logic (with all that stack unwinding crap):

" First-chance exception at 0x7c812a5b in gravitysimulator.exe: 0xC000008F: Floating-point inexact result."

-Richard

20. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
That is nothing to worry about all. There's something about some internal VB run-time stuff that does that. That exception basically means "I had to round off the result", which happens 99.9% of the time anyway. That exception should be turned off, and it usually is, but something in the VB run-time does it.

-Richard

21. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Grav,

That did it! That stopped the precession. I figured out why 'qq' was always one. I had put the expression in this form:

qq = 1 + f3*(1/6 + f3*(1/15 + f3/28) )

Note the constants in that expression are interpreted by the compiler as *integers*. Changing the expression to this fixed it:

qq = 1.0L + f3*(1.0L/6.0L + f3*(1.0L/15.0L + f3/28.0L) )

That makes every constant explicitly floating point double. What happened is the compiler interpreted "1/6" and "1/15" as integer expression, which of course, are *zero*. The result was only the f3^3 term contributed anything because the others were mulitplied by zero. :sigh:

That's why it is important to explicitly make constants floating point, because the compiler, following its rules to a T, can do things unexpected.

-Richard

22. Originally Posted by publius
The value of 'qq' is indeed "1.000000000000000" most of the time, save for a few big masses, such as Mr. Sun and Jove.
Well, that's true, it makes with the biggest difference with the largest masses, such as the sun, but the planets wouldn't affect each other much unless they were very close together, in which case this would correct for that as well. It works, I'm tellin' ya , very well. It stabilizes the orbits almost completely. What are you getting for qq for the sun and Mercury?

23. Originally Posted by publius
Grav,

That did it! That stopped the precession. I figured out why 'qq' was always one. I had put the expression in this form:

qq = 1 + f3*(1/6 + f3*(1/15 + f3/28) )

Note the constants in that expression are interpreted by the compiler as *integers*. Changing the expression to this fixed it:

qq = 1.0L + f3*(1.0L/6.0L + f3*(1.0L/15.0L + f3/28.0L) )

That makes every constant explicitly floating point double. What happened is the compiler interpreted "1/6" and "1/15" as integer expression, which of course, are *zero*. The result was only the f3^3 term contributed anything because the others were mulitplied by zero. :sigh:

That's why it is important to explicitly make constants floating point, because the compiler, following its rules to a T, can do things unexpected.

-Richard
Okay, good. Glad to hear it. Whew, now I can relax again.

24. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Grav,

I'm a bit put out with myself that I didn't notice that right off the bat. I should know that, and I did, just forgot. I didn't notice it until I was running it in the debugger, saw qq was smaller than it should be, did it with a calculator to confirm, then expanded the view to the machine code and explicitly saw what was happenning.

You may also have noticed I rearranged your qq expression a bit. Why did I do that?

How many floating point operations are in the following expression?:

x = a*b + a*c

Now, how many are in this one?

x = a*(b + c)

Floating point operations are expensive, even in hardware, although hardware is much faster than software emulation. I had a course in numerical analysis years ago, and I had stuff like this drummed in my little head.

A good optimizing compiler will re-order operations like the above, but that depends on how good the compiler authors were about programming it to recognize when it can do stuff like this.

This is all part of the "smart compiler, dumb programmer" mentality we have today. Well, knowing stuff like this is actually the difference between a "programmer", and a mere "code writer". It is apparently cheaper to spend money on the talent to write smart compilers, and then hire cheap "coders".....

-Richard

25. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Incidently, it was my rearranging to save operations that triggered that integer problem anyway. Had I left it in your form, each term would've been floating point. Had the compiler rearranged it itself during optimizing, it would've then known it had to keep results floating point.

But when I rearranged it, it assumed I meant what I typed, integer constants. So, yep, the compiler would've been smarter than me.

-Richard

26. Originally Posted by publius
You may also have noticed I rearranged your qq expression a bit.
Yes, I noticed that. It looks much more efficient now. As a matter of fact, though, that last term can probably be dropped. qq=1+f3/6 is definite, but I can't make out the next term completely. It could be f3^2/15 or f3^2/14, which kind of negates the last one altogether until I can work it out more precisely. qq=1 + f3/6 + f3^2/14 might even work better anyway, since it would account for some of the rest of the terms left off for larger time-steps, and smaller time-steps would be precise enough anyway with the first couple of terms.

27. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Grav,

I ran it at 64K time step for 2500 years, graphics on. Mercury was very tight -- no smear at all. I could see the slightest little bit, and it is still going in the wrong direction, but only the slightest little bit.

And I must say I'm very impressed that you can pull out correction terms from looking at "digit fields" like that. I think I'm going to call you "brute force". You've got talent, young man, I just wish you'd polish it by learning calculus, differential equations, vector calculus, etc. You could run with the Big Boys then, big time in the Big Boy League shore 'nuff.

However, as impressed as I am with your brute force correction, are we loosing the physics in that brute force correction. Well, you were going for what the result of a small time step would be, but I like to know what that correction really was saying in terms of the actual solution. It might be close, or it might not.

-Richard

28. Originally Posted by publius
Grav,

I ran it at 64K time step for 2500 years, graphics on. Mercury was very tight -- no smear at all. I could see the slightest little bit, and it is still going in the wrong direction, but only the slightest little bit.
Yes, I'm still working on the precise formulation for the terms of that, so it can be brought out further.

And I must say I'm very impressed that you can pull out correction terms from looking at "digit fields" like that. I think I'm going to call you "brute force". You've got talent, young man, I just wish you'd polish it by learning calculus, differential equations, vector calculus, etc. You could run with the Big Boys then, big time in the Big Boy League shore 'nuff.
Thanks. I guess brute force is my preference, but I'll keep studying other methods as well.

However, as impressed as I am with your brute force correction, are we loosing the physics in that brute force correction. Well, you were going for what the result of a small time step would be, but I like to know what that correction really was saying in terms of the actual solution. It might be close, or it might not.
The factor 'F3' varies with ti^2, so a small time-step affects the calculations very little anyway, but will still add some correction. I've been trying to get the path of the orbit in the original program a little more precise, and qq seems to help this just a little bit, becoming just slightly more precise than the original program as far as this goes too, as far as I can tell, so it can be run at larger time steps as well without going off track.

29. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Here is the final version of grav's 'qq' code. Note another thing I did was the distance cubed thing. We first have to calculate D^2 by summing squares. We want the cube. It is a waste of operations to calculate that square root and then cube it in the divisor when we already have its square.

All we have to do is multiply the square root by D^2 to get D^3, which is D3 in the code below -- I also used multiple assigments on one line, something that C allows. "=" is umambiguosly an assigment operater, the comparison operator is "==", so there is no possible confusion.

Code:
```// CodeOpGrav.cpp : N-body loop for Gravity Simulator
// Use Grav's "qq" thing

#include "stdafx.h"
#define MSIZE 800
BOOL APIENTRY DllMain( HANDLE hModule,
DWORD  ul_reason_for_call,
LPVOID lpReserved
)
{
return TRUE;
}

__declspec(dllexport) int _stdcall RunLoop(int NumObjects, double * restrict ObjMass, double * restrict objx, double * restrict objy, double * restrict objz, double * restrict objvx, double * restrict objvy, double * restrict objvz, double * restrict objsize, long * restrict CollisionObjectA, long * restrict CollisionObjectB)
{
if (NumObjects < 1) return 0; //Should be 1 or more

// static double tM[MSIZE],tM2[MSIZE];
static double tM,tM2;
// static double dx[MSIZE], dy[MSIZE], dz[MSIZE];
static double dx, dy, dz;
// static double D[MSIZE], f[MSIZE], D2[MSIZE];
static double D, f1, f2, f3, D2, D3;
static double objvxO[MSIZE], objvyO[MSIZE], objvzO[MSIZE];
double fx;
double fy;
double fz;
double qq;
int Collisions;
Collisions = 0;
for (int k=1;k<=NumObjects;k++) {
tM = ObjMass[k] * 398600440000000;
objvxO[k] = objvx[k];
objvyO[k] = objvy[k];
objvzO[k] = objvz[k];
for (int j=k+1;j<=NumObjects;j++) {

dx = objx[j] - objx[k];
dy = objy[j] - objy[k];
dz = objz[j] - objz[k];
D2 = dx*dx+dy*dy+dz*dz;
D3 = (D = sqrt(D2) )* D2;

if  (((objsize[k] + objsize[j])/2) > D) {
Collisions = Collisions + 1;
CollisionObjectA[Collisions] = j;
CollisionObjectB[Collisions] = k;
}
tM2 = ObjMass[j] * 398600440000000;
f1 = tM/D3;
f2 = tM2/D3;
f3 = f1 + f2;
qq = 1.0L + f3*(1.0L/6.0L + f3*(1.0L/15.0L + f3/28.0L) );
if (tM > 0) {
objvx[j] = objvx[j] - dx*f1*qq;
objvy[j] = objvy[j] - dy*f1*qq;
objvz[j] = objvz[j] - dz*f1*qq;
}
if (tM2 > 0) {
objvx[k] = objvx[k] + dx*f2*qq;
objvy[k] = objvy[k] + dy*f2*qq;
objvz[k] = objvz[k] + dz*f2*qq;
}

}

}
//#pragma vector always

for (int h=1;h<=NumObjects;h++) {
objx[h] = objx[h] + (objvx[h] + objvxO[h])/2.0L;
objy[h] = objy[h] + (objvy[h] + objvyO[h])/2.0L;
objz[h] = objz[h] + (objvz[h] + objvzO[h])/2.0L;
}

//End of loop computations ------
return 1; //Finished successfully

}```
-Richard

30. Order of Kilopi
Join Date
Sep 2005
Posts
9,031
Now, if you want to see the (PIII using x87 FPU) assembler code the compiler generates for the above little loop, here you go -- Ooops, there's a 15KB limit on posts, I'll split it up:

Code:
```; -- Machine type IA32
; mark_description "Intel(R) C++ Compiler for 32-bit applications, Version 9.1    Build 20061103Z %s";
; mark_description "-Qvc8 -Qlocation,link,S:\\MicrosoftVisualStudio8\\VC\\Bin -c -O3 -Ob2 -Oi -Ot -D _WINDLL -FD -EHs -EHc -MT -";
; mark_description "arch:SSE -YXStdAfx.h -FpIntelReleaseP3/CodeOptomizer.pch -FAs -FaIntelReleaseP3/ -FoIntelReleaseP3/ -W3 -nol";
; mark_description "ogo -Wp64 -Gd -Qrestrict -QxK -Qvc7.1 -Qlocation,link,S:\\MicrosoftVisualStudioNET\\Vc7\\bin";
.686P
.387
OPTION DOTNAME
ASSUME	CS:FLAT,DS:FLAT,SS:FLAT
_TEXT	SEGMENT DWORD PUBLIC FLAT 'CODE'
;	COMDAT _DllMain@12
; -- Begin  _DllMain@12
; mark_begin;
IF @Version GE 612
.MMX
MMWORD TEXTEQU <QWORD>
ENDIF
IF @Version GE 614
.XMM
XMMWORD TEXTEQU <OWORD>
ENDIF
ALIGN     2
PUBLIC _DllMain@12
_DllMain@12	PROC NEAR
; parameter 1: 8 + ebp
; parameter 2: 12 + ebp
; parameter 3: 16 + ebp
\$B1\$1:                          ; Preds \$B1\$0

;;; {

push      ebp                                           ;10.1
mov       ebp, esp                                      ;10.1
sub       esp, 3                                        ;10.1
and       esp, -8                                       ;10.1

;;;     return TRUE;

mov       eax, 1                                        ;11.12
mov       esp, ebp                                      ;11.12
pop       ebp                                           ;11.12
ret       12                                            ;11.12
ALIGN     2
; LOE
; mark_end;
_DllMain@12 ENDP
;_DllMain@12	ENDS
_TEXT	ENDS
_DATA	SEGMENT DWORD PUBLIC FLAT 'DATA'
_DATA	ENDS
; -- End  _DllMain@12
_TEXT	SEGMENT DWORD PUBLIC FLAT 'CODE'
;	COMDAT ?RunLoop@@YGHHPAN0000000PAJ1@Z
; -- Begin  ?RunLoop@@YGHHPAN0000000PAJ1@Z
; mark_begin;
ALIGN     2
PUBLIC ?RunLoop@@YGHHPAN0000000PAJ1@Z
?RunLoop@@YGHHPAN0000000PAJ1@Z	PROC NEAR
; parameter 1: 8 + ebp
; parameter 2: 12 + ebp
; parameter 3: 16 + ebp
; parameter 4: 20 + ebp
; parameter 5: 24 + ebp
; parameter 6: 28 + ebp
; parameter 7: 32 + ebp
; parameter 8: 36 + ebp
; parameter 9: 40 + ebp
; parameter 10: 44 + ebp
; parameter 11: 48 + ebp
\$B2\$1:                          ; Preds \$B2\$0

;;; {

push      ebp                                           ;15.1
mov       ebp, esp                                      ;15.1
and       esp, -8                                       ;15.1
push      edi                                           ;15.1
push      esi                                           ;15.1
push      ebx                                           ;15.1
sub       esp, 180                                      ;15.1
mov       eax, DWORD PTR [ebp+8]                        ;14.36

;;;  if (NumObjects < 1) return 0; //Should be 1 or more

test      eax, eax                                      ;16.2
jle       \$B2\$27        ; Prob 28%                      ;16.2
; LOE
\$B2\$2:                          ; Preds \$B2\$1

;;;
;;;  // static double tM[MSIZE],tM2[MSIZE];
;;;  static double tM,tM2;
;;;  // static double dx[MSIZE], dy[MSIZE], dz[MSIZE];
;;;  static double dx, dy, dz;
;;;  // static double D[MSIZE], f[MSIZE], D2[MSIZE];
;;;  static double D, f1, f2, f3, D2, D3;
;;;  static double objvxO[MSIZE], objvyO[MSIZE], objvzO[MSIZE];
;;;  double fx;
;;;  double fy;
;;;  double fz;
;;;  double qq;
;;;  int Collisions;
;;;  Collisions = 0;
;;;  for (int k=1;k<=NumObjects;k++) {

mov       eax, DWORD PTR [ebp+8]                        ;31.2
test      eax, eax                                      ;31.2
jle       \$B2\$23        ; Prob 1%                       ;31.2
; LOE
\$B2\$3:                          ; Preds \$B2\$2
fldz                                                    ;
fld1                                                    ;
fld       QWORD PTR _2il0floatpacket\$1                  ;

;;;   tM = ObjMass[k] * 398600440000000;

fld       QWORD PTR _2il0floatpacket\$2                  ;32.21
fld       QWORD PTR _2il0floatpacket\$3                  ;
fld       QWORD PTR _2il0floatpacket\$4                  ;
fld       QWORD PTR _2il0floatpacket\$5                  ;
mov       eax, DWORD PTR [ebp+36]                       ;
mov       esi, DWORD PTR [ebp+32]                       ;
mov       ebx, DWORD PTR [ebp+28]                       ;
mov       ecx, DWORD PTR [ebp+12]                       ;
fxch      st(5)                                         ;
fstp      DWORD PTR [esp+136]                           ;
xor       edx, edx                                      ;
mov       edi, 1                                        ;
mov       DWORD PTR [esp+20], edi                       ;
fstp      st(4)                                         ;
mov       DWORD PTR [esp+8], edx                        ;
mov       edx, edi                                      ;
mov       edi, 8                                        ;
lea       eax, DWORD PTR [eax+8]                        ;
lea       esi, DWORD PTR [esi+8]                        ;
mov       DWORD PTR [esp+144], eax                      ;
lea       ebx, DWORD PTR [ebx+8]                        ;
lea       ecx, DWORD PTR [ecx+8]                        ;
fxch      st(3)                                         ;
fstp      st(3)                                         ;
fstp      st(1)                                         ;
; LOE edx ecx ebx esi edi f1 f2 f3
\$B2\$4:                          ; Preds \$B2\$21 \$B2\$3
fstp      st(1)                                         ;
mov       eax, DWORD PTR [esp+144]                      ;
fld       QWORD PTR [ecx]                               ;32.21
fmul      st, st(1)                                     ;32.21
fld       QWORD PTR [ebx]                               ;
fld       QWORD PTR [esi]                               ;
fld       QWORD PTR [eax]                               ;

;;;   objvxO[k] = objvx[k];
;;;   objvyO[k] = objvy[k];
;;;   objvzO[k] = objvz[k];
;;;   for (int j=k+1;j<=NumObjects;j++) {
;;;
;;;     dx = objx[j] - objx[k];
;;;     dy = objy[j] - objy[k];
;;;     dz = objz[j] - objz[k];
;;;     D2 = dx*dx+dy*dy+dz*dz;
;;; 	D3 = (D = sqrt(D2) )* D2;
;;;
;;;     if  (((objsize[k] + objsize[j])/2) > D) {
;;; 	Collisions = Collisions + 1;
;;; 	CollisionObjectA[Collisions] = j;
;;; 	CollisionObjectB[Collisions] = k;
;;;     }
;;;     tM2 = ObjMass[j] * 398600440000000;
;;;     f1 = tM/D3;
;;;     f2 = tM2/D3;
;;;     f3 = f1 + f2;
;;;     qq = 1.0L + f3*(1.0L/6.0L + f3*(1.0L/15.0L + f3/28.0L) );
;;;     if (tM > 0) {

fxch      st(3)                                         ;54.5
fcom      st(5)                                         ;54.5
fnstsw    ax                                            ;54.5
fxch      st(2)                                         ;33.3
fst       QWORD PTR objvxO\$3006\$0\$1[edi]                ;33.3
fxch      st(1)                                         ;34.3
fst       QWORD PTR objvyO\$3006\$0\$1[edi]                ;34.3
fxch      st(3)                                         ;35.3
fst       QWORD PTR objvzO\$3006\$0\$1[edi]                ;35.3
sahf                                                    ;54.5
jbe       \$B2\$13        ; Prob 50%                      ;54.5
; LOE edx ecx ebx esi edi f1 f2 f3 f5 f6 f7
\$B2\$5:                          ; Preds \$B2\$4
fstp      st(0)                                         ;
fstp      st(2)                                         ;
fstp      st(1)                                         ;
mov       eax, DWORD PTR [ebp+8]                        ;
sub       eax, edx                                      ;
mov       DWORD PTR [esp+168], eax                      ;
test      eax, eax                                      ;36.3
jle       \$B2\$21        ; Prob 1%                       ;36.3
; LOE edx ecx ebx esi edi f1 f2 f3
\$B2\$6:                          ; Preds \$B2\$5
mov       eax, DWORD PTR [ebp+16]                       ;
fstp      QWORD PTR [esp+64]                            ;
fxch      st(1)                                         ;
fst       DWORD PTR [esp+140]                           ;
fxch      st(1)                                         ;
mov       DWORD PTR [esp+148], ecx                      ;
mov       DWORD PTR [esp+152], ebx                      ;
lea       eax, DWORD PTR [edi+eax]                      ;
fld       QWORD PTR [eax]                               ;
fst       QWORD PTR [esp+56]                            ;
mov       DWORD PTR [esp+156], eax                      ;
mov       eax, DWORD PTR [ebp+20]                       ;
mov       DWORD PTR [esp], edi                          ;
lea       eax, DWORD PTR [edi+eax]                      ;
fld       QWORD PTR [eax]                               ;
fst       QWORD PTR [esp+48]                            ;
mov       DWORD PTR [esp+164], eax                      ;
mov       eax, DWORD PTR [ebp+24]                       ;
lea       eax, DWORD PTR [edi+eax]                      ;
fld       QWORD PTR [eax]                               ;
fst       QWORD PTR [esp+40]                            ;
mov       DWORD PTR [esp+160], eax                      ;
mov       eax, DWORD PTR [ebp+40]                       ;
lea       eax, DWORD PTR [edi+eax]                      ;
fld       QWORD PTR [eax]                               ;
mov       edi, DWORD PTR [esp+8]                        ;
fst       QWORD PTR [esp+32]                            ;
mov       DWORD PTR [esp+172], eax                      ;
mov       eax, 1                                        ;
mov       DWORD PTR [esp+12], eax                       ;
mov       ebx, eax                                      ;
lea       eax, DWORD PTR [edx+1]                        ;
mov       DWORD PTR [esp+16], eax                       ;
mov       ecx, eax                                      ;
; LOE edx ecx ebx esi edi f1 f2 f4 f5 f6 f7
\$B2\$7:                          ; Preds \$B2\$11 \$B2\$6
fstp      st(0)                                         ;
fstp      st(2)                                         ;
fstp      st(0)                                         ;
fstp      st(0)                                         ;
fstp      st(0)                                         ;
fstp      st(0)                                         ;
mov       eax, DWORD PTR [esp+156]                      ;38.5
fld       QWORD PTR [eax+ebx*8]                         ;38.5
fsub      QWORD PTR [esp+56]                            ;38.20
fld       st(0)                                         ;41.13
mov       eax, DWORD PTR [esp+164]                      ;39.5
fmul      st, st(1)                                     ;41.13
fld       QWORD PTR [eax+ebx*8]                         ;39.5
fsub      QWORD PTR [esp+48]                            ;39.20
fld       st(0)                                         ;41.19
fmul      st, st(1)                                     ;41.19
mov       eax, DWORD PTR [esp+160]                      ;40.5
fld       QWORD PTR [eax+ebx*8]                         ;40.5
fsub      QWORD PTR [esp+40]                            ;40.20
fld       st(0)                                         ;41.25
fmul      st, st(1)                                     ;41.25
mov       eax, DWORD PTR [esp+172]                      ;44.5
fld       QWORD PTR [eax+ebx*8]                         ;44.5
fmul      QWORD PTR _2il0floatpacket\$1                  ;44.37
fxch      st(3)                                         ;41.25
fst       QWORD PTR [esp+88]                            ;41.25
fld       st(0)                                         ;42.12
fsqrt                                                   ;42.12
fxch      st(4)                                         ;44.5
fcomp     st(4)                                         ;44.5
fnstsw    ax                                            ;44.5
fxch      st(3)                                         ;42.12
fst       QWORD PTR [esp+72]                            ;42.12
fmulp     st(3), st                                     ;42.24
sahf                                                    ;44.5
jbe       \$B2\$9         ; Prob 50%                      ;44.5
; LOE edx ecx ebx esi edi f2 f5 f6 f7```

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts
•
here