精密星历卫星钟差插值程序

     本程序采用插值算法，来完成对精密星历钟差文件（.clk)的加密工作。如可以将30S钟差间隔的钟差文件加密到15s、10s、2s、1s等。 程序中写有四种经典的插值算法，分别是Lagrange插值算法、Newton插值算法、三次样条插值算法、分段线性插值算法等，这几种算法的正确性以得到过验证。

    首先导入精密星历钟差文件，根据导入的钟差文件进行判断，将其信息进行整理并绘制可视化图表以显示24H内每颗卫星的钟差变化趋势。

    导入文件之后，可以根据用户自身的插值需求将钟差文件加密到更小的钟差间隔，如导入的精密星历卫星钟差文件为5min钟差间隔，用户可以将其加密到30s、15s、10s以及更小的钟差间隔，并根据结果生成相应的clk文件。

插值算法：

Lagrange插值：

  Public Function Lagrange(X() As Double, Y() As Double, xx As Double)
        Dim yy As Double = 0           '最后返回的多项式值
        Dim ll(X.Length - 1) As Double '中间变量

        For i = 0 To X.Length - 1
            Dim l1 As Double = 1 '多项式系数分子
            Dim l2 As Double = 1 '多项式系数分母

            For j = 0 To X.Length - 1
                If i <> j Then
                    l1 = l1 * (xx - X(j))
                    l2 = l2 * (X(i) - X(j))
                Else
                    Continue For
                End If
            Next
            ll(i) = l1 / l2
        Next

        For i = 0 To X.Length - 1
            yy = yy + ll(i) * Y(i)
        Next
        Return yy
    End Function

Newton插值：

  Public Function Newton(X() As Double, Y() As Double, xx As Double)
        Dim yy As Double = 0                          '最后返回的多项式值
        Dim n As Integer = UBound(X)                  '节点数
        Dim A(n, n) As Double                         '差商表
        Dim dx(n) As Double                           '牛顿多项式中（x-x0）项

        '计算差商表
        For i = 0 To n
            A(i, 0) = Y(i)
        Next
        For i = 1 To n
            For j = 1 To i
                A(i, j) = (A(i, j - 1) - A(i - 1, j - 1)) / (X(i) - X(i - j))
            Next
        Next

        '计算II（x-xi)
        dx(0) = 1
        For i = 1 To n
            dx(i) = dx(i - 1) * (xx - X(i - 1))
        Next

        '计算牛顿多项式
        For i = 0 To n
            yy = yy + A(i, i) * dx(i)
        Next
        Return yy
    End Function

三次样条插值：

 Public Function CubicSpline(X As Double(), Y As Double(), xx() As Double)
        Dim yy(xx.Length - 1) As Double            '函数返回的对应xx 的插值
        Dim n As Integer = UBound(X)
        Dim S2(n, 0), m(n, 0), AA(n, n) As Double  '用于矩阵运算,即AA*M=S2
        Dim h(n - 1) As Double                     '计算步长，即X(i+1)-X(i)
        Dim dy(n - 1) As Double                    'y的微分

        For i = 0 To n - 1
            h(i) = X(i + 1) - X(i)
            dy(i) = (Y(i + 1) - Y(i)) / h(i)
        Next

        '自然边界型三次样条，二阶导数两个端点值为0
        S2(0, 0) = 0 : S2(n, 0) = 0
        For i = 1 To n - 1
            S2(i, 0) = 6.0 * (dy(i) - dy(i - 1))
        Next

        '确定系数阵AA
        AA(0, 0) = 1
        AA(n, n) = 1
        For i = 1 To n - 1
            AA(i, i - 1) = h(i - 1)
            AA(i, i) = 2 * (h(i - 1) + h(i))
            AA(i, i + 1) = h(i)
        Next

        '求二次微分阵m
        m = Multiply(Inv(AA), S2)

        '求多项式
        Dim a, b, c, d, dx As Double
        For j = 0 To xx.Length - 1
            For i = 0 To n
                If (xx(j) > X(i) Or xx(j) = X(i)) And (xx(j) < X(i + 1) Or xx(j) = X(i + 1)) Then
                    a = Y(i)
                    b = (Y(i + 1) - Y(i)) / h(i) - h(i) / 2.0 * m(i, 0) - h(i) * (m(i + 1, 0) - m(i, 0)) / 6
                    c = m(i, 0) / 2
                    d = (m(i + 1, 0) - m(i, 0)) / (6.0 * h(i))
                    dx = xx(j) - X(i)
                    yy(j) = a + b * dx + c * dx ^ 2 + d * dx ^ 3
                    Exit For
                Else
                    Continue For
                End If
            Next
        Next
        Return yy
    End Function

线性插值：

 Public Function PiecewiseLinearity(X As Double(), Y As Double(), xx As Double)
        Dim yy As Double                  '函数返回的对应xx 的插值
        For i = 1 To X.Length - 1
            If xx >= X(i - 1) And xx <= X(i) Then
                yy = (xx - X(i)) / (X(i - 1) - X(i)) * Y(i - 1) + (xx - X(i - 1)) / (X(i) - X(i - 1)) * Y(i)
                Exit For
            Else
                Continue For
            End If
        Next
        Return yy
    End Function

源程序与exe文件下载地址：精密星历卫星钟差插值程序