惯性系下航天器通过拉格朗日系数法进行轨道外推

输入t为航天器在轨道上运行的时间，R0_c、R0_t、V0_c、V0_c为两个不同航天器的位置速度信息，r_c为航天器的半通径，由已知轨道六根数计算：

self.r_c = self.a_c * (1 - self.e_c ** 2) / (1 + self.e_c * np.cos(self.f0_c))

u为引力常数，我这里设为了3.986e5, a_c为半长轴。

    def calculation_spacecraft_status(self, t):
        """
        :param t: 航天器轨道状态信息外推时间
        :return: 航天器轨道状态信息
        """
        def delta_E_equation_t(delta_E):
            # 根据当前时刻的平近点角差确定偏近点角差()
            eq = delta_E + sigma_t * (1 - np.cos(delta_E)) / np.sqrt(self.a_t) - \
                 (1 - self.r_t / self.a_t) * np.sin(delta_E) - np.sqrt(self.u / self.a_t**3) * t
            return eq

        def delta_E_equation_c(delta_E):
            # 根据当前时刻的平近点角差确定偏近点角差()
            eq = delta_E + sigma_c * (1 - np.cos(delta_E)) / np.sqrt(self.a_c) - \
                 (1 - self.r_c / self.a_c) * np.sin(delta_E) - np.sqrt(self.u / self.a_c**3) * t
            return eq

        def Lagrange_factor_c(sigma, delta_E):
            # 拉格朗日系数,外推轨道位置速度
            r = self.a_c + (self.r_c - self.a_c) * np.cos(delta_E) + sigma * np.sqrt(self.a_c) * np.sin(delta_E)
            F = 1 - self.a_c * (1 - np.cos(delta_E)) / self.r_c
            G = self.a_c * sigma * (1 - np.cos(delta_E)) / np.sqrt(self.u) + \
                self.r_c * np.sqrt(self.a_c / self.u) * np.sin(delta_E)
            F_c = -np.sqrt(self.u * self.a_c) * np.sin(delta_E) / (r * self.r_c)
            G_c = 1 - self.a_c * (1 - np.cos(delta_E)) / r

            return F, G, F_c, G_c
        # 偏近点角差

        sigma_t = np.dot(self.R0_t, self.V0_t) / np.sqrt(self.u)
        sigma_c = np.dot(self.R0_c, self.V0_c) / np.sqrt(self.u)
        # 当前时刻处的偏近点角差
        delta_E_t = fsolve(delta_E_equation_t, 0)
        delta_E_c = fsolve(delta_E_equation_c, 0)
        # print(self.R0_t, self.V0_t, self.R0_c, self.V0_c)
        # print(delta_E_t, sigma_c, self.a_c, self.r_c / self.a_c, np.sqrt(self.u / self.a_c ** 3) * t)
        #
        # 航天器位置速度外推(惯性系下)
        F, G, F_t, G_t = self.Lagrange_factor(sigma_t, delta_E_t)
        R_i_t = F * self.R0_t + G * self.V0_t
        V_i_t = F_t * self.R0_t + G_t * self.V0_t

        F, G, F_c, G_c = Lagrange_factor_c(sigma_c, delta_E_c)
        R_i_c = F * self.R0_c + G * self.V0_c
        V_i_c = F_c * self.R0_c + G_t * self.V0_c


        # 惯性系至近心点坐标系的转移矩阵
        C_Jg = self.coordinate_J_to_g()
        # 近心点坐标系下航天器位置
        R_i_t = np.dot(np.linalg.inv(C_Jg), R_i_t)
        V_i_t = np.dot(np.linalg.inv(C_Jg), V_i_t)
        R_i_c = np.dot(np.linalg.inv(C_Jg), R_i_c)
        V_i_c = np.dot(np.linalg.inv(C_Jg), V_i_c)

        return np.array([R_i_t, V_i_t]).ravel(), np.array([R_i_c, V_i_c]).ravel()

我这里是输入两个航天器的状态（位置速度），外推时间t，将状态转为为六根数后计算拉格朗日系数，最后求解在时间t后的航天器状态。

除此之外还可以通过CW方程进行轨道外推，以及龙格库塔算法等，如有需要我将会在后续的帖子中叙述。