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 perihelionOutput 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.998654Seems 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.170685Same 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;
}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.
