两点经纬度计算方位角

Sub Sample_Point2Azimuth()
    
    Dim lon_A As Double, lat_A As Double, lon_B As Double, lat_B As Double
    
    FinalRow_B = Sheets("程序示例").Cells(Rows.Count, 2).End(xlUp).Row
    Sheets("程序示例").Range("D2", "D" & FinalRow_B).ClearContents
    
    For i = 2 To FinalRow_B
        lon_A = Sheets("程序示例").Cells(1, "G").Value
        lat_A = Sheets("程序示例").Cells(1, "H").Value
        lon_B = Sheets("程序示例").Cells(i, "B").Value
        lat_B = Sheets("程序示例").Cells(i, "C").Value
        Sheets("程序示例").Cells(i, "D").Value = Point2Azimuth(lon_A, lat_A, lon_B, lat_B)
    Next i
    
End Sub
Function Point2Azimuth(lon_A As Double, lat_A As Double, lon_B As Double, lat_B As Double)
    
    ‘其中A为原始点，B为目标点。
    Dim PI As Double
    PI = 3.14159265358979
    If lon_A > 0 And lon_A < 360 And lat_A > 0 And lat_A < 360 Then
        If lon_B > 0 And lon_B < 360 And lat_B > 0 And lat_B < 360 Then
        ‘判断目标点在原始点的第几象限，或者正东南西北
            If lat_B > lat_A Then
                If lon_B > lon_A Then
                    ‘东北方
                        Point2Azimuth = 0
                        lat_A = lat_A * PI / 180
                        lon_A = lon_A * PI / 180
                        lat_B = lat_B * PI / 180
                        lon_B = lon_B * PI / 180
                        
                        Point2Azimuth = Sin(lat_A) * Sin(lat_B) + Cos(lat_A) * Cos(lat_B) * Cos(lon_B - lon_A)
                        Point2Azimuth = Sqr(1 - Point2Azimuth * Point2Azimuth)
                        Point2Azimuth = Cos(lat_B) * Sin(lon_B - lon_A) / Point2Azimuth
                        Point2Azimuth = Application.Asin(Point2Azimuth) * 180 / PI
                ElseIf lon_B < lon_A Then
                    ‘西北方
                        Point2Azimuth = 0
                        lat_A = lat_A * PI / 180
                        lon_A = lon_A * PI / 180
                        lat_B = lat_B * PI / 180
                        lon_B = lon_B * PI / 180
                        
                        Point2Azimuth = Sin(lat_A) * Sin(lat_B) + Cos(lat_A) * Cos(lat_B) * Cos(lon_B - lon_A)
                        Point2Azimuth = Sqr(1 - Point2Azimuth * Point2Azimuth)
                        Point2Azimuth = Cos(lat_B) * Sin(lon_B - lon_A) / Point2Azimuth
                        Point2Azimuth = Application.Asin(Point2Azimuth) * 180 / PI
                        Point2Azimuth = Point2Azimuth + 360
                ElseIf lon_B = lon_A Then
                    ‘正北方
                    Point2Azimuth = 0
                Else
                    ‘Nothing
                End If
            ElseIf lat_B < lat_A Then
                If lon_B > lon_A Then
                    ‘东南方
                        Point2Azimuth = 0
                        lat_A = lat_A * PI / 180
                        lon_A = lon_A * PI / 180
                        lat_B = lat_B * PI / 180
                        lon_B = lon_B * PI / 180
                        
                        Point2Azimuth = Sin(lat_A) * Sin(lat_B) + Cos(lat_A) * Cos(lat_B) * Cos(lon_B - lon_A)
                        Point2Azimuth = Sqr(1 - Point2Azimuth * Point2Azimuth)
                        Point2Azimuth = Cos(lat_B) * Sin(lon_B - lon_A) / Point2Azimuth
                        Point2Azimuth = Application.Asin(Point2Azimuth) * 180 / PI
                        Point2Azimuth = Point2Azimuth + 90
                ElseIf lon_B < lon_A Then
                    ‘西南方
                        Point2Azimuth = 0
                        lat_A = lat_A * PI / 180
                        lon_A = lon_A * PI / 180
                        lat_B = lat_B * PI / 180
                        lon_B = lon_B * PI / 180
                        
                        Point2Azimuth = Sin(lat_A) * Sin(lat_B) + Cos(lat_A) * Cos(lat_B) * Cos(lon_B - lon_A)
                        Point2Azimuth = Sqr(1 - Point2Azimuth * Point2Azimuth)
                        Point2Azimuth = Cos(lat_B) * Sin(lon_B - lon_A) / Point2Azimuth
                        Point2Azimuth = Application.Asin(Point2Azimuth) * 180 / PI
                        Point2Azimuth = Point2Azimuth + 270
                ElseIf lon_B = lon_A Then
                    ‘正南方
                    Point2Azimuth = 180
                Else
                    ‘Nothing
                End If
            ElseIf lat_B = lat_A Then
                If lon_B > lon_A Then
                    ‘正东方
                    Point2Azimuth = 90
                ElseIf lon_B < lon_A Then
                    ‘正西方
                    Point2Azimuth = 270
                ElseIf lon_B = lon_A Then
                    ‘原始点
                    Point2Azimuth = 361
                Else
                    ‘Nothing
                End If
            End If
        
        Else
            ‘目标点B经纬度无效
            Point2Azimuth = -2
        End If
    Else
        ‘原始点A经纬度无效
        Point2Azimuth = -1
                
    End If
       
End Function