I'm searching for an algorithm to convert Lambet72 coordinates to lat/long. So for example Lambert72 coordinates 148990,169450 must be converted into 50.83546802746704, 4.354415218851164 but most of the algorithm have a little offset (see attachments).
This is one of the algorithms I found which is close, but still have an error.
Someone has got a better algorithm?
static void Lambert72ToLatLong(double x, double y, ref double longitude, ref double latitude)
{
const double n = 0.77164219;
const double F = 1.81329763;
const double thetaFudge = 0.00014204;
const double e = 0.08199189;
const double a = 6378388;
const double xDiff = 150000;
const double yDiff = 5400088.44;
const 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);
longitude = (theta0 + (theta + thetaFudge) / n) * 180 / Math.PI;
latitude = 0;
for(int i=0; i<10; ++i)
{
latitude = (2 * Math.Atan(Math.Pow(F * a / rho, 1 / n) * Math.Pow((1 + e * Math.Sin(latitude)) / (1 - e * Math.Sin(latitude)), e / 2))) - Math.PI / 2;
}
latitude *= 180 / Math.PI;
}
It seems you are trying to convert Lambert72 (Belgium Dattum 72 LCC 3p) to WGS84 GPS coordinates. Thus, after converting to spherical latitude and longitude coordinates, an additional step has to be performed to end up with the WGS84 coordinates. Here is an alternative code in C# with the full conversion of Lambert72 to WGS74 which complies with the required decimal accurary, as you easily can verify for the coordinates presented in your quoted images.
static void Lambert72toWGS84latlong(double X, double Y)
{
double LongRef = 0.076042943; //
double nLamb = 0.7716421928;
double aCarre = Math.Pow(6378388.0,2.0);
double bLamb = 6378388.0 * (1.0 - (1.0 / 297.0));
double eCarre = (aCarre - Math.Pow(bLamb, 2.0)) / aCarre;
double KLamb = 11565915.812935;
double eLamb = Math.Sqrt(eCarre);
double eSur2 = eLamb / 2.0;
double Tan1 = (X - 150000.012) / (5400088.437 - Y);
double Lambda = LongRef + (1.0 / nLamb) * (0.000142043 + Math.Atan(Tan1));
double RLamb = Math.Sqrt(Math.Pow((X - 150000.012) , 2.0) + Math.Pow ((5400088.437 - Y) ,2.0));
double TanZDemi = Math.Pow((RLamb / KLamb),(1.0 / nLamb));
double Lati1 = 2.0 * Math.Atan(TanZDemi);
double eSin;
double Mult1, Mult2, Mult;
double LatiN, Diff;
double lat, lng ;
int i=0;
do
{
eSin = eLamb * Math.Sin(Lati1);
Mult1 = 1.0 - eSin;
Mult2 = 1.0 + eSin;
Mult = Math.Pow((Mult1 / Mult2) , (eLamb / 2.0));
LatiN = (Math.PI / 2.0) - (2.0 * (Math.Atan(TanZDemi * Mult)));
Diff = LatiN - Lati1;
Lati1 = LatiN;
i++;
} while (Math.Abs(Diff)> 0.0000000277777);
lat=LatiN;
lng=Lambda;
double SinLat = Math.Sin(lat);
double SinLng = Math.Sin(lng);
double CoSinLat = Math.Cos(lat);
double CoSinLng = Math.Cos(lng);
double dx = -125.8;
double dy = 79.9;
double dz = -100.5;
double da = -251.0;
double df = -0.000014192702;
double LWf = 1.0 / 297.0;
double LWa = 6378388.0;
double LWb = (1 - LWf) * LWa;
double LWe2 = (2.0 * LWf) - (LWf * LWf);
double Adb = 1.0 / (1.0 - LWf);
double Rn = LWa / Math.Sqrt(1.0 - LWe2 * SinLat * SinLat);
double Rm = LWa * (1 - LWe2) /Math.Pow((1.0 - LWe2 * lat * lat) ,1.5);
double DLat = -dx * SinLat * CoSinLng - dy * SinLat * SinLng + dz * CoSinLat;
DLat = DLat + da * (Rn * LWe2 * SinLat * CoSinLat) / LWa;
DLat = DLat + df * (Rm * Adb + Rn / Adb) * SinLat * CoSinLat;
DLat = DLat / (Rm + 0.0);
double DLng = (-dx * SinLng + dy * CoSinLng) / ((Rn + 0.0) * CoSinLat);
double Dh = dx * CoSinLat * CoSinLng + dy * CoSinLat * SinLng + dz * SinLat;
Dh = Dh - da * LWa / Rn + df * Rn * lat * lat / Adb;
double LatWGS84 = ((lat + DLat) * 180.0) / Math.PI;
double LngWGS84 = ((lng + DLng) * 180.0) / Math.PI;
MessageBox.Show("WGS84-Latitude=" + LatWGS84.ToString("###.######") +
"--WGS84 Longitude=" + LngWGS84.ToString("###.######"));
}
Hope these are helpful.