Add .travis.yml to branch/4.9

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

.travis.yml

@@ -0,0 +1,78 @@
+# This is the config file for building proj.4 and running its test suite
+# with Travis-ci.org
+
+language: c
+
+compiler:
+  - gcc
+  - clang
+
+before_install:
+ - sudo pip install git+git://github.com/tbonfort/cpp-coveralls.git@extensions
+
+install:
+# cmake build
+  - mkdir build_cmake
+  - cd build_cmake
+  - cmake .. -DCMAKE_INSTALL_PREFIX=/tmp/proj_cmake_install
+  - make -j3
+  - make install
+  - find /tmp/proj_cmake_install
+  - cd ..
+# autoconf build
+  - mkdir build_autoconf
+  - cd build_autoconf
+  - ../configure --prefix=/tmp/proj_autoconf_install
+  - make -j3
+  - make install
+  - find /tmp/proj_autoconf_install
+  - make check
+  - cd ..
+# Now with grids
+  - wget http://download.osgeo.org/proj/proj-datumgrid-1.5.zip
+  - cd nad
+  - unzip -o ../proj-datumgrid-1.5.zip
+  - cd ..
+# cmake build with grids
+  - mkdir build_cmake_nad
+  - cd build_cmake_nad
+  - cmake .. -DCMAKE_INSTALL_PREFIX=/tmp/proj_cmake_install_nad
+  - make -j3
+  - make install
+  - find /tmp/proj_cmake_install_nad
+  - cd ..
+# autoconf build with grids
+  - mkdir build_autoconf_nad
+  - cd build_autoconf_nad
+  - ../configure --prefix=/tmp/proj_autoconf_install_nad
+  - make -j3
+  - make install
+  - find /tmp/proj_autoconf_install_nad
+  - make check
+  - cd src
+  - make multistresstest
+  - make test228
+  - cd ..
+  - PROJ_LIB=../nad src/multistresstest
+  - cd ..
+# autoconf build with grids and coverage
+  - CFLAGS="--coverage" LDFLAGS="-lgcov" ./configure
+  - make -j3
+  - make check
+  - mv src/.libs/*.gc* src
+
+script:
+  - echo "done"
+
+after_success:
+  - coveralls --extension .c
+
+notifications:
+  #email:
+  #  recipients:
+  #    - gdal-commits@lists.osgeo.org
+
+  irc:
+    channels:
+      - "irc.freenode.org#gdal"
+    use_notice: true

ChangeLog

@@ -1,3 +1,8 @@
+2014-09-17 Even Rouault <even.rouault@spatialys.com>
+
+	* src/pj_datums.c, src/pj_ellps.c: Add clrk80ign ellipsoid and use it
+	in carthage datum def (#245)
+
 2014-09-13  Frank Warmerdam  <warmerdam@pobox.com>
 
 	* Generate 4.9.0 release.

man/man3/geodesic.3

@@ -1,6 +1,6 @@
 .\" @(#)geodesic.3
 .nr LL 7.0i
-.TH GEODESIC 3 "2013/07/11 Rel. 4.9.0" 
+.TH GEODESIC 3 "2014/12/17 Rel. 4.9.0" 
 .ad b
 .hy 1
 .SH NAME
@@ -38,7 +38,7 @@ measure angles (latitudes, longitudes, and azimuths) in degrees, unlike
 the rest of the \fBproj\fR library, which uses radians.  The
 documentation for this library is included in geodesic.h.  A formatted
 version of the documentation is available at
-http://geographiclib.sf.net/1.32/C
+http://geographiclib.sf.net/1.40/C
 .SH EXAMPLE
 The following program reads in lines with the coordinates for two points
 in decimal degrees (\fIlat1\fR, \fIlon1\fR, \fIlat2\fR, \fIlon2\fR) and
@@ -72,7 +72,7 @@ libproj.a \- library of projections and support procedures
 .SH SEE ALSO
 Full online documentation for \fBgeodesic(3)\fR,
 .br
-http://geographiclib.sf.net/1.32/C
+http://geographiclib.sf.net/1.40/C
 .PP
 .B geod(1)
 .PP

nad/testvarious

@@ -681,6 +681,22 @@ $EXE -f '%.12f' \
  -E >>${OUT} <<EOF
 987122.4183302754536 -14429896.5395309105515
 EOF
+echo "##############################################################" >> ${OUT}
+echo "Test bug 245 (use +datum=carthage)" >> ${OUT}
+#
+$EXE -f '%.12f' \
+     +proj=longlat +datum=WGS84 +nodefs +to +proj=utm +zone=32 +datum=carthage +nodefs \
+ -E >>${OUT} <<EOF
+10 34
+EOF
+echo "##############################################################" >> ${OUT}
+echo "Test bug 245 (use expension of +datum=carthage)" >> ${OUT}
+#
+$EXE -f '%.12f' \
+     +proj=longlat +datum=WGS84 +nodefs +to +proj=utm +zone=32 +a=6378249.2 +b=6356515 +towgs84=-263.0,6.0,431.0 +nodefs \
+ -E >>${OUT} <<EOF
+10 34
+EOF
 ##############################################################################
 # Done!
 # do 'diff' with distribution results

nad/tv_out.dist

@@ -324,3 +324,9 @@ Test bug 244
 ##############################################################
 Test bug 244 (2)
 987122.4183302754536 -14429896.5395309105515	-140.100000000000	-87.000000000000 0.000000000000
+##############################################################
+Test bug 245 (use +datum=carthage)
+10 34	592302.981946153101	3762148.734061029274 -30.311017032713
+##############################################################
+Test bug 245 (use expension of +datum=carthage)
+10 34	592302.981946153101	3762148.734061030205 -30.311017032713

src/geodesic.c

@@ -13,12 +13,12 @@
  *   C. F. F. Karney,
  *   Algorithms for geodesics,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
- *   http://dx.doi.org/10.1007/s00190-012-0578-z
+ *   https://dx.doi.org/10.1007/s00190-012-0578-z
  *   Addenda: http://geographiclib.sf.net/geod-addenda.html
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  */
@@ -222,6 +222,7 @@ static void C1pf(real eps, real c[]);
 static real A2m1f(real eps);
 static void C2f(real eps, real c[]);
 static int transit(real lon1, real lon2);
