Zamiana oraz transformacja współrzędnych GPS na współrzędne płaskie [csharp]

Bardzo długo szukałem podpowiedzi na nurtujące mnie pytanie odnośnie zamiany współrzędnych GPS na współrzędne płaskie by móc obliczyć np długość od jednego pkt do drugiego lub też przeliczyć pole powierzchni obszaru. No i znalazłem na stronie AGH imienia Stanisława Staszica w Krakowie dość proste rozwiązanie. Po niedługim czasie linki wygasły więc postanowiłem je zachować. Link do skryptu w csharp znajduje się na wkej.org.

using System;

using System.Collections.Generic;
using System.Text;


namespace FreeMierniczy
{ 
    public class GeodGPS
    {
        public const double U92_L_ZERO = 19.0;    //Poludnik srodkowy
        public const double U92_FE = 500000.0;    // False Eting
        public const double U92_FN = -5300000.0;  // False Northing
        public const double U92_M_ZERO = 0.9993;  // Wsp. skali
        public const double U92_E = 0.0818191910428;
        // pozostale stale
        public const double PI = 3.14159265358979;
        public const double RO = 6367449.14577;
        public const double A2 = 0.0008377318247344;
        public const double A4 = 0.0000007608527788826;
        public const double A6 = 0.000000001197638019173;
        public const double A8 = 0.00000000000244337624251;
        /* Obliczenie współrzędnej X na podstawie współrzędnych geodezyjnych B,L
        Brad = B * Pi / 180
        Lrad = L * Pi / 180
        Lorad = Lo * Pi / 180
        k = ((1 - E * Sin(Brad)) / (1 + E * Sin(Brad))) ^ (E / 2)
        C = k * Tan((Brad / 2) + (Pi / 4))
        firad = (2 * Atn(C)) - (Pi / 2)
        Xmer = Atn(Sin(firad) / (Cos(firad) * Cos(Lrad - Lorad)))
        Ymer = 0.5 * Log((1 + Cos(firad) * Sin(Lrad - Lorad)) / (1 - Cos(firad) * Sin(Lrad - Lorad)))
        Xgk = Ro * (Xmer + (a2 * Sin(2 * Xmer) *Cosh(2 * Ymer)) + (a4 * Sin(4 * Xmer) * Cosh(4 * Ymer)) + (a6 * Sin(6 * Xmer) *Cosh(6 * Ymer)) + (a8 * Sin(8 * Xmer) *Cosh(8 * Ymer)))
        Współrzędna X Wgs84 = mo * Xgk + FN */

        public double getXU92(double longitude, double latitude)
        {
            
            double brad, lrad, lorad, k, c, firad, xmer, ymer, xgk;
            brad = latitude * PI / 180;
            lrad = longitude * PI / 180;
            lorad = U92_L_ZERO * PI / 180;
            k = Math.Pow(((1 - U92_E * Math.Sin(brad)) / (1 + U92_E * Math.Sin(brad))), U92_E / 2);
            c = k * Math.Tan(brad / 2 + PI / 4);
            firad = 2 * Math.Atan(c) - PI / 2;
            xmer = Math.Atan(Math.Sin(firad) / (Math.Cos(firad) * Math.Cos(lrad - lorad)));
            ymer = 0.5 * Math.Log((1 + Math.Cos(firad) * Math.Sin(lrad - lorad)) / (1 - Math.Cos(firad) * Math.Sin(lrad - lorad)));
            xgk = RO * (xmer + A2 * Math.Sin(2 * xmer) * Math.Cosh(2 * ymer) + 
                A4 * Math.Sin(4 * xmer) * Math.Cosh(4 * ymer) + 
                A6 * Math.Sin(6 * xmer) * Math.Cosh(6 * ymer) + 
                A8 * Math.Sin(8 * xmer) * Math.Cosh(8 * ymer));
            return U92_M_ZERO * xgk + U92_FN;

        }
        /* Obliczenie współrzędnej Y na podstawie współrzędnych geodezyjnych B,L
         Brad = B * Pi / 180
         Lrad = L * Pi / 180
         Lorad = Lo * Pi / 180
         k = ((1 - E * Sin(Brad)) / (1 + E * Sin(Brad))) ^ (E / 2)
         C = k * Tan((Brad / 2) + (Pi / 4))
         firad = (2 * Atn(C)) - (Pi / 2)
         Xmer = Atn(Sin(firad) / (Cos(firad) * Cos(Lrad - Lorad)))
         Ymer = 0.5 * Log((1 + Cos(firad) * Sin(Lrad - Lorad)) / (1 - Cos(firad) * Sin(Lrad - Lorad)))
         Ygk = Ro * (Ymer + (a2 * Cos(2 * Xmer) * Sinh(2 * Ymer)) + (a4 * Cos(4 * Xmer) * Sinh(4 * Ymer)) + (a6 * Cos(6 * Xmer) *Sinh(6 * Ymer)) + (a8 * Cos(8 * Xmer) * Sinh(8 * Ymer)))
         Współrzędna Y Wgs84 = mo * Ygk + FE */
        public double getYU92(double longitude, double latitude)
        {
            double brad, lrad, lorad, k, c, firad, xmer, ymer, ygk;
            brad = latitude * PI / 180;
            lrad = longitude * PI / 180;
            lorad = U92_L_ZERO * PI / 180;
            k = Math.Pow(((1 - U92_E * Math.Sin(brad)) / (1 + U92_E * Math.Sin(brad))), U92_E / 2);
            c = k * Math.Tan(brad / 2 + PI / 4);
            firad = 2 * Math.Atan(c) - PI / 2;
            xmer = Math.Atan(Math.Sin(firad) / (Math.Cos(firad) * Math.Cos(lrad - lorad)));
            ymer = 0.5 * Math.Log((1 + Math.Cos(firad) * Math.Sin(lrad - lorad)) / (1 - Math.Cos(firad) * Math.Sin(lrad - lorad)));
            ygk = RO * (ymer + A2 * Math.Cos(2 * xmer) * Math.Sinh(2 * ymer) +
                A4 * Math.Cos(4 * xmer) * Math.Sinh(4 * ymer) +
                A6 * Math.Cos(6 * xmer) * Math.Sinh(6 * ymer) +
                A8 * Math.Cos(8 * xmer) * Math.Sinh(8 * ymer));
            return U92_M_ZERO * ygk + U92_FE;
        }

       
    }
}