No subject

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


Looking at how those parameters are derived:

=09P->dlam2 =3D adjlon(lam_2 - lam_1);

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

=09P->z02 =3D aacos(P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2));=

=09P->hz0 =3D .5 * P->z02;

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

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

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

=09P->dlam2 ~=3D0
=09P->cp1 ~=3D P->cp2
=09P->sp1 ~=3D P->sp2
=09P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2) ~=3D sin=B2(phi) +=
 cos=B2(phi) ~=3D 1
=09cos(P->dlam2) ~=3D cos(0) ~=3D 1
=09P->z02 ~=3D aacos(1) ~=3D 0

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

Also: cz1 ~=3D 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.

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