关于经纬度算法

  public Angle se(Angle latA, Angle lonA, TKsoft.Earth.Angle latB, Angle lonB)
        {
            //知道 两点经纬度，求角度的方法
            double cosLatB = Math.Cos(latB.Radians);
            Angle tcA = TKsoft.Earth.Angle.FromRadians(Math.Atan2(
                Math.Sin(lonA.Radians - lonB.Radians) * cosLatB,
                Math.Cos(latA.Radians) * Math.Sin(latB.Radians) -
                Math.Sin(latA.Radians) * cosLatB *
                Math.Cos(lonA.Radians - lonB.Radians)));
            if (tcA.Radians < 0)
                tcA.Radians = tcA.Radians + Math.PI * 2;
            tcA.Radians = Math.PI * 2 - tcA.Radians;
            return tcA;
        }

以上代码因项目问题，angle是一个类，各位以此代码 做为参考应该能猜出个大概来（Angle 前缀引用 被 我了）

 /**
         *根据一点、角点(正北0度，顺生针增大)、距离(米)计算另一点
        */
        public static System.Drawing.PointF CalPoint(double x, double y, double angel, double distance)
        {
            //
            System.Drawing.PointF endPoint = new System.Drawing.PointF();
            //角度换算成正东0度，逆时针增大
            angel = -angel - 360+90;
            //米换算成经纬度
            var dis = distance / 1000 / 111.7;
            //角点换成弧度
            angel = angel * Math.PI / 180;
            //计算y坐标
            endPoint.Y = (float)(Math.Sin(angel) * dis + y);
            //计算x坐标
            endPoint.X = (float)(Math.Cos(angel) * dis + x);
            return endPoint;
        }

    /// <summary>
        /// 计算两个经纬度之间的经纬度差
        /// </summary>
        /// <param name="pt1"></param>
        /// <param name="pt2"></param>
        /// <returns></returns>
        public static double Distance(IPoint pt1, IPoint pt2)
        {
            return Distance(pt1.X, pt1.Y, pt2.X, pt2.Y);
        }
        /// <summary>
        /// 计算两个坐标点之间的距离
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="tx"></param>
        /// <param name="ty"></param>
        /// <returns></returns>
        public static double Distance(double x, double y, double tx, double ty)
        {
            double a = ty - y;
            double b = tx - x;
            return Math.Sqrt(a * a + b * b);
        }

        /// <summary>
        /// 计算起始点和结束点之间的地理方向
        /// </summary>
        /// <param name="start"></param>
        /// <param name="end"></param>
        /// <returns></returns>
        public static double CalcDirection(IPoint start, IPoint end)
        {
            return CalcDirection(start.X, start.Y, end.X, end.Y);
        }


    public static double CalcDirection(double bX, double bY, double eX, double eY)
        {
            double degree = 0;
            double dx = eX - bX;
            double dy = eY - bY;
            if (dx == 0)
            {
                if (bY < eY) degree = 0;
                else degree = 180;
            }
            else
            {
                degree = Math.Atan(dy / dx) * 180 / Math.PI;
                if (degree < 0 && bY < eY) degree = 180 + degree;
                degree = 90 - degree;
                if (bY > eY && bX > eX) degree += 180;
                if (degree == 90 && bX > eX) degree += 180;
            }
            return degree;
        }

    /// <summary>
        /// 获取指定纬度线上纬线圈的半径
        /// </summary>
        /// <param name="lat"></param>
        /// <returns></returns>
        public static double CalcLatRadius(double lat)
        {
            double arc = lat / 180 * Math.PI;
            return Math.Cos(arc) * DrawArgs.EquatorialRadius;
        }

    /// <summary>
        /// 求出指定长度的线段在指定纬度圈上的经度跨度
        /// </summary>
        /// <param name="length">长度(或者距离，单位:米)</param>
        /// <param name="lat">纬度</param>
        /// <returns></returns>
        public static double CalcLonSpan(double length, double lat)
        {
            double c = 2 * Math.PI * CalcLatRadius(lat);
            return length / c * 360;
        }

/// <summary>
        /// 计算指定长度的线段在纬度跨度
        /// </summary>
        /// <param name="length"></param>
        /// <returns></returns>
        public static double CalcLatSpan(double length)
        {
            double c = 2 * Math.PI * DrawArgs.EquatorialRadius;
            return length / c * 360;
        }

