墨卡托投影C#实现

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 = 10; //反向转换程序中的迭代次数
        static public double __IterativeValue = 0; //反向转换程序中的迭代初始值
        static public double __A = 6378137.0; //椭球体长轴,米
        static public double __B = 6356752.314; //椭球体短轴,米
        static public double __B0 = 0; //标准纬度,弧度
        static public double __L0 = 0; //原点经度,弧度

        //角度到弧度的转换
        static public double DegreeToRad(double degree)
        {
            return Math.PI * ((double)degree / (double)180);
        }

        //弧度到角度的转换
        static public double RadToDegree(double rad)
        {
            return (180 * rad) / Math.PI;
        }

        //以上参数的设定由如下几个public函数完成
        //设定__A与__B
        static public void SetAB(double a, double b)
        {
            if (a <= 0 || b <= 0)
            {
                return;
            }
            __A = a;
            __B = b;
        }
        //设定__B0
        static public void SetB0(double b0)
        {
            b0 = DegreeToRad(b0);
            if (b0 < -Math.PI / 2 || b0 > Math.PI / 2)
            {
                return;
            }
            __B0 = b0;
        }

        //设定__L0
        static public void SetL0(double l0)
        {
            l0 = DegreeToRad(l0);
            if (l0 < -Math.PI || l0 > Math.PI)
            {
                return;
            }
            __L0 = l0;
        }

        /*******************************************
        投影正向转换程序
        double B: 纬度,角度
        double L: 经度,角度
        double& X: 纵向直角坐标
        double& Y: 横向直角坐标
        *******************************************/
        //static public int ToProj(double B, double L, ref double X, ref double Y)
        static public PointXY LBToXY(PointLB pmtLonLat0, PointLB pmtLonLat1)
        {
            SetB0(pmtLonLat0.lat);
            SetL0(pmtLonLat0.lon);

            double B = DegreeToRad(pmtLonLat1.lat);
            double L = DegreeToRad(pmtLonLat1.lon);

            PointXY xy = new PointXY();
            xy.x = 0; xy.y = 0;

            double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp;
            double E = Math.Exp(1);
            if (L < -Math.PI || L > Math.PI || B < -Math.PI / 2 || B > Math.PI / 2)
            {
                return xy;
            }
            if (__A <= 0 || __B <= 0)
            {
                return xy;
            }
            f = (__A - __B) / __A;
            dtemp = 1 - (__B / __A) * (__B / __A);
            if (dtemp < 0)
            {
                return xy;
            }
            e = Math.Sqrt(dtemp);
            dtemp = (__A / __B) * (__A / __B) - 1;
            if (dtemp < 0)
            {
                return xy;
            }
            e_ = Math.Sqrt(dtemp);
            NB0 = ((__A * __A) / __B) / Math.Sqrt(1 + 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 / 4 + B / 2) * Math.Pow((1 - e * Math.Sin(B)) / (1 + e * Math.Sin(B)), e / 2));
            
            return xy;
        }

        /*******************************************
        投影反向转换程序
        double X: 纵向直角坐标
        double Y: 横向直角坐标
        double& B: 维度,弧度
        double& L: 经度,弧度
        *******************************************/
        //static public int FromProj(double X, double Y, ref double B, ref double L)
        static public PointLB XYtoLB(PointLB pmtLB0, PointXY pmtXY)
        {
            SetB0(pmtLB0.lat);
            SetL0(pmtLB0.lon);

            double X = pmtXY.x;
            double Y = pmtXY.y;

            double B = 0, L = 0;

            PointLB lb = new PointLB();
            lb.lat = 0; lb.lon = 0;

            double f/*扁率*/, e/*第一偏心率*/, e_/*第二偏心率*/, NB0/*卯酉圈曲率半径*/, K, dtemp;
            double E = Math.Exp(1);
            if (__A <= 0 || __B <= 0)
            {
                return lb;
            }
            f = (__A - __B) / __A;
            dtemp = 1 - (__B / __A) * (__B / __A);
            if (dtemp < 0)
            {
                return lb;
            }
            e = Math.Sqrt(dtemp);
            dtemp = (__A / __B) * (__A / __B) - 1;
            if (dtemp < 0)
            {
                return lb;
            }
            e_ = Math.Sqrt(dtemp);
            NB0 = ((__A * __A) / __B) / Math.Sqrt(1 + e_ * e_ * Math.Cos(__B0) * Math.Cos(__B0));
            K = NB0 * Math.Cos(__B0);

            L = X / K + __L0;
            B = __IterativeValue;

            for (int i = 0; i < __IterativeTimes; i++)
            {
                B = Math.PI / 2 - 2 * Math.Atan(Math.Pow(E, (-Y / K)) * Math.Pow(E, (e / 2) * Math.Log((1 - e * Math.Sin(B)) / (1 + e * Math.Sin(B)))));
            }

            lb.lon = RadToDegree(L);
            lb.lat = RadToDegree(B);

            return lb;
        }
    }
}