c#求大地坐标正算

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

一、程序简介

大地坐标正算是根据给定的地球椭球体参数和起点经纬度、方位角和距离计算出终点的经纬度坐标

二、程序实现逻辑

定义地球椭球体参数
定义WGS84椭球体长半轴a=6378137.0，短半轴b=6356752.3142，第一偏心率e=(a^2-b2)/a=0.0818191908426。

输入起点的大地坐标和方位角距离
输入起点的大地坐标(lat1, lon1)和方位角bear和距离s。

计算辅助量
先将起点纬度换算为弧度，计算出扁率f=(a-b)/a，以及纬度圆的半径rho=a*(1-e^2)/(1-e2sin^2(lat1))1.5，扁率圆的半径nu=a/sqrt(1-e^2sin2(lat1))，然后计算出由起点沿方位角bear前进距离s到达的点的纬度lat2和初始子午线弧长度A=0。

迭代计算
根据公式计算纬度lat2、经差dLon、子午线弧曲率半径C、球面距离s和初始子午线弧长A，直到dLon小于1e-12。

计算终点经度lon2
根据公式lon2=lon1+dLon。

三、程序代码

以下是C#代码实现大地坐标正算：

using System;

namespace GeodeticCalculation
{
   
    public class Ellipsoid
    {
   
        public const double WGS84_A = 6378137.0;
        public const double WGS84_B = 6356752.3142;
        public const double WGS84_F = (WGS84_A - WGS84_B) / WGS84_A;

        public double a;
        public double b;
        public double f;

        public Ellipsoid(double a, double b)
        {
   
            this.a = a;
            this.b = b;
            f = (a - b) / a;
        }

        public double GetRho(double lat)
        {
   
            double sinLat = Math.Sin(lat);
            return a * (1 - f) / Math.Pow(1 - f * sinLat * sinLat, 1.5);
        }

        public double GetNu(double lat)
        {
   
            double sinLat = Math.Sin(lat);
            return a / Math.Sqrt(1 - f * sinLat * sinLat);
        }
    }

    public class GeodeticCalculator
    {
   
        public static