/// <summary>
        /// 计算num的pwo次幂
        /// </summary>
        /// <param name="num"></param>
        /// <param name="pow"></param>
        /// <returns></returns>
        public static double Pow(double num, int pow)
        {
            double result = Convert.ToDouble(num);
            for (int i = 1; i < pow; i++)
            {
                result *= num;
            }
            return result;
        }

   //定义角度转化为弧度的函数
        public double angle_radian(double[] angle)                                       //角度转化为弧度函数
        {
            double radian;
            radian = (angle[0] + angle[1] / 60 + angle[2] / 3600) / 180 * PI;
            return radian;
        }

        //定义弧度转化为角度的函数
        public double[] radian_angle(double radian)                                      //弧度转化为角度函数
        {
            double degree, minute, second;
            double[] angle = { 0, 0, 0 };
            degree = System.Math.Truncate(radian /PI * 180);
            minute = System.Math.Truncate((radian / PI * 180 - degree) * 60);
            second = System.Math.Abs(((radian / PI * 180 - degree) * 60 - minute) * 60);
            minute = System.Math.Abs(minute);
            angle[0] = degree;
            angle[1] = minute;
            angle[2] = second;
            return angle;
        }

 public static double E = 0.081813334,PI = Math.PI;              //定义常变量
         public static double B1, L1, A1, B2, L2, A2, S;                 //定义主变量
         public static double A,B,C,α,β;                              //定义正算变量
         public static double a1, a2, b1, b2;                           //定义反算变量


//大地主题正算
   //大地主题正算
        private void button1_Click(object sender, EventArgs e)
        { 
             double[] Bre1={0,0,0},Lat1={0,0,0},Ang1={0,0,0}, Bre2={0,0,0},Lat2={0,0,0},Ang2={0,0,0};
             double W1, sinu1, sinu2, cosu1, sinA0, cotσ1, sin2σ1, cos2σ1,couAo, sin2σ1plusσ0, cos2σ1plusσ0, σ, σ0, λ, δ;

             Bre1[0] = Convert.ToDouble(textBox1.Text); Bre1[1] = Convert.ToDouble(textBox2.Text);Bre1[2] = Convert.ToDouble(textBox3.Text);
             Lat1[0] = Convert.ToDouble(textBox4.Text); Lat1[1] = Convert.ToDouble(textBox5.Text);Lat1[2] = Convert.ToDouble(textBox6.Text);
             Ang1[0] = Convert.ToDouble(textBox7.Text); Ang1[1] = Convert.ToDouble(textBox8.Text); Ang1[2] = Convert.ToDouble(textBox9.Text);
             S= Convert.ToDouble(textBox10.Text);
             B1= angle_radian(Bre1);L1= angle_radian(Lat1);A1= angle_radian(Ang1);

            //计算起点规划纬度
             W1 = Math.Sqrt(1.0 - Math.Pow(E *Math.Sin(B1), 2));
             sinu1 = (Math.Sin(B1) * Math.Sqrt(1 - E * E)) / W1;
             cosu1 = Math.Cos(B1) / W1;

            //计算辅助函数值
             sinA0 = cosu1 * Math.Sin(A1);
             cotσ1 = cosu1 * Math.Cos(A1) / sinu1;
             sin2σ1 = cotσ1 * 2 / (cotσ1 * cotσ1 + 1);
             cos2σ1 = (cotσ1 * cotσ1 - 1) / (cotσ1 * cotσ1 + 1);

            //计算必要系数值
             couAo = Math.Sqrt(1 - sinA0 * sinA0);
             A = 6356863.020 + (10708.949 - 13.474 * couAo * couAo) * couAo * couAo;
             B = (5354.469 - 8.978 * couAo * couAo) * couAo * couAo;
             C = (2.238 * couAo * couAo) * couAo * couAo + 0.006;
             α = (33523299 - (28189 - 70 * couAo * couAo) * couAo * couAo) * Math.Pow(10, -10);
             β= (14093.5 - 48.5 * couAo * couAo) * couAo * couAo * System.Math.Pow(10, -10);

             //计算球面长度
             σ0 = (S - (B + C * cos2σ1) * sin2σ1) / A;
             sin2σ1plusσ0 = sin2σ1 * Math.Cos(2 * σ0) + cos2σ1 * Math.Sin(2 * σ0);
             cos2σ1plusσ0 = cos2σ1 * Math.Cos(2 * σ0) - sin2σ1 * Math.Sin(2 * σ0);
             σ = σ0 + (B + 5 * C * cos2σ1plusσ0) * sin2σ1plusσ0 / A;

             //计算经度差改正数
             δ = (α * σ + β * (sin2σ1plusσ0 - sin2σ1)) * sinA0;

            //计算终点大地坐标及坐标方位角
            sinu2 = sinu1 * Math.Cos(σ) + cosu1 * Math.Cos(A1) * Math.Sin(σ);
            B2 = Math.Atan(sinu2 / (Math.Sqrt(1 - E * E) * Math.Sqrt(1 - sinu2 * sinu2)));
            λ = Math.Atan(Math.Sin(A1) * Math.Sin(σ) / (cosu1 * Math.Cos(σ) - sinu1 * Math.Sin(σ) * Math.Cos(A1)));

            if (Math.Sin(A1) > 0)
            {
                if (Math.Tan(λ) > 0)
                {
                    λ = Math.Abs(λ);
                }
                else
                {
                    λ = PI - Math.Abs(λ);
                }
            }
            else
            {
                if (Math.Tan(λ) < 0)
                {
                    λ = -Math.Abs(λ);
                }
                else
                {
                    λ = Math.Abs(λ) - PI;
                }
            }
            L2 = L1 + λ - δ;
            A2 =Math.Atan(cosu1 * Math.Sin(A1) / (cosu1 * Math.Cos(σ) *Math.Cos(A1) - sinu1 *Math.Sin(σ)));
            if (Math.Sin(A1) < 0)
            {
                if (Math.Tan(A2) > 0)
                {
                    A2 = Math.Abs(A2);
                }
                else
                {
                    A2 = PI - Math.Abs(A2);
                }
            }
            else
            {
                if (Math.Tan(A2) > 0)
                {
                    A2 = PI + Math.Abs(A2);
                }
                else
                {
                    A2 = 2 * PI - Math.Abs(A2);
                }
            }
            //将弧度制角化为角度制角                                                                                              
            Bre2 = radian_angle(B2);
            Lat2 = radian_angle(L2);
            Ang2 = radian_angle(A2);


            //通过语句在程序界面显示最后结果
            textBox11.Text = Convert.ToString(Bre2[0]); textBox12.Text = Convert.ToString(Bre2[1]); textBox13.Text = Convert.ToString(Bre2[2]);
            textBox14.Text = Convert.ToString(Lat2[0]); textBox15.Text = Convert.ToString(Lat2[1]); textBox16.Text = Convert.ToString(Lat2[2]);
            textBox17.Text = Convert.ToString(Ang2[0]); textBox18.Text = Convert.ToString(Ang2[1]); textBox19.Text = Convert.ToString(Ang2[2]);

        
        }

