写一段python代码，根据给定的导航数据文件（RINEX，brdc0010.21m）格式，计算每颗卫星在所给参考历元前后共一个小时内每隔30s卫星在地固坐标系中的位置。

好的，以下是一段用Python实现的代码，可以计算每个卫星在指定时间范围内每隔30秒的位置信息：

import pandas as pd
import numpy as np
from datetime import datetime, timedelta
from math import pi, sin, cos, sqrt

# 读取导航数据文件
data = pd.read_csv('brdc0010.21m', skiprows=8, header=None, delim_whitespace=True)

# 定义常数
GM = 3.986005e14  # 地球引力常数
OMEGA_DOT_E = 7.2921151467e-5  # 地球自转角速度
C = 2.99792458e8  # 光速

# 定义函数，计算卫星在指定时间的位置
def calc_sat_pos(t, sat):
    week = int(sat['week'])
    toe = sat['toe']
    M0 = sat['m0']
    e = sat['ecc']
    A = sat['sqrta'] ** 2
    n0 = sqrt(GM / A ** 3)
    n = n0 + sat['delta_n']
    tk = t.timestamp() - (week * 604800 + toe)
    Mk = M0 + n * tk
    E = Mk
    for i in range(10):
        E_old = E
        E = Mk + e * sin(E_old)
        if abs(E - E_old) < 1e-12:
            break
    v = np.arctan2(sqrt(1 - e ** 2) * sin(E), cos(E) - e)
    phi = v + sat['omega']
    delta_u = sat['cus'] * sin(2 * phi) + sat['cuc'] * cos(2 * phi]
    delta_r = sat['crs'] * sin(2 * phi) + sat['crc'] * cos(2 * phi)]
    delta_i = sat['cis'] * sin(2 * phi) + sat['cic'] * cos(2 * phi)]
    u = phi + delta_u
    r = A * (1 - e * cos(E)) + delta_r
    i = sat['inc'] + delta_i + sat['idot'] * tk
    x = r * cos(u)
    y = r * sin(u)
    Omega = sat['omegadot'] - OMEGA_DOT_E
    X = x * cos(Omega) - y * cos(i) * sin(Omega)
    Y = x * sin(Omega) + y * cos(i) * cos(Omega)
    Z = y * sin(i)
    return np.array([X, Y, Z])

# 定义函数，计算两个向量之间的夹角
def calc_angle(v1, v2):
    return np.arccos(np.dot(v1, v2) / np.sqrt(np.dot(v1, v1) * np.dot(v2, v2)))

# 读取参考历元
ref_epoch = datetime(2021, 1, 1, 0, 0, 0)

# 计算时间范围
start_time = ref_epoch
end_time = ref_epoch + timedelta(hours=1)

# 初始化结果列表
result = []

# 遍历每个卫星
for i in range(0, len(data), 8):
    # 读取卫星信息
    sat = {
        'prn': data[0][i],
        'week': data[1][i],
        'toe': data[2][i],
        'm0': data[3][i],
        'delta_n': data[4][i],
        'ecc': data[5][i],
        'sqrta': data[6][i],
        'omega': data[7][i],
        'inc': data[8][i],
        'omegadot': data[9][i],
        'idot': data[10][i],
        'crc': data[11][i],
        'crs': data[12][i],
        'cuc': data[13][i],
        'cus': data[14][i],
        'cic': data[15][i],
        'cis': data[16][i],
        'toe_toc': data[17][i],
        'af0': data[18][i],
        'af1': data[19][i],
        'af2': data[20][i],
    }
    # 计算卫星在参考历元的位置
    sat_pos_ref = calc_sat_pos(ref_epoch, sat)
    # 遍历时间范围内的每个30秒
    t = start_time
    while t <= end_time:
        # 计算卫星在当前时间的位置
        sat_pos = calc_sat_pos(t, sat)
        # 计算卫星位置变化量
        delta_pos = sat_pos - sat_pos_ref
        # 计算卫星速度向量
        vel = delta_pos / (30 * C)
        # 计算卫星加速度向量
        acc = vel / (30 * C)
        # 计算卫星速度和加速度向量的夹角
        angle = calc_angle(vel, acc)
        # 将结果添加到结果列表中
        result.append({
            'prn': sat['prn'],
            'time': t,
            'pos': sat_pos,
            'vel': vel,
            'acc': acc,
            'angle': angle
        })
        # 增加30秒
        t += timedelta(seconds=30)

# 将结果保存到文件
df = pd.DataFrame(result)
df.to_csv('sat_pos.csv', index=False)

这段代码做了以下几件事情：

1. 读取导航数据文件，提取每个卫星的信息。
2. 定义函数 `calc_sat_pos`，计算指定时间卫星在地固坐标系中的位置。
3. 定义函数 `calc_angle`，计算两个向量之间的夹角。
4. 遍历每个卫星，计算在指定时间范围内每隔30秒卫星在地固坐标系中的位置、速度、加速度和速度、加速度之间的夹角。
5. 将结果保存到文件中。