Geodetic location helper, for translating RT90 coordinates from/To WGS84 coordinates.

GeoLocationHelper.cs

public static class GeoLocationHelper
{
    public static GeoCoordinate ToWGS84(double x, double y)
    {
        var axis = 6378137.0;
        var flattening = 1.0 / 298.257222101;
        var centralMeridian = 15.0 + 48.0 / 60.0 + 22.624306 / 3600.0;
        var scale = 1.00000561024;
        var falseNorthing = -667.711;
        var falseEasting = 1500064.274;

        return GridToGeodetic(x, y,
                              axis,
                              flattening,
                              centralMeridian,
                              scale,
                              falseNorthing,
                              falseEasting);
    }

    public static RT90Coordinate ToRT90(double latitude, double longitude)
    {
        var axis = 6378137.0;
        var flattening = 1.0 / 298.257222101;
        var centralMeridian = 15.0 + 48.0 / 60.0 + 22.624306 / 3600.0;
        var scale = 1.00000561024;
        var falseNorthing = -667.711;
        var falseEasting = 1500064.274;

        return GeodeticToGrid(latitude,
                              longitude,
                              axis,
                              flattening,
                              centralMeridian,
                              scale,
                              falseNorthing,
                              falseEasting);
    }

    internal static GeoCoordinate GridToGeodetic(double x,
                                                 double y,
                                                 double axis,
                                                 double flattening,
                                                 double centralMeridian,
                                                 double scale,
                                                 double falseNorthing,
                                                 double falseEasting)
    {
        double e2 = flattening * (2.0 - flattening);
        double n = flattening / (2.0 - flattening);
        double a_roof = axis / (1.0 + n) * (1.0 + n * n / 4.0 + n * n * n * n / 64.0);
        double delta1 = n / 2.0 - 2.0 * n * n / 3.0 + 37.0 * n * n * n / 96.0 - n * n * n * n / 360.0;
        double delta2 = n * n / 48.0 + n * n * n / 15.0 - 437.0 * n * n * n * n / 1440.0;
        double delta3 = 17.0 * n * n * n / 480.0 - 37 * n * n * n * n / 840.0;
        double delta4 = 4397.0 * n * n * n * n / 161280.0;

        double Astar = e2 + e2 * e2 + e2 * e2 * e2 + e2 * e2 * e2 * e2;
        double Bstar = -(7.0 * e2 * e2 + 17.0 * e2 * e2 * e2 + 30.0 * e2 * e2 * e2 * e2) / 6.0;
        double Cstar = (224.0 * e2 * e2 * e2 + 889.0 * e2 * e2 * e2 * e2) / 120.0;
        double Dstar = -(4279.0 * e2 * e2 * e2 * e2) / 1260.0;

        double deg_to_rad = Math.PI / 180;
        double lambda_zero = centralMeridian * deg_to_rad;
        double xi = (x - falseNorthing) / (scale * a_roof);
        double eta = (y - falseEasting) / (scale * a_roof);

        double xi_prim = xi -
                        delta1 * Math.Sin(2.0 * xi) * CosH(2.0 * eta) -
                        delta2 * Math.Sin(4.0 * xi) * CosH(4.0 * eta) -
                        delta3 * Math.Sin(6.0 * xi) * CosH(6.0 * eta) -
                        delta4 * Math.Sin(8.0 * xi) * CosH(8.0 * eta);

        double eta_prim = eta -
                        delta1 * Math.Cos(2.0 * xi) * SinH(2.0 * eta) -
                        delta2 * Math.Cos(4.0 * xi) * SinH(4.0 * eta) -
                        delta3 * Math.Cos(6.0 * xi) * SinH(6.0 * eta) -
                        delta4 * Math.Cos(8.0 * xi) * SinH(8.0 * eta);

        double phi_star = Math.Asin(Math.Sin(xi_prim) / CosH(eta_prim));
        double delta_lambda = Math.Atan(SinH(eta_prim) / Math.Cos(xi_prim));
        double lon_radian = lambda_zero + delta_lambda;

        double lat_radian = phi_star + Math.Sin(phi_star) * Math.Cos(phi_star) *
                        (Astar +
                         Bstar * Math.Pow(Math.Sin(phi_star), 2) +
                         Cstar * Math.Pow(Math.Sin(phi_star), 4) +
                         Dstar * Math.Pow(Math.Sin(phi_star), 6));

        var newLatitude = lat_radian * 180.0 / Math.PI;
        var newLongitude = lon_radian * 180.0 / Math.PI;

        return new GeoCoordinate(newLatitude, newLongitude);
    }

