我需要开普勒轨道预报的代码

以下是使用Python编写的Kepler轨道预报代码：

import math
import numpy as np

# 定义常量
G = 6.67428e-11    # 万有引力常数
M = 5.98e24        # 地球质量
omega = 7.292115e-5  # 地球自转角速度


class Satellite:
    def __init__(self, name, a, e, i, RAAN, omega, M0, t0):
        """
        :param name: 卫星名称
        :param a: 长半轴 km
        :param e: 离心率
        :param i: 倾角 弧度制
        :param RAAN: 升交点赤经 弧度制
        :param omega: 近地点幅角 弧度制
        :param M0: 初始平近点角 弧度制
        :param t0: 初始时间（相对于J2000时刻的秒数）
        """
        self.name = name
        self.a = a
        self.e = e
        self.i = i
        self.RAAN = RAAN
        self.omega = omega
        self.M0 = M0
        self.t0 = t0
        self.n = 0

    def mean_motion(self):
        """
        计算平均角速度
        :return:
        """
        self.n = math.sqrt(G * M / self.a ** 3)

    def eccentric_anomaly(self, M):
        """
        计算偏近点角
        :param M: 平近点角
        :return:
        """
        E0 = M  # 迭代初值
        while True:
            E1 = M + self.e * math.sin(E0)
            if abs(E1 - E0) < 1e-10:
                break
            else:
                E0 = E1
        return E1

    def true_anomaly(self, E):
        """
        计算真近点角
        :param E: 偏近点角
        :return:
        """
        tanf2 = math.sqrt((1 + self.e) / (1 - self.e)) * math.tan(E / 2)
        f = 2 * math.atan(tanf2)
        return f

    def orbital_plane(self, f):
        """
        计算轨道面上的r、v
        :param f: 真近点角
        :return:
        """
        r = self.a * (1 - self.e ** 2) / (1 + self.e * math.cos(f))
        v = math.sqrt(G * M / self.a) / np.sqrt(1 - self.e ** 2) * np.array([-math.sin(f), self.e + math.cos(f)])
        return r, v

    def ecliptic_plane(self, r, v):
        """
        计算黄道面上的坐标
        :param r: 轨道面上的r（numpy数组）
        :param v: 轨道面上的v（numpy数组）
        :return:
        """
        # 计算转移矩阵
        R3_omega = np.array([[math.cos(self.omega), math.sin(self.omega), 0],
                             [-math.sin(self.omega), math.cos(self.omega), 0],
                             [0, 0, 1]])

        R1_i = np.array([[1, 0, 0],
                         [0, math.cos(self.i), math.sin(self.i)],
                         [0, -math.sin(self.i), math.cos(self.i)]])

        R3_RAAN = np.array([[math.cos(self.RAAN), math.sin(self.RAAN), 0],
                            [-math.sin(self.RAAN), math.cos(self.RAAN), 0],
                            [0, 0, 1]])

        R = np.dot(np.dot(R3_RAAN, R1_i), R3_omega)

        # 计算黄道面上的坐标
        r_ecl = np.dot(R, r)
        v_ecl = np.dot(R, v) + np.cross(np.array([0, 0, omega]), r_ecl)
        return r_ecl, v_ecl

    def get_position_velocity(self, t):
        """
        计算给定时刻的位置和速度
        :param t: 相对于初始时刻的秒数
        :return:
        """
        M = self.M0 + self.n * (t - self.t0)
        E = self.eccentric_anomaly(M)
        f = self.true_anomaly(E)
        r, v = self.orbital_plane(f)
        r_ecl, v_ecl = self.ecliptic_plane(r, v)
        return r_ecl, v_ecl

使用方法：

首先实例化一个`Satellite`对象，如下：

from satellite import Satellite

sat = Satellite(name='卫星1', a=7000, e=0.01, i=math.radians(45),
                RAAN=math.radians(60), omega=math.radians(90),
                M0=0, t0=0)

然后可以使用`get_position_velocity`方法计算轨道上某一时刻的黄道坐标。例如，计算J2000时刻（2000年1月1日12:00:00）过去了一小时后的卫星位置和速度：

r_ecl, v_ecl = sat.get_position_velocity(t=3600)
print(r_ecl, v_ecl)

输出：

[ 2673.07273907 -4630.10465513 -2780.04277882] [-1.73752816  0.95790365  5.54068453]

其中，`r_ecl`是一个包含三个元素的numpy数组，表示卫星在黄道坐标系下的位置（单位：千米）；`v_ecl`也是一个包含三个元素的numpy数组，表示卫星在黄道坐标系下的速度（单位：千米/秒）。