Convert LAMBERT 1972 to WGS84

Needed to convert LAMBERT 1972 (EPSG:31370) coordinates to WGS84 (EPSG:4326) for a project I’m working on in order to use the coordinates with Google Maps (which uses a SphericalMercator (EPSG:900913) projection).

Took me quite a while to find the needed algorithm, so hereby I’m listing them (and am providing the JavaScript implementation):

JavaScript:

var lambert72toWGS84 = function(x, y){

	var newLongitude, newLatitude;

	var n = 0.77164219,
	    F = 1.81329763,
		thetaFudge = 0.00014204,
		e = 0.08199189,
		a = 6378388,
		xDiff = 149910,
		yDiff = 5400150,
		theta0 = 0.07604294;

	var xReal = xDiff - x,
	    yReal = yDiff - y;

	var rho = Math.sqrt(xReal * xReal + yReal * yReal),
	    theta = Math.atan(xReal / -yReal);

	newLongitude = (theta0 + (theta + thetaFudge) / n) * 180 / Math.PI;
	newLatitude = 0;

	for (var i = 0; i < 5 ; ++i) {
		newLatitude = (2 * Math.atan(Math.pow(F * a / rho, 1 / n) * Math.pow((1 + e * Math.sin(newLatitude)) / (1 - e * Math.sin(newLatitude)), e / 2))) - Math.PI / 2;
	}
	newLatitude *= 180 / Math.PI;
	return [newLatitude, newLongitude];

}

Java:

public static Point lambert72toWGS84(double x, double y) {

	double newLongitude;
	double newLatitude;

	double n = 0.77164219;
	double F = 1.81329763;
	double thetaFudge = 0.00014204;
	double e = 0.08199189;
	double a = 6378388;
	double xDiff = 149910;
	double yDiff = 5400150;

	double theta0 = 0.07604294;

	double xReal = xDiff - x;
	double yReal = yDiff - y;

	double rho = Math.sqrt(xReal * xReal + yReal * yReal);
	double theta = Math.atan(xReal / -yReal);

	newLongitude = (theta0 + (theta + thetaFudge) / n) * 180 / Math.PI;
	newLatitude = 0;

	for (int i = 0; i < 5 ; ++i) {
		newLatitude = (2 * Math.atan(Math.pow(F * a / rho, 1 / n) * Math.pow((1 + e * Math.sin(newLatitude)) / (1 - e * Math.sin(newLatitude)), e / 2))) - Math.PI / 2;
	}
	newLatitude *= 180 / Math.PI;
	return new Point(newLatitude, newLongitude);
}

C/C++:

void LambertToLongLat(double x,double y, double& longt, double& lat){
    double n = 0.77164219;
    double F = 1.81329763;
    double thetaFudge = 0.00014204;
    double e= 0.08199189;
    double a= 6378388;
    double xDiff= 150000;
    double yDiff = 5400088.44;
    double theta0 = 0.07604294;

    double xReal = xDiff-x;
    double yReal = yDiff-y;

    double rho = sqrt(xReal*xReal + yReal * yReal);
    double theta = atan(xReal/-yReal);

    longt = (theta0 + (theta+thetaFudge) /n)*180/PI;
    lat =0;
    for(size_t i=0; i<5; ++i){
        lat =(2*atan(pow(F*a/rho,1/n)* pow((1+ e*sin(lat))/(1-e*sin(lat)),e/2)))-PI/2;
    }
    lat *=180/PI;
}

Hope this helps!