using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Threading.Tasks;

namespace MRP

{

    public class ClassMct

    {

        static public int __IterativeTimes = ; //反向转换程序中的迭代次数

        static public double __IterativeValue = ; //反向转换程序中的迭代初始值

        static public double __A = 6378.137; //椭球体长轴,千米

        static public double __B = 6356.752314; //椭球体短轴,千米

        static public double __B0 = ; //标准纬度,弧度

        static public double __L0 = ; //原点经度,弧度

        //角度到弧度的转换

        static public double DegreeToRad(double degree)

        {

            return Math.PI * degree / ;

        }

        //弧度到角度的转换

        static public double RadToDegree(double rad)

        {

            return ( * rad) / Math.PI;

        }

        //设定__A与__B

        static public void SetAB(double a, double b)

        {

            if (a <=  || b <= )

            {

                return;

            }

            __A = a;

            __B = b;

        }

        //设定__B0

        static public void SetLB0(double pmtL0, double pmtB0)

        {

            double l0 = DegreeToRad(pmtL0);

            if (l0 < -Math.PI || l0 > Math.PI)

            {

                return;

            }

            __L0 = l0;

            double b0 = DegreeToRad(pmtB0);

            if (b0 < -Math.PI /  || b0 > Math.PI / )

            {

                return;

            }

            __B0 = b0;

        }

        /*******************************************

        经纬度转XY坐标

        pmtLB0: 参考点经纬度

        pmtLB1: 要转换的经纬度

        返回值: 直角坐标，单位：公里

        *******************************************/

        static public PointXY LBToXY(PointLB pmtLB0, PointLB pmtLB1)

        {

            SetLB0(pmtLB0.lon, pmtLB0.lat);

            double B = DegreeToRad(pmtLB1.lat);

            double L = DegreeToRad(pmtLB1.lon);

            PointXY xy = new PointXY();

            xy.x = ; xy.y = ;

            double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp;

            double E = Math.Exp();

            if (L < -Math.PI || L > Math.PI || B < -Math.PI /  || B > Math.PI / )

            {

                return xy;

            }

            if (__A <=  || __B <= )

            {

                return xy;

            }

            f = (__A - __B) / __A;

            dtemp =  - (__B / __A) * (__B / __A);

            if (dtemp < )

            {

                return xy;

            }

            e = Math.Sqrt(dtemp);

            dtemp = (__A / __B) * (__A / __B) - ;

            if (dtemp < )

            {

                return xy;

            }

            e_ = Math.Sqrt(dtemp);

            NB0 = ((__A * __A) / __B) / Math.Sqrt( + e_ * e_ * Math.Cos(__B0) * Math.Cos(__B0));

            K = NB0 * Math.Cos(__B0);

            xy.x = K * (L - __L0);

            xy.y = K * Math.Log(Math.Tan(Math.PI /  + (B) / ) * Math.Pow(( - e * Math.Sin(B)) / ( + e * Math.Sin(B)), e / ));

            double y0 = K * Math.Log(Math.Tan(Math.PI /  + (__B0) / ) * Math.Pow(( - e * Math.Sin(__B0)) / ( + e * Math.Sin(__B0)), e / ));

            xy.y = xy.y - y0;

            xy.y = -xy.y;//正常的Y坐标系（向上）转程序的Y坐标系（向下）

            return xy;

        }

        /*******************************************

        XY坐标转经纬度

        pmtLB0: 参考点经纬度

        pmtXY: 要转换的XY坐标，单位：公里

        返回值: 经纬度

        *******************************************/

        static public PointLB XYtoLB(PointLB pmtLB0, PointXY pmtXY)

        {

            SetLB0(pmtLB0.lon, pmtLB0.lat);

            double X = pmtXY.x;

            double Y = -pmtXY.y;//程序的Y坐标系（向下）转正常的Y坐标系（向上）    

            double B = , L = ;

            PointLB lb = new PointLB();

            lb.lat = ; lb.lon = ;

            double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp;

            double E = Math.Exp();

            if (__A <=  || __B <= )

            {

                return lb;

            }

            f = (__A - __B) / __A;

            dtemp =  - (__B / __A) * (__B / __A);

            if (dtemp < )

            {

                return lb;

            }

            e = Math.Sqrt(dtemp);

            dtemp = (__A / __B) * (__A / __B) - ;

            if (dtemp < )

            {

                return lb;

            }

            e_ = Math.Sqrt(dtemp);

            NB0 = ((__A * __A) / __B) / Math.Sqrt( + e_ * e_ * Math.Cos(__B0) * Math.Cos(__B0));

            K = NB0 * Math.Cos(__B0);

            double y0 = K * Math.Log(Math.Tan(Math.PI /  + (__B0) / ) * Math.Pow(( - e * Math.Sin(__B0)) / ( + e * Math.Sin(__B0)), e / ));

            Y = Y + y0;

            L = X / K + __L0;

            B = __IterativeValue;

            for (int i = ; i < __IterativeTimes; i++)

            {

                B = Math.PI /  -  * Math.Atan(Math.Pow(E, (-Y / K)) * Math.Pow(E, (e / ) * Math.Log(( - e * Math.Sin(B)) / ( + e * Math.Sin(B)))));

            }

            lb.lon = RadToDegree(L);

            lb.lat = RadToDegree(B);

            return lb;

        }

    }

}