墨卡托投影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 = 6378.137; //椭球体长轴,千米
        static public double __B = 6356.752314; //椭球体短轴,千米
        static public double __B0 = 0; //标准纬度,弧度
        static public double __L0 = 0; //原点经度,弧度

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

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

        //设定__A与__B
        static public void SetAB(double a, double b)
        {
            if (a <= 0 || b <= 0)
            {
                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 / 2 || b0 > Math.PI / 2)
            {
                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 = 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));
            double y0 = K * Math.Log(Math.Tan(Math.PI / 4 + (__B0) / 2) * Math.Pow((1 - e * Math.Sin(__B0)) / (1 + e * Math.Sin(__B0)), e / 2));
            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 = 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);

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

            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;
        }
    }
}