No subject

Tue Dec 20 00:28:44 EST 2011

=09lp.phi =3D aacos(hypot(P->thz0 * s, d) * P->rhshz0);
=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

Glynn Clements <glynn at>

More information about the Proj mailing list