基本比例尺分幅、图幅编号、高斯投影

近日工程项目涉及基本比例尺分幅，用了CASS11的相关功能，发现有些缺陷，就打算自己整一个，期间遇到些问题，记载下来备查。

cass11见上图，不妥之处有以下几点：
1.界面比较乱，条理不清。
2.“设置中央线”没有必要，按规范即可。
3.缺失10万、25万、50万比例尺。
4.无大地坐标输出方案。

重新设计后的程序界面见下图：

分为参数设置、数据录入两个区，参数设置全部用选择框，界面简洁、清晰。程序编写过程中遇到的问题及处理方式：
1.当输入的生成范围跨多个带区时，采用大地坐标没有影响，采用投影坐标需分别计算每个带区的中央子午线及其经区范围，注意处理范围左边缘带区的最小经度、右边缘带区的最大经度范围。另外3度带的50万，100万比例尺的中央线也必须按6度带进行计算，这样中央子午线两边的投影变形一致。
2.高斯投影正算就不想浪费时间了，在CSDN检索到“氧气女神”的代码，程序结构很棒，BUT,计算结果是错的，私信了也不回话，这妹子估计比较忙。这是用她的数据及她的程序计算结果。“https://blog.csdn.net/qq_44718128/article/details/117874196”

后来检索到“卤雅少年”的代码“https://blog.csdn.net/Zj1638/article/details/125380593”，将两人的代码整合了一下，得到高斯投影正算代码，精度为0.0005，先将就用吧，哪位兄台有更高精度的代码请发一个，谢谢！我的邮箱：huaheroys@126.com

       #region 高斯投影正算 纬度，经度，中央子午线，代号，是否加带号大数 单位度  
        public static double[] BLtoXY(double B,double L,double L0,double dh,bool jdh)        {

            Editor ed = Autodesk.AutoCAD.ApplicationServices.Application.DocumentManager.MdiActiveDocument.Editor;
            //CGCS2000椭球参数
            double a = 6378137.0;
            double f = 0.00335281066475;
            double e1 = 0.00669437999014631;
            double e2 = 0.00673949674228149;
            double m0, m2, m4, m6, m8;
            double a0, a2, a4, a6, a8;       
         
           //m0-m8的值
            m0 = a * (1 - e1);
            m2 = 1.5 * m0 * e1;
            m4 = 1.25 * m2 * e1;
            m6 = (7.0 / 6) * m4 * e1;
            m8 = (9.0 / 8) * m6 * e1;
            //a0-a8的值
            a0 = m0 + m2 / 2.0 + (3.0 / 8) * m4 + (5.0 / 16) * m6 + (35.0 / 128) * m8;
            a2 = (1.0 / 2) * (m2 + m4) + (15.0 / 32) * m6 + (7.0 / 16) * m8;
            a4 = (1.0 / 8) * m4 + (3.0 / 16) * m6 + (7.0 / 32) * m8;
            a6 = (1.0 / 32) * m6 + (1.0 / 16) * m8;
            a8 = (1.0 / 128) * m8;

            //求椭球面上某点到赤道的子午线弧长，精度至0.001m
            double radB = (Math.PI / 180) * B; //将纬度转换为弧度
            double Δ1, Δ2, Δ3, Δ4, Δ5;
            Δ1 = a0 * radB;
            Δ2 = (a2 / 2) * Math.Sin(2 * radB);
            Δ3 = (a4 / 4) * Math.Sin(4 * radB);
            Δ4 = (a6 / 6) * Math.Sin(6 * radB);
            Δ5 = (a8 / 8) * Math.Sin(8 * radB);
            double X = Math.Round(Δ1 - Δ2 + Δ3 - Δ4 + Δ5, 8);
            //卯酉圈曲率半径N
            double N = a / Math.Sqrt(1 - e1 * Math.Sin(radB) * Math.Sin(radB));
            // tan(B)的值
            double t = Math.Tan(radB);
            // η的值
            double η = Math.Sqrt(e2) * Math.Cos(radB);
            //l的值
            double ρm = 180 * 60 * 60 / Math.PI;
            double l = ((L - L0) * 3600) / ρm;
            //高斯坐标中的x值和y值
            double x = Math.Round(X + (N / 2) * t * Math.Cos(radB) * Math.Cos(radB) * l * l + (N / 24) * t * (5 - t * t + 9 * η * η + 4 * Math.Pow(η, 4)) * Math.Pow(Math.Cos(radB), 4) * Math.Pow(l, 4) + (N / 720) * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(l, 6) * (Math.Pow(Math.Cos(radB), 6)), 4);
            double y = 500000 + Math.Round(N * Math.Cos(radB) * l + (N / 6.0) * (1 - t * t + η * η) * Math.Pow(Math.Cos(radB), 3) * Math.Pow(l, 3) + (N / 120.0) * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * η * η - 58 * η * η * t * t) * Math.Pow(Math.Cos(radB), 5) * Math.Pow(l, 5), 4);

            if (jdh) { y = y + dh * 1000000; }         
            return new double[] { y, x };
        }

        #endregion

程序运行效果-大地坐标

程序运行效果-投影坐标

欢迎批评指导，QQ群：861650360.