+static int transitdirect(real lon1, real lon2);
 static void accini(real s[]);
 static void acccopy(const real s[], real t[]);
 static void accadd(real s[], real y);
@@ -271,18 +272,16 @@ void geod_lineinit(struct geod_geodesicline* l,
   l->caps = (caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE) |
     GEOD_LATITUDE | GEOD_AZIMUTH; /* Always allow latitude and azimuth */
 
-  /* Guard against underflow in salp0 */
-  azi1 = AngRound(AngNormalize(azi1));
-  lon1 = AngNormalize(lon1);
   l->lat1 = lat1;
   l->lon1 = lon1;
-  l->azi1 = azi1;
+  /* Guard against underflow in salp0 */
+  l->azi1 = AngRound(AngNormalize(azi1));
   /* alp1 is in [0, pi] */
-  alp1 = azi1 * degree;
+  alp1 = l->azi1 * degree;
   /* Enforce sin(pi) == 0 and cos(pi/2) == 0.  Better to face the ensuing
    * problems directly than to skirt them. */
-  l->salp1 =      azi1  == -180 ? 0 : sin(alp1);
-  l->calp1 = fabs(azi1) ==   90 ? 0 : cos(alp1);
+  l->salp1 =      l->azi1  == -180 ? 0 : sin(alp1);
+  l->calp1 = fabs(l->azi1) ==   90 ? 0 : cos(alp1);
   phi = lat1 * degree;
   /* Ensure cbet1 = +epsilon at poles */
   sbet1 = l->f1 * sin(phi);
@@ -349,7 +348,7 @@ void geod_lineinit(struct geod_geodesicline* l,
 }
 
 real geod_genposition(const struct geod_geodesicline* l,
-                      boolx arcmode, real s12_a12,
+                      unsigned flags, real s12_a12,
                       real* plat2, real* plon2, real* pazi2,
                       real* ps12, real* pm12,
                       real* pM12, real* pM21,
@@ -371,11 +370,11 @@ real geod_genposition(const struct geod_geodesicline* l,
 
   outmask &= l->caps & OUT_ALL;
   if (!( TRUE /*Init()*/ &&
-         (arcmode || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
+         (flags & GEOD_ARCMODE || (l->caps & GEOD_DISTANCE_IN & OUT_ALL)) ))
     /* Uninitialized or impossible distance calculation requested */
     return NaN;
 
-  if (arcmode) {
+  if (flags & GEOD_ARCMODE) {
     real s12a;
     /* Interpret s12_a12 as spherical arc length */
     sig12 = s12_a12 * degree;
@@ -435,7 +434,7 @@ real geod_genposition(const struct geod_geodesicline* l,
   csig2 = l->csig1 * csig12 - l->ssig1 * ssig12;
   dn2 = sqrt(1 + l->k2 * sq(ssig2));
   if (outmask & (GEOD_DISTANCE | GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE)) {
-    if (arcmode || fabs(l->f) > 0.01)
+    if (flags & GEOD_ARCMODE || fabs(l->f) > 0.01)
       B12 = SinCosSeries(TRUE, ssig2, csig2, l->C1a, nC1);
     AB1 = (1 + l->A1m1) * (B12 - l->B11);
   }
@@ -446,26 +445,29 @@ real geod_genposition(const struct geod_geodesicline* l,
   if (cbet2 == 0)
     /* I.e., salp0 = 0, csig2 = 0.  Break the degeneracy in this case */
     cbet2 = csig2 = tiny;
-  /* tan(omg2) = sin(alp0) * tan(sig2) */
-  somg2 = l->salp0 * ssig2; comg2 = csig2;  /* No need to normalize */
   /* tan(alp0) = cos(sig2)*tan(alp2) */
   salp2 = l->salp0; calp2 = l->calp0 * csig2; /* No need to normalize */
-  /* omg12 = omg2 - omg1 */
-  omg12 = atan2(somg2 * l->comg1 - comg2 * l->somg1,
-                comg2 * l->comg1 + somg2 * l->somg1);
 
   if (outmask & GEOD_DISTANCE)
-    s12 = arcmode ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
+    s12 = flags & GEOD_ARCMODE ? l->b * ((1 + l->A1m1) * sig12 + AB1) : s12_a12;
 
   if (outmask & GEOD_LONGITUDE) {
+    /* tan(omg2) = sin(alp0) * tan(sig2) */
+    somg2 = l->salp0 * ssig2; comg2 = csig2;  /* No need to normalize */
+    /* omg12 = omg2 - omg1 */
+    omg12 = flags & GEOD_LONG_NOWRAP ? sig12
+      - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1))
+      + (atan2(somg2, comg2) - atan2(l->somg1, l->comg1))
+      : atan2(somg2 * l->comg1 - comg2 * l->somg1,
+              comg2 * l->comg1 + somg2 * l->somg1);
     lam12 = omg12 + l->A3c *
       ( sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3-1)
                  - l->B31));
     lon12 = lam12 / degree;
     /* Use AngNormalize2 because longitude might have wrapped multiple
      * times. */
-    lon12 = AngNormalize2(lon12);
-    lon2 = AngNormalize(l->lon1 + lon12);
+    lon2 = flags & GEOD_LONG_NOWRAP ? l->lon1 + lon12 :
+      AngNormalize(AngNormalize(l->lon1) + AngNormalize2(lon12));
   }
 
   if (outmask & GEOD_LATITUDE)
@@ -542,7 +544,7 @@ real geod_genposition(const struct geod_geodesicline* l,
   if (outmask & GEOD_AREA)
     *pS12 = S12;
 
-  return arcmode ? s12_a12 : sig12 / degree;
+  return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree;
 }
 
 void geod_position(const struct geod_geodesicline* l, real s12,
@@ -552,7 +554,7 @@ void geod_position(const struct geod_geodesicline* l, real s12,
 
 real geod_gendirect(const struct geod_geodesic* g,
                     real lat1, real lon1, real azi1,
-                    boolx arcmode, real s12_a12,
+                    unsigned flags, real s12_a12,
                     real* plat2, real* plon2, real* pazi2,
                     real* ps12, real* pm12, real* pM12, real* pM21,
                     real* pS12) {
@@ -568,8 +570,9 @@ real geod_gendirect(const struct geod_geodesic* g,
 
   geod_lineinit(&l, g, lat1, lon1, azi1,
                 /* Automatically supply GEOD_DISTANCE_IN if necessary */
-                outmask | (arcmode ? GEOD_NONE : GEOD_DISTANCE_IN));
-  return geod_genposition(&l, arcmode, s12_a12,
+                outmask |
+                (flags & GEOD_ARCMODE ? GEOD_NONE : GEOD_DISTANCE_IN));
+  return geod_genposition(&l, flags, s12_a12,
                           plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
 }
 
@@ -577,7 +580,7 @@ void geod_direct(const struct geod_geodesic* g,
                  real lat1, real lon1, real azi1,
                  real s12,
                  real* plat2, real* plon2, real* pazi2) {
-  geod_gendirect(g, lat1, lon1, azi1, FALSE, s12, plat2, plon2, pazi2,
+  geod_gendirect(g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, plat2, plon2, pazi2,
                  0, 0, 0, 0, 0);
 }
 
@@ -700,7 +703,7 @@ real geod_geninverse(const struct geod_geodesic* g,
     /* Add the check for sig12 since zero length geodesics might yield m12 <
      * 0.  Test case was
      *
-     *    echo 20.001 0 20.001 0 | Geod -i
+     *    echo 20.001 0 20.001 0 | GeodSolve -i
      *
      * In fact, we will have sig12 > pi/2 for meridional geodesic which is
      * not a shortest path. */
@@ -1214,7 +1217,8 @@ real InverseStart(const struct geod_geodesic* g,
       calp1 = sbet12a - cbet2 * sbet1 * sq(somg12) / (1 - comg12);
     }
   }
-  if (salp1 > 0)              /* Sanity check on starting guess */
+  /* Sanity check on starting guess.  Backwards check allows NaN through. */
+  if (!(salp1 <= 0))
     SinCosNorm(&salp1, &calp1);
   else {
     salp1 = 1; calp1 = 0;
@@ -1516,6 +1520,13 @@ int transit(real lon1, real lon2) {
     (lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
 }
 
+int transitdirect(real lon1, real lon2) {
+  lon1 = fmod(lon1, (real)(720));
+  lon2 = fmod(lon2, (real)(720));
+  return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
+           ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
+}
+
 void accini(real s[]) {
   /* Initialize an accumulator; this is an array with two elements. */
   s[0] = s[1] = 0;
@@ -1583,13 +1594,13 @@ void geod_polygon_addedge(const struct geod_geodesic* g,
                           real azi, real s) {
   if (p->num) {                 /* Do nothing is num is zero */
     real lat, lon, S12;
-    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, p->polyline ? 0 : &S12);
     accadd(p->P, s);
     if (!p->polyline) {
       accadd(p->A, S12);
-      p->crossings += transit(p->lon, lon);
+      p->crossings += transitdirect(p->lon, lon);
     }
     p->lat = lat; p->lon = lon;
     ++p->num;
@@ -1720,11 +1731,11 @@ unsigned geod_polygon_testedge(const struct geod_geodesic* g,
   crossings = p->crossings;
   {
     real lat, lon, s12, S12;
-    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+    geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_NOWRAP, s,
                    &lat, &lon, 0,
                    0, 0, 0, 0, &S12);
     tempsum += S12;
-    crossings += transit(p->lon, lon);
+    crossings += transitdirect(p->lon, lon);
     geod_geninverse(g, lat,  lon, p->lat0,  p->lon0,
                     &s12, 0, 0, 0, 0, 0, &S12);
     perimeter += s12;

src/geodesic.h

@@ -4,10 +4,10 @@
  *
  * This an implementation in C of the geodesic algorithms described in
  * - C. F. F. Karney,
- *   <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
+ *   <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
  *   Algorithms for geodesics</a>,
  *   J. Geodesy <b>87</b>, 43--55 (2013);
- *   DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
+ *   DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
  *   10.1007/s00190-012-0578-z</a>;
  *   addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
  *   geod-addenda.html</a>.
@@ -75,18 +75,18 @@
  * (obviously) uniquely defined.  However, in a few special cases there are
  * multiple azimuths which yield the same shortest distance.  Here is a
  * catalog of those cases:
- * - \e lat1 = &minus;\e lat2 (with neither at a pole).  If \e azi1 = \e
- *   azi2, the geodesic is unique.  Otherwise there are two geodesics
+ * - \e lat1 = &minus;\e lat2 (with neither point at a pole).  If \e azi1 =
+ *   \e azi2, the geodesic is unique.  Otherwise there are two geodesics
  *   and the second one is obtained by setting [\e azi1, \e azi2] = [\e
  *   azi2, \e azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 =
  *   &minus;\e S12.  (This occurs when the longitude difference is near
  *   &plusmn;180&deg; for oblate ellipsoids.)
- * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither at a pole).  If
- *   \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
+ * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole).
+ *   If \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
  *   Otherwise there are two geodesics and the second one is obtained by
- *   setting [\e azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e
- *   S12 = &minus;\e S12.  (This occurs when the \e lat2 is near
- *   &minus;\e lat1 for prolate ellipsoids.)
+ *   setting [\e azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e S12
+ *   = &minus;\e S12.  (This occurs when \e lat2 is near &minus;\e lat1 for
+ *   prolate ellipsoids.)
  * - Points 1 and 2 at opposite poles.  There are infinitely many
  *   geodesics which can be generated by setting [\e azi1, \e azi2] =
  *   [\e azi1, \e azi2] + [\e d, &minus;\e d], for arbitrary \e d.  (For
@@ -103,17 +103,17 @@
  * to the member functions.  Most of the internal comments have been retained.
  * However, in the process of transcription some documentation has been lost
  * and the documentation for the C++ classes, GeographicLib::Geodesic,
- * GeographicLib::GeodesicLine, and GeographicLib::PolygonArea, should be
+ * GeographicLib::GeodesicLine, and GeographicLib::PolygonAreaT, should be
  * consulted.  The C++ code remains the "reference implementation".  Think
  * twice about restructuring the internals of the C code since this may make
  * porting fixes from the C++ code more difficult.
  *
- * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.32.
+ * <a href="../index.html">GeographicLib</a> 1.40.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
@@ -128,7 +128,7 @@
  * The minor version of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_MINOR 32
+#define GEODESIC_VERSION_MINOR 40
 /**
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
@@ -381,11 +381,13 @@ extern "C" {
    * @param[in] lat1 latitude of point 1 (degrees).
    * @param[in] lon1 longitude of point 1 (degrees).
    * @param[in] azi1 azimuth at point 1 (degrees).
-   * @param[in] arcmode flag determining the meaning of the \e
-   *   s12_a12.
-   * @param[in] s12_a12 if \e arcmode is 0, this is the distance between
-   *   point 1 and point 2 (meters); otherwise it is the arc length between
-   *   point 1 and point 2 (degrees); it can be negative.
+   * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
+   *   GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
+   *   GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into
+   *   the range [&minus;180&deg;, 180&deg;).
+   * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the distance
+   *   between point 1 and point 2 (meters); otherwise it is the arc length
+   *   between point 1 and point 2 (degrees); it can be negative.
    * @param[out] plat2 pointer to the latitude of point 2 (degrees).
    * @param[out] plon2 pointer to the longitude of point 2 (degrees).
    * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
@@ -402,14 +404,21 @@ extern "C" {
    *
    * \e g must have been initialized with a call to geod_init().  \e lat1
    * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
-   * should be in the range [&minus;540&deg;, 540&deg;).  The function value \e
-   * a12 equals \e s12_a12 is \e arcmode is non-zero.  Any of the "return"
-   * arguments \e plat2, etc., may be replaced by 0, if you do not need some
-   * quantities computed.
+   * should be in the range [&minus;540&deg;, 540&deg;).  The function
+   * value \e a12 equals \e s12_a12 if \e flags & GEOD_ARCMODE.  Any of the
+   * "return" arguments \e plat2, etc., may be replaced by 0, if you do not
+   * need some quantities computed.
+   *
+   * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 &minus;
+   * \e lon1 indicates how many times the geodesic wrapped around the
+   * ellipsoid.  Because \e lon2 might be outside the normal allowed range
+   * for longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it,
+   * e.g., with fmod(\e lon2, 360.0) before using it in subsequent
+   * calculations
    **********************************************************************/
   double geod_gendirect(const struct geod_geodesic* g,
                         double lat1, double lon1, double azi1,
-                        int arcmode, double s12_a12,
+                        unsigned flags, double s12_a12,
                         double* plat2, double* plon2, double* pazi2,
                         double* ps12, double* pm12, double* pM12, double* pM21,
                         double* pS12);
@@ -453,12 +462,16 @@ extern "C" {
    *
    * @param[in] l a pointer to the geod_geodesicline object specifying the
    *   geodesic line.
-   * @param[in] arcmode flag determining the meaning of the second parameter;
-   *   if arcmode is 0, then \e l must have been initialized with \e caps |=
+   * @param[in] flags bitor'ed combination of geod_flags(); \e flags &
+   *   GEOD_ARCMODE determines the meaning of \e s12_a12 and \e flags &
+   *   GEOD_LONG_NOWRAP prevents the value of \e lon2 being wrapped into
+   *   the range [&minus;180&deg;, 180&deg;); if \e flags & GEOD_ARCMODE is
+   *   0, then \e l must have been initialized with \e caps |=
    *   GEOD_DISTANCE_IN.
-   * @param[in] s12_a12 if \e arcmode is 0, this is the distance between
-   *   point 1 and point 2 (meters); otherwise it is the arc length between
-   *   point 1 and point 2 (degrees); it can be negative.
+   * @param[in] s12_a12 if \e flags & GEOD_ARCMODE is 0, this is the
+   *   distance between point 1 and point 2 (meters); otherwise it is the
+   *   arc length between point 1 and point 2 (degrees); it can be
+   *   negative.
    * @param[out] plat2 pointer to the latitude of point 2 (degrees).
    * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
    *   that \e l was initialized with \e caps |= GEOD_LONGITUDE.
@@ -480,11 +493,18 @@ extern "C" {
    * @return \e a12 arc length of between point 1 and point 2 (degrees).
    *
    * \e l must have been initialized with a call to geod_lineinit() with \e
-   * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
-   * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
-   * \e plat2, etc., may be replaced by 0, if you do not need some quantities
-   * computed.  Requesting a value which \e l is not capable of computing is
-   * not an error; the corresponding argument will not be altered.
+   * caps |= GEOD_DISTANCE_IN.  The value \e azi2 returned is in the range
+   * [&minus;180&deg;, 180&deg;).  Any of the "return" arguments \e plat2,
+   * etc., may be replaced by 0, if you do not need some quantities
+   * computed.  Requesting a value which \e l is not capable of computing
+   * is not an error; the corresponding argument will not be altered.
+   *
+   * With \e flags & GEOD_LONG_NOWRAP bit set, the quantity \e lon2 &minus;
+   * \e lon1 indicates how many times the geodesic wrapped around the
+   * ellipsoid.  Because \e lon2 might be outside the normal allowed range
+   * for longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it,
+   * e.g., with fmod(\e lon2, 360.0) before using it in subsequent
+   * calculations
    *
    * Example, compute way points between JFK and Singapore Changi Airport
    * using geod_genposition().  In this example, the points are evenly space in
@@ -494,7 +514,7 @@ extern "C" {
    @code
    struct geod_geodesic g;
    struct geod_geodesicline l;
-   double a12, azi1, lat[101],lon[101];
+   double a12, azi1, lat[101], lon[101];
    int i;
    geod_init(&g, 6378137, 1/298.257223563);
    a12 = geod_geninverse(&g, 40.64, -73.78, 1.36, 103.99,
@@ -508,7 +528,7 @@ extern "C" {
    @endcode
    **********************************************************************/
   double geod_genposition(const struct geod_geodesicline* l,
-                          int arcmode, double s12_a12,
+                          unsigned flags, double s12_a12,
                           double* plat2, double* plon2, double* pazi2,
                           double* ps12, double* pm12,
                           double* pM12, double* pM21,
@@ -526,6 +546,10 @@ extern "C" {
    * polylinep is non-zero, then the vertices and edges define a polyline and
    * only the perimeter is returned by geod_polygon_compute().
    *
+   * The area and perimeter are accumulated at two times the standard floating
+   * point precision to guard against the loss of accuracy with many-sided
+   * polygons.  At any point you can ask for the perimeter and area so far.
+   *
    * An example of the use of this function is given in the documentation for
    * geod_polygon_compute().
    **********************************************************************/
@@ -592,10 +616,12 @@ extern "C" {
    *   polyline (meters).
    * @return the number of points.
    *
-   * Only simple polygons (which are not self-intersecting) are allowed.
-   * There's no need to "close" the polygon by repeating the first vertex.  Set
-   * \e pA or \e pP to zero, if you do not want the corresponding quantity
-   * returned.
+   * The area and perimeter are accumulated at two times the standard floating
+   * point precision to guard against the loss of accuracy with many-sided
+   * polygons.  Only simple polygons (which are not self-intersecting) are
+   * allowed.  There's no need to "close" the polygon by repeating the first
+   * vertex.  Set \e pA or \e pP to zero, if you do not want the corresponding
+   * quantity returned.
    *
    * Example, compute the perimeter and area of the geodesic triangle with
    * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
@@ -704,7 +730,7 @@ extern "C" {
    * area returned is signed with counter-clockwise traversal being treated as
    * positive.
    *
-   * Example, compute the area of Antarctic:
+   * Example, compute the area of Antarctica:
    @code
    double
      lats[] = {-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
@@ -723,7 +749,7 @@ extern "C" {
                         double* pA, double* pP);
 
   /**
-   * mask values for the the \e caps argument to geod_lineinit().
+   * mask values for the \e caps argument to geod_lineinit().
    **********************************************************************/
   enum geod_mask {
     GEOD_NONE         = 0U,                     /**< Calculate nothing */
@@ -738,6 +764,16 @@ extern "C" {
     GEOD_ALL          = 0x7F80U| 0x1FU          /**< Calculate everything */
   };
 
+  /**
+   * flag values for the \e flags argument to geod_gendirect() and
+   * geod_genposition()
+   **********************************************************************/
+  enum geod_flags {
+    GEOD_NOFLAGS      = 0U,     /**< No flags */
+    GEOD_ARCMODE      = 1U<<0,  /**< Position given in terms of arc distance */
+    GEOD_LONG_NOWRAP  = 1U<<15  /**< Don't wrap longitude */
+  };
+
 #if defined(__cplusplus)
 }
 #endif

src/pj_datums.c

@@ -49,7 +49,7 @@ C_NAMESPACE_VAR struct PJ_DATUMS pj_datums[] = {
                                                      "clrk66", 
 				"North_American_Datum_1927",
 "potsdam",  "towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7",  "bessel",  "Potsdam Rauenberg 1950 DHDN",
-"carthage",  "towgs84=-263.0,6.0,431.0",  "clrk80",  "Carthage 1934 Tunisia",
+"carthage",  "towgs84=-263.0,6.0,431.0",  "clrk80ign",  "Carthage 1934 Tunisia",
 "hermannskogel", "towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232",  "bessel",  "Hermannskogel",
 "ire65",  "towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15",  "mod_airy",  "Ireland 1965",
 "nzgd49",    "towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", "intl", "New Zealand Geodetic Datum 1949",

src/pj_ellps.c

@@ -19,6 +19,7 @@ pj_ellps[] = {
 "bess_nam",	"a=6377483.865", "rf=299.1528128", "Bessel 1841 (Namibia)",
 "clrk66",	"a=6378206.4", "b=6356583.8", "Clarke 1866",
 "clrk80",	"a=6378249.145", "rf=293.4663", "Clarke 1880 mod.",
+"clrk80ign","a=6378249.2", "rf=293.4660212936269", "Clarke 1880 (IGN).",
 "CPM",  	"a=6375738.7", "rf=334.29", "Comm. des Poids et Mesures 1799",
 "delmbr",	"a=6376428.",  "rf=311.5", "Delambre 1810 (Belgium)",
 "engelis",	"a=6378136.05", "rf=298.2566", "Engelis 1985",

src/pj_transform.c

@@ -212,7 +212,7 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
 /*      Do we need to translate from geoid to ellipsoidal vertical      */
 /*      datum?                                                          */
 /* -------------------------------------------------------------------- */
-    if( srcdefn->has_geoid_vgrids )
+    if( srcdefn->has_geoid_vgrids && z != NULL )
     {
         if( pj_apply_vgridshift( srcdefn, "sgeoidgrids", 
                                  &(srcdefn->vgridlist_geoid), 
@@ -237,7 +237,7 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
 /*      Do we need to translate from geoid to ellipsoidal vertical      */
 /*      datum?                                                          */
 /* -------------------------------------------------------------------- */
-    if( dstdefn->has_geoid_vgrids )
+    if( dstdefn->has_geoid_vgrids && z != NULL )
     {
         if( pj_apply_vgridshift( dstdefn, "sgeoidgrids", 
                                  &(dstdefn->vgridlist_geoid),