taby said:
Angle is the angle between two subsequent aphelions.

Wtf is an aphelion? Some alien species, maybe green unicorns with tentacles i would guess?

Please speak in geometrical terms. :D

Yep, your dot is always extremely close to one when you measure, so you get noise and nans all the time.

Seems no nans happen. I've tried it with enabled fp exceptions, but did not trigger anything. :(
Which is strange, since frob said he got exceptions iirc.

Here is the code i use to enable it, stolen from the physics engine:

#include "SystemTools.h"
#include <stdio.h>
#include <direct.h>
#include <float.h>
#include <pmmintrin.h>
#include <mutex>
#include <algorithm>
#include <windows.h>
#include <Mmsystem.h>

void SystemTools::EnableFPExceptions ()
{
	#if defined (_MSC_VER)
		uint32_t mask = (_EM_ZERODIVIDE | _EM_INVALID | _EM_DENORMAL);
	#else
		uint32_t mask = D_FLOAT_EXCEPTIONS_MASK	(FE_DIVBYZERO | FE_INVALID | FE_INEXACT);
	#endif

	#if defined (_MSC_VER)
		_clearfp();
		uint32_t m_floatMask = _controlfp(0, 0);
		_controlfp(m_floatMask & ~mask, _MCW_EM);
	#endif
	
	#if (defined(_M_IX86) || defined(__x86_64__) || defined(_M_X64))
		uint32_t m_simdMask = _mm_getcsr();
		_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
		_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
	#endif
}

How would I use atan2 to emulate acos? Sorry, I'm still a trigonometric newbie.

The aphelion is the furthest point away from the Sun in the orbit. The perihelion is the closest point.

This works better. Gives 34 arcseconds per Earth century, where I want 43 arcseconds per Earth century.

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 = numeric_limits<float>::epsilon();// powf(2.0f, -23.0f);
	const float remainder = fmodf(d, epsilon);

	d += remainder;

	return d;
}

Do you see now that as we inch toward a snap to float solution, the precession inches toward being 43 arcseconds per Earth century?

This one gives an answer of about 47 arcseconds per Earth century.

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 = numeric_limits<float>::epsilon();// powf(2.0f, -23.0f);
	const double remainder = fmod(d, epsilon);

	d = nexttowardf(d, d - remainder);

	return d;
}

Where:

return static_cast<double>(static_cast<float>(d));

the answer is 47 arc seconds per century as well.

The answers are practically identical.

taby said:
How would I use atan2 to emulate acos?

It requires two orthonormal basis vectors to define the coordinate system of measurements.
For an arbitrary vector: float angle = atan2(basis1.Dot(vector), basis2.Dot(vector)).
The vector does not need to be unit length, eliminating some noise. And we get a signed angle result.
If the coords system is world space, it's simply angle = atan2(vector.x, vector.y).

taby said:
The aphelion is the furthest point away from the Sun in the orbit. The perihelion is the closest point.

So you you measure from last aphelion to current aphelion? Or from aphelion to perihelion?

The former would be just 0 or 360 degrees all the time, since the point is always at the same spot.
Only the latter would give a robust measurement.

It would make sense to measure time between last aphelion and current aphelion,
but it does not make sense to measure angle.

taby said:
Do you see now that as we inch toward a snap to float solution, the precession inches toward being 43 arcseconds per Earth century?

All i see is that nature does not snap or quantize at astronomic scales.
But arcseconds implies it's about a duration. Could you use time to confirm your results?

Thanks again for your diligence!!!

JoeJ said:

It requires two orthonormal basis vectors to define the coordinate system of measurements.
For an arbitrary vector: float angle = atan2(basis1.Dot(vector), basis2.Dot(vector)).
The vector does not need to be unit length, eliminating some noise. And we get a signed angle result.
If the coords system is world space, it's simply angle = atan2(vector.x, vector.y).

OK, that's great to know.

So you you measure from last aphelion to current aphelion? Or from aphelion to perihelion?

I measure from last aphelion to current aphelion, yes.

The former would be just 0 or 360 degrees all the time, since the point is always at the same spot.
Only the latter would give a robust measurement.

Not really. The only time you get zero arcseconds per Earth century is when alpha = beta = 1 and dt is small enough, like 0.01 when using the Euler integration.

It would make sense to measure time between last aphelion and current aphelion,
but it does not make sense to measure angle.

I think there's a misunderstanding here, and it's my fault: arcseconds is a measure of angle, not time. It's analogous to radians.

taby said:
Do you see now that as we inch toward a snap to float solution, the precession inches toward being 43 arcseconds per Earth century?

All i see is that nature does not snap or quantize at astronomic scales.
But arcseconds implies it's about a duration. Could you use time to confirm your results?

I suppose you could note the duration, if you like. The duration is like 88 Earth days.

P.S. Here is an image… delta_p is the angle from one aphelion to the next. It's from the paper https://www.techrxiv.org/doi/full/10.36227/techrxiv.171470734.46233134/v1​

The effect here is grossly exaggerated.

P.S. The snapping-to occurs on the microscopic scale (which we have set to 0.01 when using Euler integration). The snapping-to occurs at every step. It's a very tiny amount, yes, but at the macroscopic scale it all turns out to be pretty much exactly what is observed (47 versus 43 arcseconds per Earth century).

Note that 43 arcseconds is like 0.0002 radian… a macrosliver… made up of even tinier microslivers (that is, 88 Earth days / 0.01 seconds = 760320000 microslivers in total). For all of the values possible, it just happens to luck out? No, it's a systemic property, I used google to do these arithmetic. :)

Note that the symplectic and Euler integration give roughly the same amount. Using Euler is not a weakness – it serves as a reference.