taby said:
This doesn't work. I checked the value of max, and it is 2^64 - 1, as expected.

Indeed. Likely i had tomatoes on my eyes. Running your code again shows the same results as mine.

So it works. Why do you say it does not? Different results in the simulation after how many integration steps?

Maybe you have denormal / infinite / nan numbers sometimes?

I would check this with std::isfinite() to be sure. That's the kind of stuff where enabling floating point exceptions is helpful, as you get a debugbreak without any checks.

Running your code millions of times…

Output with the cast:

[120079438] hit perihelion
[0] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 45.734582
 angle 45.734582
 delta 42.936958
 avg   45.734582
[360238317] hit perihelion
[1] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 90.843000
 angle 45.108419
 delta 42.936958
 avg   45.421500
[600397195] hit perihelion
[2] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 140.822554
 angle 49.979553
 delta 42.936958
 avg   46.940851
[840556073] hit perihelion
[3] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 186.503801
 angle 45.681247
 delta 42.936958
 avg   46.625950
[1080714949] hit perihelion
[4] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 231.594201
 angle 45.090401
 delta 42.936958
 avg   46.318840
[1320873825] hit perihelion
[5] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 281.541227
 angle 49.947025
 delta 42.936958
 avg   46.923538
[1561032701] hit perihelion
[6] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 330.816692
 angle 49.275466
 delta 42.936958
 avg   47.259527
[1801191575] hit perihelion
[7] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 375.780764
 angle 44.964071
 delta 42.936958
 avg   46.972595
[2041350450] hit perihelion

Output with masked bits:

[120079431] hit perihelion
[0] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 3.605790
 angle 3.605790
 delta 42.936958
 avg   3.605790
[360238296] hit perihelion
[1] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 12.062117
 angle 8.456327
 delta 42.936958
 avg   6.031059
[600397161] hit perihelion
[2] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 16.093513
 angle 4.031396
 delta 42.936958
 avg   5.364504
[840556023] hit perihelion
[3] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 25.017388
 angle 8.923875
 delta 42.936958
 avg   6.254347
[1080714886] hit perihelion
[4] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 31.391584
 angle 6.374196
 delta 42.936958
 avg   6.278317
[1320873748] hit perihelion
[5] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 39.250222
 angle 7.858637
 delta 42.936958
 avg   6.541704
[1561032609] hit perihelion
[6] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 48.443218
 angle 9.192997
 delta 42.936958
 avg   6.920460
[1801191471] hit perihelion
[7] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 51.565924
 angle 3.122706
 delta 42.936958
 avg   6.445741
[2041350330] hit perihelion
[8] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 59.320469
 angle 7.754545
 delta 42.936958
 avg   6.591163
[-2013458106] hit perihelion
[9] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 69.986539
 angle 10.666071
 delta 42.936958
 avg   6.998654

Seems the same!

Running without without messing on beta at all:

[120079432] hit perihelion
[0] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 12.684491
 angle 12.684491
 delta 42.936958
 avg   12.684491
[360238298] hit perihelion
[1] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 26.591345
 angle 13.906855
 delta 42.936958
 avg   13.295673
[600397162] hit perihelion
[2] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 39.901058
 angle 13.309713
 delta 42.936958
 avg   13.300353
[840556027] hit perihelion
[3] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 54.380449
 angle 14.479391
 delta 42.936958
 avg   13.595112
[1080714890] hit perihelion
[4] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 70.045937
 angle 15.665487
 delta 42.936958
 avg   14.009187
[1320873753] hit perihelion
[5] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 85.183814
 angle 15.137878
 delta 42.936958
 avg   14.197302
[1561032615] hit perihelion
[6] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 97.804080
 angle 12.620265
 delta 42.936958
 avg   13.972011
[1801191477] hit perihelion
[7] hit aphelion
 orbit 0.000000
 dot   1.000000
 total 113.365476
 angle 15.561397
 delta 42.936958
 avg   14.170685

Same results as well.
Comparing numbers for delta on [7], all methods give the same result of 42.936958, 42.936958, 42.936958
Is this any surprise to you? What do i miss?

(I have accidentally swapped framecount ([ ]) and orbitcount for my output, so those printed numbers are garbage, but gave a nice index to compare)

Angle is the variable in question. The delta is constant.

Only the cast answer gives the proper result for angle. Why? I’m not sure.

for fun, try setting alpha and beta to 1. That should give you a precession of zero degrees, especially as the time step is lowered.

Ok here’s another way to look at it…

is there a way to calculate the number of mantissa bits needed for a certain epsilon like numericlimits<float>::epsilon()? I mean, of course there’s a solution but I am not sure what it is exactly. It involves log that’s For sure.

Actually, I found the right math at: https://en.wikipedia.org/wiki/Machine_epsilon​

So, for our purposes, epsilon == 2^(-mantissa_bits) where mantissa_bits is 23 for floats.

taby said:
Angle is the variable in question. The delta is constant.

Ah ok. But what does angle mean?
Why is it measured and accumulated to total only when the planet comes closer to the sun?
The whole math here seems not robust. Errors will multiply and accumulate, and i still think the ‘correct’ result is random luck. There should be a better way to measure this, but i do not yet see what it tries to measure at all.

I added code to render the trajectory, and have drawn one orbit for each of the 3 methods in question over each other. It looks the same.

This code gives an answer of about 20. That's closer than the 43 given by the analytical solution, but obviously not quite right.

This is what I'm trying:

double truncate_normalized_double(double d)
{
	if (d <= 0.0)
		return 0.0;
	else if (d >= 1.0)
		return 1.0;

	static const double epsilon = pow(2.0, -23.0); 
	double remainder = fmod(d, epsilon);

	d += remainder;

	return d;
}

@JoeJ Angle is the angle between two subsequent aphelions.

I see you calculate angle from a dot product of unit vectors.
But you do not clamp the dot to be within (-1,1), so you'll get nan from acos half of the time.
The other half will be just noise, since the vectors are so close together.

But i may imagine wrong when the previous direction is updated. I'm assuming you do it each step, which likely is wrong.
But in general the clamping is always necessary with acos. So for high accuracy i would use atan2.