[Proj] Computing the Lambert conformal conic projection accurately

[Charles Karney]

I've posted GeographicLib 1.5 on SourceForge.  See

  http://sf.net/projects/geographiclib/files/distrib/

This includes a couple of features of note:

(a) A Matlab interface to the UTMUPS, MGRS, Geoid, Geodesic classes has
    been included.

(b) I've reformulated the Lambert conical projection to improve its
    numerical accuracy.

To follow up on point (b)...  In the 2 parallel case, the Lambert
projection depends on a parameter n, given by (see Snyder, 1987, p. 108)

  n = (log(m1) - log(m2)) / (log(t1) - log(t2))
    = (log(sec(beta2)) - log(sec(beta1))) / 
	     (asinh(tan(chi2)) - asinh(tan(chi1)))

where {phi,beta,chi}{1,2} are the {geographic,reduced,conformal}
latitudes of the standard parallels.  Obviously, this expression is
inaccurate if the standard parallels are close to one another.  The
normal way of dealing with this (proj4, geotrans) is to substitute the
limiting expression

  n = sin(phi1)

if abs(phi2 - phi1) is "sufficiently" small.  However this has the
disadvantage that there's a range of values of (phi2 - phi1) where
neither the limiting form nor the "exact" formula gives accurate results
numerically.

However, there's a systematic method for dealing with numerical problems
like this, namely the calculus of finite differences.  This method
transforms the expression for n to an equivalent one which can be
evaluated accurately for any phi2 != phi1.  The expression is too
lengthy to include here.  However, I've included fairly extensive
comments within the code for the LambertConformalConic class.  The
technique can also be applied to the projection itself and the code then
naturally deals with the limiting cases of the Mercator and polar
stereographic projections.

A useful resource for divided differences is

  W. M. Kahan and R. J. Fateman,
  Symbolic computation of divided differences,
  SIGSAM Bull. 33(3), 7-28 (1999)
  http://doi.acm.org/10.1145/334714.334716
  http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf

(Kahan, incidentally, the "father" of the IEEE floating point standard.)

-- 
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662