Consider "parabola fit" for crossings, also, double-crossings

[lfborjas]

See:

swiss-ephemeris/csrc/swevents.c

Lines 1676 to 1722 in 169e3df

	/* y00, y11, y2 are values for -2*dx, -dx, 0.
	* find zero points of parabola.
	* return: 0 if none
	* 1 if one zero in [-dx.. 0[
	* 1 if both zeros in [-dx.. 0[
	*/
	static int find_zero(double y00, double y11, double y2, double dx,
	double *dxret, double *dxret2)
	{
	double a, b, c, x1, x2;
	c = y11;
	b = (y2 - y00) / 2.0;
	a = (y2 + y00) / 2.0 - c;
	if (b * b - 4 * a * c < 0)
	return 0;
	if (fabs(a) < 1e-100) return 0;
	x1 = (-b + sqrt(b * b - 4 * a * c)) / 2 / a;
	x2 = (-b - sqrt(b * b - 4 * a * c)) / 2 / a;
	// up to here the calcuation was made as if the x-values were -1, 0, 1.
	// This is why below they are shifted by -1
	if (x1 == x2) {
	*dxret = (x1 - 1) * dx;
	*dxret2 = (x1 - 1) * dx;
	return 1;
	}
	if (x1 >=0 && x1 < 1 && x2 >= 0 && x2 < 1) {
	if (x1 > x2) { // two zeroes, order return values
	*dxret = (x2 - 1) * dx;
	*dxret2 = (x1 - 1) * dx;
	} else {
	*dxret = (x1 - 1) * dx;
	*dxret2 = (x2 - 1) * dx;
	}
	return 2;
	}
	if (x1 >=0 && x1 < 1) {
	*dxret = (x1 - 1) * dx;
	*dxret2 = (x2 - 1) * dx; // set this value just in case, should not be used.
	return 1;
	}
	if (x2 >=0 && x2 < 1) {
	*dxret = (x2 - 1) * dx;
	*dxret2 = (x1 - 1) * dx;
	return 1;
	}
	return 0; // should not happen!
	}

A "parabola fit" is often mentioned in the swisseph forum, and I believe it's using the above. Unlike interpolation, it requires bracketing and some cajoling of values into it to account for "circle jumps," and out of it to actually find the equivalent date value; the gain being that it doesn't require as many ephemeris IO events as regular interpolation. This method is also able to detect double crossings, which we currently don't do (though the tools are there! See the comments for crossingBetween)