Lambert93 to WGS84, IGN algorithms, longitude problem

Question

I'm looking for a Java-based solution to convert coordinates from Lambert93 to WGS84.

I found this, "lambert-java", which is an implementation of the algorithms (PDF) published by IGN (French government agency in charge of geographic information).

However, when testing it with data from North of France, the resulting longitude is off by something like 50km. Here's an example:

• Source data (in Lambert93): X = 668832.5384, Y = 6950138.7285
• Geofree conversion to WGS84: lat = 49.64961, lng = 2.56865 (this is correct)
• Output of the above Java program: lat = 49.65262, lng = 1.90503 (this is incorrect)

As far as I can tell, the program is "faithful" in that it faithfully implements the algorithms from IGN; so the fault either lies in the algorithms themselves, or in some parameter or some constant that needs to be adjusted.

How can one sort this out?

(There are other Java-based solutions capable of doing all sorts of GIS manipulations, but since I only need to do this conversion (on lots of data), I'd rather use a very simple, very light package, and "lambert-java" fits this bill nicely; however of course it also needs to be correct...)

Edit 2014-03: this was solved by the original author five months ago and has been working wonderfully ever since.

Answer 1

I've made a lot of simplification. My need is to convert from Lambert93 (IGN) to WGS84. For the region of Bordeaux (south/west of France), the result is exactly the same as Geofree.

import static java.lang.Math.abs;
import static java.lang.Math.atan;
import static java.lang.Math.exp;
import static java.lang.Math.log;
import static java.lang.Math.pow;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;

import java.awt.geom.Point2D;

public class Lambert93 {

   private final static double M_PI_2         = Math.PI / 2.0;
   private final static double DEFAULT_EPS    = 1e-10;
   private final static double E_WGS84        = 0.08181919106;
   private final static double E2             = E_WGS84 / 2.0;
   private final static double LON_MERID_IERS = 3.0 * Math.PI / 180.0;
   private final static double N              = 0.7256077650;
   private final static double C              = 11_754_255.426;
   private final static double XS             =    700_000.000;
   private final static double YS             = 12_655_612.050;

   private static double latitudeFromLatitudeISO( final double latISo ) {
      double phi0 = 2 * atan( exp( latISo ) ) - M_PI_2;
      double phiI = 2
         * atan( pow(( 1 + E_WGS84 * sin( phi0 )) /
                     ( 1 - E_WGS84 * sin( phi0 )), E2) * exp( latISo ) )
         - M_PI_2;
      double delta = abs( phiI - phi0 );
      while( delta > DEFAULT_EPS ){
         phi0 = phiI;
         phiI =
            2 * atan( pow(
                  ( 1 + E_WGS84 * sin( phi0 ) ) /
                  ( 1 - E_WGS84 * sin( phi0 ) )  , E2)
               * exp( latISo ) ) - M_PI_2;
         delta = abs( phiI - phi0 );
      }
      return phiI;
   }

   public static void toLatLon( double x, double y, Point2D.Double out ) {
      final double dX     = x - XS;
      final double dY     = y - YS;
      final double R      = sqrt( dX * dX + dY * dY );
      final double gamma  = atan( dX / -dY );
      final double latIso = -1 / N * log( abs( R / C ) );
      out.x = Math.toDegrees( LON_MERID_IERS + gamma / N );
      out.y = Math.toDegrees( latitudeFromLatitudeISO( latIso ));
   }
}