//大地主题 反算

  private void button2_Click(object sender, EventArgs e)
        {
             double B1, L1, A1, B2, L2, A2, S;
             double[] Bre1 = { 0, 0, 0 }, Lat1 = { 0, 0, 0 }, Ang1 = { 0, 0, 0 }, Bre2 = { 0, 0, 0 }, Lat2 = { 0, 0, 0 }, Ang2 = { 0, 0, 0 };
             double W1, W2, sinu1, sinu2, cosu1, cosu2, sinA0,cosA0,  sinσ, cosσ, σ,  λ, L, δ, δ1,p, q, x, y, β1,B11,C11;

             Bre1[0] = Convert.ToDouble(textBox25.Text); Bre1[1] = Convert.ToDouble(textBox24.Text); Bre1[2] = Convert.ToDouble(textBox23.Text);
             Lat1[0] = Convert.ToDouble(textBox20.Text); Lat1[1] = Convert.ToDouble(textBox21.Text); Lat1[2] = Convert.ToDouble(textBox22.Text);
             Bre2[0] = Convert.ToDouble(textBox26.Text); Bre2[1] = Convert.ToDouble(textBox27.Text); Bre2[2] = Convert.ToDouble(textBox28.Text);
             Lat2[0] = Convert.ToDouble(textBox29.Text); Lat2[1] = Convert.ToDouble(textBox30.Text); Lat2[2] = Convert.ToDouble(textBox31.Text);
             B1 = angle_radian(Bre1); L1 = angle_radian(Lat1); 
             B2 = angle_radian(Bre2); L2 = angle_radian(Lat2);

             W1 = Math.Sqrt(1 - E * E * Math.Sin(B1) * Math.Sin(B1));
             W2 = Math.Sqrt(1 - E * E * Math.Sin(B2) * Math.Sin(B2));
             sinu1 = (Math.Sin(B1) * Math.Sqrt(1 - E * E)) / W1;
             sinu2 = (Math.Sin(B2) * Math.Sqrt(1 - E * E)) / W2;
             cosu1 = Math.Cos(B1) / W1;
             cosu2 = Math.Cos(B2) / W2;
             L = L2 - L1;
             a1 = sinu1 * sinu2;
             a2 = cosu1 * cosu2;
             b1 = cosu1 * sinu2;
             b2 = sinu1 * cosu2;

             //用逐次趋近法计算大地起点大地方位角、球面长度及经差λ=L+δ：
             δ = 0; λ = 1.0; δ1 = 0;
             do
             {
                 δ = δ1;
                 p = cosu2 * Math.Sin(λ);
                 q = b1 - b2 * Math.Cos(λ);
                 A1 = Math.Atan(p / q);
                 if (p > 0)
                 {
                     if (q > 0)
                     {
                         A1 = Math.Abs(A1);
                     }
                     else
                     {
                         A1 = PI - Math.Abs(A1);
                     }
                 }
                 else
                 {
                     if (p < 0)
                     {
                         A1 =PI + Math.Abs(A1);
                     }
                     else
                     {
                         A1 = 2 * PI - Math.Abs(A1);
                     }
                 }
                 sinσ = p * Math.Sin(A1) + q * Math.Cos(A1);
                 cosσ = a1 + a2 * Math.Cos(λ);
                 σ = Math.Atan(sinσ / cosσ);

                 if (cosσ > 0)   
                      { σ = Math.Abs(σ); } 
                 else
                      { σ = PI - Math.Abs(σ); }

                 sinA0 = cosu1 * Math.Sin(A1);
                 cosA0 = Math.Sqrt(1 - sinA0 * sinA0);
                 x = 2 * a1 - Math.Pow(cosA0, 2) * cosσ;
                α = (33523299 - (28189 - 70 * cosA0 * cosA0) * cosA0 * cosA0) * Math.Pow(10, -10);
                 β1 = (28189 - 94 * cosA0 * cosA0) * Math.Pow(10, -10);
                 δ1 = (α * σ - β1 * x *Math.Sin(σ)) * sinA0;
                 λ = L + δ1;
             } while (Math.Abs(δ - δ1) >= Math.Pow(10, -9));
             δ = δ1;

             //计算某些系数
             A = 6356863.020 + (10708.949 - 13.474 * cosA0 * cosA0) * cosA0 * cosA0;
             B11 = 10708.938 - 17.956 * cosA0 * cosA0;
             C11 = 4.487;

             //计算大地线长度S
             y = (Math.Pow(cosA0, 4) - 2 * x * x) * cosσ;
             S = A * σ + (B11 * x + C11 * y) * sinσ;

             //计算反方位角A2
             p = cosu1 * System.Math.Sin(λ); q = b1 * System.Math.Cos(λ) - b2;
             A2 = System.Math.Atan(p / q);

             if (p < 0)
             {
                 if (q < 0)
                 {
                     A2 = Math.Abs(A2);
                 }
                 else
                 {
                     A2 =PI - Math.Abs(A2);
                 }
             }
             else
             {
                 if (p > 0)
                 {
                     A2 = PI + Math.Abs(A2);
                 }
                 else
                 {
                     A2 = 2 * PI -Math.Abs(A2);
                 }
             }

             //将弧度制角转化为角度制角
             Ang1 = radian_angle(A1);
             Ang2 = radian_angle(A2);

             textBox33.Text = Convert.ToString(Ang1[0]); textBox34.Text = Convert.ToString(Ang1[1]); textBox35.Text = Convert.ToString(Ang1[2]);
             textBox36.Text = Convert.ToString(Ang2[0]); textBox37.Text = Convert.ToString(Ang2[1]); textBox38.Text = Convert.ToString(Ang2[2]);
             textBox32.Text = Convert.ToString(S);


        }

    /// <summary>
        /// 计算两个经纬度之间的球面距离
        /// </summary>
        /// <param Name="p1">p1的经纬度</param>
        /// <param Name="p2">p2的经纬度</param>
        /// <returns>距离(单位是米)</returns>
        public static double Distance(PointF p1, PointF p2)
        {
            double _nRslt;
            double x1 = p1.X * Math.PI / 180.0;
            double y1 = p1.Y * Math.PI / 180.0;
            double x2 = p2.X * Math.PI / 180.0;
            double y2 = p2.Y * Math.PI / 180.0;
            _nRslt = Math.Sin(y1) * Math.Sin(y2) + Math.Cos(y1) * Math.Cos(y2) * Math.Cos(Math.Abs(x2 - x1));
            if (_nRslt < 1)
            {
                _nRslt = Math.Atan(-_nRslt / Math.Sqrt(-_nRslt * _nRslt + 1)) + 2.0 * Math.Atan(1);
            }
            else
            {
                _nRslt = 0.0;
            }

            _nRslt = _nRslt * DrawArgs.EquatorialRadius;
            return _nRslt;
        }