Commit 391842bd authored by ChenHao's avatar ChenHao

模块一提交

parent 2100f28a
import math
#import matplotlib.pyplot as plt # 作图
def TLE_read(path): # 读取TLE两行轨道数据,以列表形式返回
with open(path, "r") as file: # 打开文件,数据存入列表 data_list
data_list = file.readlines()
num = int(len(data_list) / 3) # 卫星数量,每个卫星的LTE数据为3行
# 初始化列表
name = [0] * num # 卫星名称
serial_num = [0] * num # 卫星编号
t_1 = [0] * num # TLE历时(年份)
t_0 = [0] * num # TLE历时(天数)
delta = [0] * num # 轨道倾角
OMEGA = [0] * num # 升交点赤经
M_0 = [0] * num # 平近点角
turn = [0] * num # 每天环绕地球的圈数
e = [0] * num # 轨道偏心率
e_0 = 0 # 轨道偏心率小数部分(TLE文件中只有小数部分)
e_l = 0 # 轨道偏心率长度
for i in range(num): # 将有效数据分别存入列表,参考 https://blog.csdn.net/haochajin/article/details/83095447
# 参考 https://www.cnblogs.com/rainbow70626/p/6012162.html
name[i] = data_list[i * 3 + 0] # 卫星名称在第 0 行
serial_num[i] = data_list[i * 3 + 1][2:7] # 第一行
t_1[i] = int(data_list[i * 3 + 1][18:20]) # 第一行(float将字符串转化为浮点数)
t_0[i] = float(data_list[i * 3 + 1][20:32]) # 第一行
delta[i] = float(data_list[i * 3 + 2][8:16]) / 180 # 第二行(除 180 表示弧度)
OMEGA[i] = float(data_list[i * 3 + 2][17:25]) / 180 # 第二行
M_0[i] = float(data_list[i * 3 + 2][43:51]) / 180 # 第二行
turn[i] = float(data_list[i * 3 + 2][52:63]) # 第二行
e_0 = float(data_list[i * 3 + 2][26:33]) # 小数部分
e_l = len(str(int(data_list[i * 3 + 2][26:33])))
# 取其整数部分长度,将其化为纯小数(0->1),str 前面有多余的零所以要 int 一下
e[i] = e_0 / 10 ** e_l # 化为0->1的小数
for i in range(num): # 截取有效卫星名称
k = 9 # 至少有9位有效数据
while name[i][k] != " ": # 到空格为止
k += 1
else:
name[i] = name[i][:k]
return num, name, serial_num, t_1, t_0, delta, OMEGA, M_0, turn, e # 返回读取的原始数据
TLE_num, TLE_name, TLE_serial_number, TLE_t_1, TLE_t_0, TLE_delta, TLE_OMEGA, TLE_M_0, TLE_turn, TLE_e = TLE_read(
r"D:\py test\iridium-NEXT.txt")
'''
# 子函数输出无误
for i in range(1):
print(TLE_num,'\n')
print(len(TLE_name),TLE_name,'\n')
print(len(TLE_serial_number),TLE_serial_number,'\n')
print(len(TLE_t_1),TLE_t_1,'\n')
print(len(TLE_t_0),TLE_t_0,'\n')
print(len(TLE_delta),TLE_delta,'\n')
print(len(TLE_OMEGA),TLE_OMEGA,'\n')
print(len(TLE_M_0),TLE_M_0,'\n')
print(len(TLE_turn),TLE_turn,'\n')
print(len(TLE_e),TLE_e,'\n')
plt.plot(TLE_e)
plt.show()
'''
def data_prepare_TI(num,turn): # 对不能直接使用的TLE数据进行转换计算(与时间无关的参数)
G = 6.67428*10**-11 # 万有引力长常量(N·m²/kg²)2006 年国际推荐数值
M = 5.965*10**24 # 地球重量(kg)
# K = #开普勒常数
pai = math.pi # 圆周率
angle_v = [0]*num # 卫星角速度
a = [0]*num # 轨道半长轴
for i in range(num):
angle_v[i] = 2*pai / (86400/turn[i]) # ω = 2pai / T
a[i] = ( (G*M / angle_v[i]**2) )**(1/3) # a = (GM / angle_v ** 2)^(1/3)
return angle_v,a
TLE_angle_v,TLE_a = data_prepare_TI(TLE_num,TLE_turn)
'''
# 子函数输出无误
for i in range(1):
print(len(TLE_angle_v),TLE_angle_v,'\n')
print(len(TLE_a),TLE_a,'\n')
for j in range(len(TLE_a)): # 检查半长轴
if (TLE_a[j] < 6371000):
print("TLE_a[j]半长轴小于地球半径!")
TLE_b = [0]*TLE_num # 半短轴
for i in range(TLE_num):
TLE_b[i] = TLE_a[i] * (1-TLE_e[i]**2)**(1/2)
print(TLE_b)
plt.plot(TLE_a)
plt.plot(TLE_b)
plt.show()
'''
def elapsed_time(num, time_start, t_0, turn, ST): # 运行时间,根据仿真开始时间,得到所有计算平近点角的时间间隔
T = 0
for i in range(num):
T += 86400 / turn[i] # 卫星周期
T_ave = T / num # 平均卫星周期
TS_num = int(T_ave // ST) # 时隙数量,时隙长为 60s
tm_yday = 0 # 初始化
t_tsx = [[0] * TS_num for i in range(num)] # 到每个时隙的间隔
month = [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # 2020 年每月天数
for i in range(time_start[1] - 1):
tm_yday += month[i]
tm_yday = tm_yday + time_start[2] # 2020 年第几天
tm_hour = time_start[3]
tm_min = time_start[4]
tm_sec = time_start[5]
decimals = (tm_hour * 3600 + tm_min * 60 + tm_sec) / 86400 # 将时分秒化为天数
start_time = tm_yday + decimals
for i in range(num): # 每颗卫星
for j in range(TS_num): # 每个时隙
t_tsx[i][j] = (start_time - t_0[i]) * 86400 + j * 60
return TS_num, t_tsx
ST = 60
TS_start = [2020, 12, 25, 8, 0, 0]
TS_num, TLEt_tsx = elapsed_time(TLE_num, TS_start, TLE_t_0, TLE_turn, ST)
'''
# 子函数输出无误
print(TS_num)
for i in range(1):
print(TLEt_tsx[0])
print(len(TLEt_tsx),len(TLEt_tsx[0]))
for i in range(TLE_num):
plt.plot(TLEt_tsx[i])
plt.show()
'''
def solve_E(M, e): # 解方程,M = E - e*sin(E)
pai = math.pi
E1 = 0
E2 = 2 * pai
E = pai # 二分查找法
while (abs(E - M - e * math.sin(E)) > 0.001): # f(E) = E - M - esinE 函数单增,误差小于0.001
if (E - M - e * math.sin(E) > 0):
E2 = E
E = (E1 + E2) / 2
elif (E - M - e * math.sin(E) < 0):
E1 = E
E = (E1 + E2) / 2
E = (E1 + E2) / 2
return E
# E = solve_E(1.5472,0.1937) # 弧度值
# print(E)
def data_prepare_TD(num, M_0, e, angle_v, a, TS_num, t_tsx): # 对不能直接使用的TLE数据进行转换计算(与时间有关的参数)
pai = math.pi # 圆周率
M = [[0] * TS_num for i in range(num)] # 所有时刻的平近点角,弧度值,与时间是线性关系,周期为 2π
E = [[0] * TS_num for i in range(num)] # 所有时刻的偏近点角,弧度值
r = [[0] * TS_num for i in range(num)] # 所有时刻卫星离地心的距离
f = [[0] * TS_num for i in range(num)] # 所有时刻的真近点角,弧度值
for i in range(num): # 每颗卫星
for j in range(TS_num): # 每个时隙
M[i][j] = (M_0[i] + angle_v[i] * t_tsx[i][j]) % (2 * pai) # M = Mo + angle_v * (t - to)
E[i][j] = solve_E(M[i][j], e[i])
# print(1 - e[i] * math.cos(E[i][j]))
r[i][j] = a[i] * (1 - e[i] * math.cos(E[i][j])) # r = a *(1 - e * cosE)
f[i][j] = math.atan(((1 - e[i] ** 2) ** (1 / 2) * math.sin(E[i][j])) / (math.cos(E[i][j]) - e[i]))
# f = arctan( (1-e^2)^(1/2)*sinE / (cosE-e) )
return M, E, r, f
TLE_M, TLE_E, TLE_r, TLE_f = data_prepare_TD(TLE_num, TLE_M_0, TLE_e, TLE_angle_v, TLE_a, TS_num, TLEt_tsx)
'''
for i in range(1):
# print(len(TLE_M),len(TLE_M[0]),'\n')
# print(TLE_M[0],'\n')
# print(len(TLE_E),len(TLE_E[0]),'\n')
# print(TLE_E[0],'\n')
# print(len(TLE_r),len(TLE_r[0]),'\n')
# print(TLE_r[0],'\n')
# print(len(TLE_f),len(TLE_f[0]),'\n')
# print(TLE_f[1],'\n')
# for i in range(3): # 输出无误
# plt.plot(TLE_M[i])
# plt.grid(True) # 增加格点
# plt.show()
# for i in range(3): # 输出无误
# plt.plot(TLE_E[i])
# plt.grid(True) # 增加格点
# plt.show()
X = [6371000]*TLE_num
x0 = list(range(0,TLE_num))
plt.plot(TLE_a)
plt.plot(X)
r = [[0]* TLE_num for i in range(TS_num)]
for i in range(TLE_num):
for j in range(TS_num):
r[j][i] = TLE_r[i][j]
for i in range(1): # 输出不应该比半径还小
#plt.plot(r[i])
plt.scatter(x0,r[i],marker='o')
for i in range(1): # 输出不应该比半径还小
plt.plot(TLE_r[i])
plt.grid(True) # 增加格点
plt.show()
# for i in range(3): # 输出无误
# plt.plot(TLE_f[i])
# plt.grid(True) # 增加格点
# plt.show()
'''
def coordinates(num,delta,OMEGA,TS_num,r,f): # 输入相关参数,输出直角坐标
x = [[0]* TS_num for i in range(num)] # x坐标
y = [[0]* TS_num for i in range(num)] # y坐标
z = [[0]* TS_num for i in range(num)] # z坐标
for i in range(num): # 计算并储存卫星的直角坐标
for j in range(TS_num):
x[i][j] = r[i][j] * math.cos(f[i][j]) * math.cos(OMEGA[i]) \
- r[i][j] * math.sin(f[i][j]) * math.sin(OMEGA[i]) * math.cos(delta[i])
y[i][j] = r[i][j] * math.cos(f[i][j]) * math.sin(OMEGA[i]) \
+ r[i][j] * math.cos(f[i][j]) * math.cos(OMEGA[i]) * math.cos(delta[i])
z[i][j] = r[i][j] * math.sin(f[i][j]) * math.sin(delta[i])
return x,y,z # 返回直角坐标到列表
x,y,z = coordinates(TLE_num,TLE_delta,TLE_OMEGA,TS_num,TLE_r,TLE_f)
'''
for i in range(1):
plt.plot(x[0])
plt.show()
plt.plot(y[0])
plt.show()
plt.plot(z[0])
plt.show()
'''
def distance(num, TS_num, x, y, z): # 计算并储存任意两颗卫星之间的距离
dis_mat = [[[0] * num for i in range(num)] for j in range(TS_num)] # 以三维列表储存任意两颗卫星之间的距离,第一层为不同时隙
for i in range(TS_num):
for j in range(num):
for k in range(j + 1, num): # 上三角矩阵
dis_mat[i][j][k] = \
((x[j][i] - x[k][i]) ** 2 + \
(y[j][i] - y[k][i]) ** 2 + \
(z[j][i] - z[k][i]) ** 2) ** (1 / 2)
return dis_mat
dis_mat = distance(TLE_num, TS_num, x, y, z)
'''
for i in range(1): # 没想好怎么验证准确性
print(type(dis_mat))
print(len(dis_mat))
print(len(dis_mat[0]))
print(len(dis_mat[0][0]))
x0 = list(range(0,75))
for i in range(0,75):
plt.scatter(x0,dis_mat[0][i],marker='o') # 以散点图的形式展示这些距离
plt.show()
'''
def judge_visual(num, r, dis_mat): # 判断并储存任意两颗卫星之间的可视条件
Re = 6378124 # 地球半径
h = 50000 # 电离层高度,铱星Next的星间链路采用23GHz微波频率,假设在电离层进行通信
gamma = 2.498 # 铱星波束扫描角度范围(弧度),根据其轨道高度与地面覆盖范围(约数)得出
visual_mat = [[[0] * num for i in range(num)] for j in range(TS_num)] # 以三维列表储通断矩阵
R = Re + h
ga = gamma / 2
for i in range(TS_num):
for j in range(num):
for k in range(j + 1, num): # 上三角矩阵
dis_min = r[j][i] * math.cos(ga) + r[k][i] * math.cos(ga)
# dis_max = (r[j][i]**2-R**2)**(1/2) + (r[k][i]**2-R**2)**(1/2)
dis_max = abs((r[j][i] ** 2 - R ** 2)) ** (1 / 2) + abs((r[k][i] ** 2 - R ** 2)) ** (1 / 2)
if (dis_mat[i][j][k] >= dis_min and dis_mat[i][j][k] <= dis_max):
visual_mat[i][j][k] = 1
else:
visual_mat[i][j][k] = 0
return visual_mat
visual_mat = judge_visual(TLE_num, TLE_r, dis_mat)
# for i in range(len(visual_mat[0])):
# print(visual_mat[0][i][:40])
def self_mat(num): # 相邻时隙间的同一颗卫星是联通的
matrix = [[0] * num for i in range(num)] # 以二维储时通断矩阵
for i in range(num):
matrix[i][i] = 1 # 相邻时隙间卫星连通
return matrix
mat_self = self_mat(TLE_num)
'''
for i in range(1):
print(len(mat_self))
print(len(mat_self[0]))
for i in range(len(mat_self)):
print(i,mat_self[i])
'''
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment