//-----------------------------------------------------------------------
// <copyright file="RdLatLngConverter.cs" company="Deltares">
//  Copyright (c) 2010 Deltares. All rights reserved.
// </copyright>
// <author>B.S.T. The</author>
// <email>tom.the@deltares.nl</email>
// <date>04-11-2010</date>
// <summary>n.a.</summary>
//-----------------------------------------------------------------------
using Deltares.Geometry;

namespace Deltares.Dam.Data
{
    using System;

    public class RdLatLngConverter
    {
        private const double A01 = 3236.0331637;
        private const double A02 = -0.2472814;
        private const double A03 = -0.0655238;
        private const double A04 = 0.0000371;
        private const double A20 = -32.5915821;
        private const double A21 = -0.8501341;
        private const double A22 = -0.0171137;
        private const double A23 = -0.0003859;
        private const double A24 = -0.0000090;
        private const double A40 = 0.0052771;
        private const double A41 = 0.0003314;
        private const double A42 = 0.0000143;

        private const double B10 = 5261.3028966;
        private const double B11 = 105.9780241;
        private const double B12 = 2.4576469;
        private const double B13 = 0.0560089;
        private const double B14 = 0.0012770;
        private const double B15 = 0.0000291;
        private const double B30 = -0.8192156;
        private const double B31 = -0.0560092;
        private const double B32 = -0.0025614;
        private const double B33 = -0.0000973;
        private const double B50 = 0.0002574;
        private const double B51 = 0.0000293;

        private const double C01 = 190066.98903;
        private const double C03 = -32.38360;
        private const double C05 = -0.00661;
        private const double C11 = -11830.85831;
        private const double C13 = -0.60639;
        private const double C21 = -114.19754;
        private const double C23 = 0.15774;
        private const double C31 = -2.34078;
        private const double C41 = -0.04158;
        private const double D02 = 3638.36193;
        private const double D04 = 0.09351;
        private const double D10 = 309020.31810;
        private const double D12 = -157.95222;
        private const double D14 = -0.05419;
        private const double D20 = 72.97141;
        private const double D22 = -6.43481;
        private const double D30 = 59.79734;
        private const double D32 = -0.07379;
        private const double D40 = -0.03444;
        private const double Lam0 = 5.387638889;
        private const double Phi0 = 52.156160556;

        private const double X0 = 155000.00;
        private const double Y0 = 463000.00;

        public static GeometryPoint RdToLatLng(double xRd, double yRd)
        {
            var phi = GetPhi(xRd, yRd);
            var lam = GetLam(xRd, yRd);
            var lat = GetLat(lam, phi);
            var lng = GetLng(lam, phi);

            return new GeometryPoint() {X = lat, Y = lng};
        }

        public static GeometryPoint LatLngToRd(double latitude, double longitude)
        {
            var x = GetX(latitude, longitude);
            var y = GetY(latitude, longitude);
            
            return new GeometryPoint() {X = x, Y = y};
        }


        private static double GetLng(double lam, double phi)
        {
            return lam + (-37.902 + 0.329*(phi - 52) - 14.667*(lam - 5))/100000;
        }

        private static double GetLat(double lam, double phi)
        {
            return phi + (-96.862 - 11.714*(phi - 52) - 0.125*(lam - 5))/100000;
        }

        private static double GetPhi(double x, double y)
        {
            var dX = (x - X0)*Math.Pow(10, -5);
            var dY = (y - Y0)*Math.Pow(10, -5);
            var dPhi = A01*dY + A20*Math.Pow(dX, 2) + A02*Math.Pow(dY, 2) + A21*Math.Pow(dX, 2)*dY +
                          A03*Math.Pow(dY, 3) + A40*Math.Pow(dX, 4) +
                          A22*Math.Pow(dX, 2)*Math.Pow(dY, 2) + A04*Math.Pow(dY, 4) + A41*Math.Pow(dX, 4)*dY +
                          A23*Math.Pow(dX, 2)*Math.Pow(dY, 3) +
                          A42*Math.Pow(dX, 4)*Math.Pow(dY, 2) + A24*Math.Pow(dX, 2)*Math.Pow(dY, 4);
            return Phi0 + (dPhi/3600.00);
        }

        private static double GetLam(double x, double y)
        {
            var dX = (x - X0)*Math.Pow(10, -5);
            var dY = (y - Y0)*Math.Pow(10, -5);
            var dLam = B10*dX + B11*dX*dY + B30*Math.Pow(dX, 3) + B12*dX*Math.Pow(dY, 2) + B31*Math.Pow(dX, 3)*dY +
                          B13*dX*Math.Pow(dY, 3) +
                          B50*Math.Pow(dX, 5) + B32*Math.Pow(dX, 3)*Math.Pow(dY, 2) + B14*dX*Math.Pow(dY, 4) +
                          B51*Math.Pow(dX, 5)*dY +
                          B33*Math.Pow(dX, 3)*Math.Pow(dY, 3) + B15*dX*Math.Pow(dY, 5);
            return Lam0 + (dLam/3600.00);
        }

        private static double GetX(double lat, double lng)
        {
            var dPhi = (lat - Phi0)*0.36;
            var dLam = (lng - Lam0)*0.36;
            var dX = C01*dLam + C11*dPhi*dLam + C21*Math.Pow(dPhi, 2)*dLam + C03*Math.Pow(dLam, 3) +
                        C31*Math.Pow(dPhi, 3)*dLam + C13*dPhi*Math.Pow(dLam, 3) +
                        C23*Math.Pow(dPhi, 2)*Math.Pow(dLam, 3) +
                        C41*Math.Pow(dPhi, 4)*dLam + C05*Math.Pow(dLam, 5);
            return X0 + dX;
        }

        private static double GetY(double lat, double lng)
        {
            var dPhi = (lat - Phi0)*0.36;
            var dLam = (lng - Lam0)*0.36;
            var dY = D10*dPhi + D20*Math.Pow(dPhi, 2) + D02*Math.Pow(dLam, 2) + D12*dPhi*Math.Pow(dLam, 2) +
                        D30*Math.Pow(dPhi, 3) + D22*Math.Pow(dPhi, 2)*Math.Pow(dLam, 2) + D40*Math.Pow(dPhi, 4) +
                        D04*Math.Pow(dLam, 4) +
                        D32*Math.Pow(dPhi, 3)*Math.Pow(dLam, 2) + D14*dPhi*Math.Pow(dLam, 4);
            return Y0 + dY;
        }
    }
}