    internal static RT90Coordinate GeodeticToGrid(double latitude,
                                                  double longitude,
                                                  double axis,
                                                  double flattening,
                                                  double centralMeridian,
                                                  double scale,
                                                  double falseNorthing,
                                                  double falseEasting)
    {
        double e2 = flattening * (2.0 - flattening);
        double n = flattening / (2.0 - flattening);
        double a_roof = axis / (1.0 + n) * (1.0 + n * n / 4.0 + n * n * n * n / 64.0);
        double A = e2;
        double B = (5.0 * e2 * e2 - e2 * e2 * e2) / 6.0;
        double C = (104.0 * e2 * e2 * e2 - 45.0 * e2 * e2 * e2 * e2) / 120.0;
        double D = (1237.0 * e2 * e2 * e2 * e2) / 1260.0;
        double beta1 = n / 2.0 - 2.0 * n * n / 3.0 + 5.0 * n * n * n / 16.0 + 41.0 * n * n * n * n / 180.0;
        double beta2 = 13.0 * n * n / 48.0 - 3.0 * n * n * n / 5.0 + 557.0 * n * n * n * n / 1440.0;
        double beta3 = 61.0 * n * n * n / 240.0 - 103.0 * n * n * n * n / 140.0;
        double beta4 = 49561.0 * n * n * n * n / 161280.0;

        double deg_to_rad = Math.PI / 180.0;
        double phi = latitude * deg_to_rad;
        double lambda = longitude * deg_to_rad;
        double lambda_zero = centralMeridian * deg_to_rad;

        double phi_star = phi - Math.Sin(phi) * Math.Cos(phi)
                        * (A + B * Math.Pow(Math.Sin(phi), 2)
                             + C * Math.Pow(Math.Sin(phi), 4)
                             + D * Math.Pow(Math.Sin(phi), 6));

        double delta_lambda = lambda - lambda_zero;
        double xi_prim = Math.Atan(Math.Tan(phi_star) / Math.Cos(delta_lambda));
        double eta_prim = ATanH(Math.Cos(phi_star) * Math.Sin(delta_lambda));

        double x = scale * a_roof * (xi_prim +
                        beta1 * Math.Sin(2.0 * xi_prim) * CosH(2.0 * eta_prim) +
                        beta2 * Math.Sin(4.0 * xi_prim) * CosH(4.0 * eta_prim) +
                        beta3 * Math.Sin(6.0 * xi_prim) * CosH(6.0 * eta_prim) +
                        beta4 * Math.Sin(8.0 * xi_prim) * CosH(8.0 * eta_prim)) +
                        falseNorthing;

        double y = scale * a_roof * (eta_prim +
                        beta1 * Math.Cos(2.0 * xi_prim) * SinH(2.0 * eta_prim) +
                        beta2 * Math.Cos(4.0 * xi_prim) * SinH(4.0 * eta_prim) +
                        beta3 * Math.Cos(6.0 * xi_prim) * SinH(6.0 * eta_prim) +
                        beta4 * Math.Cos(8.0 * xi_prim) * SinH(8.0 * eta_prim)) +
                        falseEasting;

        var newX = Math.Round(x * 1000.0) / 1000.0;
        var newY = Math.Round(y * 1000.0) / 1000.0;

        return new RT90Coordinate(newX, newY);
    }

    internal static double SinH(double angle)
    {
        return 0.5 * (Math.Exp(angle) - Math.Exp(-angle));
    }

    internal static double CosH(double angle)
    {
        return 0.5 * (Math.Exp(angle) + Math.Exp(-angle));
    }

    internal static double ATanH(double angle)
    {
        return 0.5 * Math.Log((1.0 + angle) / (1.0 - angle));
    }
}