[Proj] Problems with /Op option used to compile proj in Microsoft Visual Studio

[Glynn Clements]

Calogero Mauceri wrote:

> > Ultimately, I suspect that the calculation for the two-point
> > equidistant projection just wasn't designed for the case where the two
> > points are around one arc-second apart (approx. 15 metres; would be
> > approx. 40 metres on Earth).
> 
> Glynn, so you are suggesting that more than in the float number 
> approximations, the problem is in why the number passed to the aacos is 
> so close to 1.
> Well, then my question is, which are the cases where the two point 
> equidistant projection could not behave well? Is it a problem of the two 
> point equidistant definition or a problem of the way it is implemented 
> in proj lib?

It's presumably possible to reformulate the calculation to be more
numerically stable, the question is whether it's worthwhile. AFAICT,
the projection seems to be more common at larger scales.

>From the information you provided, I'm assuming that it's failing at
either:

	lp.phi = aacos(hypot(P->thz0 * s, d) * P->rhshz0);
or:
	P->z02 = aacos(P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2));

Looking at how those parameters are derived:

	P->dlam2 = adjlon(lam_2 - lam_1);

	P->cp1 = cos(phi_1);
	P->cp2 = cos(phi_2);
	P->sp1 = sin(phi_1);
	P->sp2 = sin(phi_2);

	P->z02 = aacos(P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2));
	P->hz0 = .5 * P->z02;

	P->thz0 = tan(P->hz0);
	P->rhshz0 = .5 / sin(P->hz0);

	cz1 = cos(hypot(xy.y, xy.x + P->hz0));
	cz2 = cos(hypot(xy.y, xy.x - P->hz0));
	s = cz1 + cz2;
	d = cz1 - cz2;

When lam_1 ~= lam_2 and phi_1 ~= phi_2, then:

	P->dlam2 ~=0
	P->cp1 ~= P->cp2
	P->sp1 ~= P->sp2
	P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2) ~= sin²(phi) + cos²(phi) ~= 1
	cos(P->dlam2) ~= cos(0) ~= 1
	P->z02 ~= aacos(1) ~= 0

So the calculation of P->z02 is likely to be numerically unstable.

Also: cz1 ~= cz2, so d will have poor relative accuracy.

IOW, the formulation used is going to be unstable when the angle
between the reference points is small.

More generally, any arccosine calculation should be seen as a warning
sign regarding numerical stability, as the derivative approaches
infinity as the result approaches zero, meaning that you get the
greatest error magnification at the point where you can least afford
it.

-- 
Glynn Clements <glynn at gclements.plus.com>