replacement of geodesic engine with one from Charles Karney (#197)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2297 4e78687f-474d-0410-85f9-8d5e500ac6b2

man/man1/geod.1

@@ -2,7 +2,7 @@
 .nr LL 5.5i
 .ad b
 .hy 1
-.TH GEOD 1 "2000/03/21 Rel. 4.4" 
+.TH GEOD 1 "2012/11/29 Rel. 4.8"
 .SH NAME
 geod \- direct geodesic computations
 .br
@@ -45,7 +45,7 @@ The following runline control parameters can appear in any order:
 .B \-I
 Specifies that the inverse geodesic computation is to be performed.
 May be used with execution of
-.B goed
+.B geod
 as an alternative to
 .B invgeod
 execution.
@@ -179,7 +179,7 @@ U.S. stature miles from Boston, MA, to Portland, OR:
 .RE
 which gives the results:
 .RS 5
- \f(CW-66d31'50.141"   75d39'13.083"   2587.504
+ \f(CW-66d31'50.141" 75d39'13.083" 2587.504
 .RE
 where the first two values are the
 azimuth from Boston to Portland, the back azimuth from Portland to
@@ -194,13 +194,14 @@ Portland's location by azimuth and distance:
 .RE
 which gives:
 .RS 5
- \f(CW45d31'0.003"N   123d40'59.985"W 75d39'13.094"\fR
+ \f(CW45d31'0.003"N 123d40'59.985"W 75d39'13.094"\fR
 .RE
 Note: lack of precision in the distance value compromises
 the precision of the Portland location.
 .SH SEE ALSO
-Thomas, P.D., 1970,
-.I "Spheroidal Geodesics, Reference Systems & Local Geometry:"
-U.S. Naval Oceanographic  Office, S-138.
+C. F. F. Karney,
+.I "Algorithms for Geodesics," J. Geodesy (2012);
+DOI: 10.1007/s00190-012-0578-z,
+http://geographiclib.sf.net/geod-addenda.html.
 .SH HOME PAGE
 http://www.remotesensing.org/proj

src/Makefile.am

@@ -11,7 +11,8 @@ EXTRA_DIST = makefile.vc proj.def
 proj_SOURCES = proj.c gen_cheb.c p_series.c
 cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
 nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_for.c geod_inv.c geodesic.h
+geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h \
+	geodesic.c geodesic.h
 multistresstest_SOURCES = multistresstest.c
 
 proj_LDADD = libproj.la

src/geod.c

@@ -1,6 +1,6 @@
 /* <<<< Geodesic filter program >>>> */
 # include "projects.h"
-# include "geodesic.h"
+# include "geod_interface.h"
 # include "emess.h"
 # include <ctype.h>
 # include <stdio.h>

src/geod_interface.c

@@ -0,0 +1,34 @@
+#include "projects.h"
+#include "geod_interface.h"
+
+void geod_ini(void) {
+  GeodesicInit(&GlobalGeodesic, geod_a, geod_f);
+}
+
+void geod_pre(void) {
+  double
+    degree = PI/180,
+    lat1 = phi1 / degree, lon1 = lam1 /degree, azi1 = al12 / degree;
+  GeodesicLineInit(&GlobalGeodesicLine, &GlobalGeodesic,
+		   lat1, lon1, azi1, 0U);
+}
+
+void geod_for(void) {
+  double degree = PI/180, s12 = geod_S, lat2, lon2, azi2;
+  Position(&GlobalGeodesicLine, s12, &lat2, &lon2, &azi2);
+  azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
+  phi2 = lat2 * degree;
+  lam2 = lon2 * degree;
+  al21 = azi2 * degree;
+}
+
+void geod_inv(void) {
+  double
+    degree = PI / 180,
+    lat1 = phi1 / degree, lon1 = lam1 / degree,
+    lat2 = phi2 / degree, lon2 = lam2 / degree,
+    azi1, azi2, s12;
+  Inverse(&GlobalGeodesic, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2);
+  azi2 += azi2 >= 0 ? -180 : 180; /* Compute back azimuth */
+  al12 = azi1 * degree; al21 = azi2 * degree; geod_S = s12;
+}

src/geod_interface.h

@@ -0,0 +1,45 @@
+#if !defined(GEOD_INTERFACE_H)
+#define GEOD_INTERFACE_H
+
+#include "geodesic.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#ifndef _IN_GEOD_SET
+#  define GEOD_EXTERN extern
+#else
+#  define GEOD_EXTERN
+#endif
+
+GEOD_EXTERN struct geodesic {
+  double A, FLAT, LAM1, PHI1, ALPHA12, LAM2, PHI2, ALPHA21, DIST;
+} GEODESIC;
+
+# define geod_a	GEODESIC.A
+# define geod_f	GEODESIC.FLAT
+# define lam1	GEODESIC.LAM1
+# define phi1	GEODESIC.PHI1
+# define al12	GEODESIC.ALPHA12
+# define lam2	GEODESIC.LAM2
+# define phi2	GEODESIC.PHI2
+# define al21	GEODESIC.ALPHA21
+# define geod_S	GEODESIC.DIST
+
+GEOD_EXTERN struct Geodesic GlobalGeodesic;
+GEOD_EXTERN struct GeodesicLine GlobalGeodesicLine;
+GEOD_EXTERN int n_alpha, n_S;
+GEOD_EXTERN double to_meter, fr_meter, del_alpha;
+
+void geod_set(int, char **);
+void geod_ini(void);
+void geod_pre(void);
+void geod_for(void);
+void geod_inv(void);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

src/geod_set.c

@@ -3,7 +3,7 @@
 
 #include <string.h>
 #include "projects.h"
-#include "geodesic.h"
+#include "geod_interface.h"
 #include "emess.h"
 	void
 geod_set(int argc, char **argv) {
@@ -32,16 +32,8 @@ geod_set(int argc, char **argv) {
 		fr_meter = 1. / (to_meter = atof(unit_list[i].to_meter));
 	} else
 		to_meter = fr_meter = 1.;
-	if ((ellipse = es) != 0.) {
-		onef = sqrt(1. - es);
-		geod_f = 1 - onef;
-		f2 = geod_f/2;
-		f4 = geod_f/4;
-		f64 = geod_f*geod_f/64;
-	} else {
-		onef = 1.;
-		geod_f = f2 = f4 = f64 = 0.;
-	}
+	geod_f = es/(1 + sqrt(1 - es));
+	geod_ini();
 	/* check if line or arc mode */
 	if (pj_param(NULL,start, "tlat_1").i) {
 		double del_S;