GPS广播星历批量计算卫星位置Python

import numpy as np
from math import *
#算法函数
def GPS(year,month,day,hour,minute,second,a0,a1,a2,IodeIssueofData,Crs,deta_n,M0,Cuc,
        e,Cus,a_sqrt,t0e,Cic,W0,Cis,i0,Crc,w,W_,I_):
    #1.计算卫星运行的平均角速度n
    u=3.986004418*e+14;
    n0=pow(u,0.4)*pow( a_sqrt,3);
    n=n0+deta_n;
    #2.计算归化时间tk
    s1=int((1461*(year+4800+ int((month-14)/12 )))/4 );
    s2= int(367*(month-2-int((month-14)/12 )*12)/12 );
    s3=int(3*int((year+4900+int((month-14)/12) )/100)/4 );
    jd=day-32074+s1+s2-s3-0.4+hour/23+minute/1440+second/86400;
    z=int((jd-2342736.4)/7);
    t_=((jd-2342736.4)/7-z)*604800;
    t0c=t_;
    deta_t=a0+a1*(t_-t0c)+a2*(t_-t0c);
    t=t_-deta_t;
    tk=t-t0e;
    #3.观测时刻卫星平近点角Mk的计算
    Mk=M0+n*tk;
    #4.计算偏近点角Ek
    E=1;
    Ek=Mk;
    while E>pow(10,-6):
        E=Ek;
        Ek=Mk+e*sin(degrees(Ek)); 
        E=abs(Ek-E);
    #4.真近点角Vk的计算
    Vk=atan(degrees((sqrt(1-e*e)*sin(degrees(Ek)))/(cos(degrees(Ek))-e)));
    #6.升交距角Fk的计算
    Fk=Vk+w;
    #7.摄动改正项su,sr,si的计算
    su=Cuc*cos(2*degrees(Fk))+Cus*sin(2*degrees(Fk));
    sr=Crc*cos(2*degrees(Fk))+Crs*sin(2*degrees(Fk));
    si=Cic*cos(2*degrees(Fk))+Cis*sin(2*degrees(Fk));
    #8.计算经过改正的升交距角uk、卫星矢径rk和轨道倾角ik
    uk=Fk+su;
    rk=a_sqrt*a_sqrt*(1-e*cos(degrees(Ek)))+sr;
    ik=i0+si+I_*tk;
    #9.计算卫星在轨道平面直角坐标系的坐标
    xk=rk*cos(degrees(uk));
    yk=rk*sin(degrees(uk));
    #10.观测时刻升交点经度Wk的计算
    we=7.29211467*e-4;
    Wk=W0+(W_-we)*tk-we*t0e;
    #11.计算卫星在地心固定坐标系中的直角坐标
    Xk=xk*cos(degrees(Wk))-yk*cos(degrees(ik))*sin(degrees(Wk));
    Yk=xk*sin(degrees(Wk))+yk*cos(degrees(ik))*cos(degrees(Wk));
    Zk=yk*sin(degrees(ik));
    po=[Xk,Yk,Zk]
    return po
#计算文件内所有卫星的位置，自行删除文件头
#打开文件
Gf=open('你的文件路径','r');
#文件行数
line=-1;
position=np.zeros((500,10))
po=np.zeros(3)
class Sa:
    PRN=0
    Year=0
    Month=0
    Day=0
    Hour=0
    Minute=0
    Second=0
    ClockBias=0
    ClockDirft=0
    ClockDirftRate=0
    IodeIssueofData=0
    Crs=0
    Delta_n=0
    M0=0
    Cuc=0
    e=0
    Cus=0
    sqrt_A=0
    toe=0
    Cic=0
    OMEGA=0
    Cis=0
    i0=0
    Crc=0
    omega=0
    OMEGADOT=0
    IDOT=0
    Codes0nL2Channel=0
    GPSWeek=0
    L2_p_data_flag=0
    SV_accuracy=0
    SV_health=0
    TGD=0
    IODC_Issue_of_Data=0
    Transmission_time_of_message=0
    Fit_interval=0
satellite=[Sa() for i in range(500)]
#读取数据
i=0;#卫星数量
while True:  
    str =Gf.readline()  
    if not str:  
        break  
    line=line+1;
    #直接使用动态结构体satellite    
    if line%8==0:
        satellite[i].PRN=int(str[1:3]);
        satellite[i].Year=int(str[4:8]);
        satellite[i].Month=int(str[9:11]);
        satellite[i].Day=int(str[12:14]);
        satellite[i].Hour=int(str[14:17]);
        satellite[i].Minute=int(str[18:20]);
        satellite[i].Second=float(str[21:23]);
        satellite[i].ClockBias=float(str[23:42]);
        satellite[i].ClockDirft=float(str[42:61]);
        satellite[i].ClockDirftRate=float(str[61:80]);      
    if line%8==1:
        satellite[i].IodeIssueofData=float(str[4:23]);
        satellite[i].Crs=float(str[23:42]);
        satellite[i].Delta_n=float(str[42:61]);
        satellite[i].M0=float(str[61:80]);
    
    if line%8==2:
        satellite[i].Cuc=float((str[4:23]));
        satellite[i].e=float((str[23:42]));
        satellite[i].Cus=float((str[42:61]));
        satellite[i].sqrt_A=float((str[61:80]));
    
    if line%8==3:
        satellite[i].toe=float((str[4:23]));
        satellite[i].Cic=float((str[23:42]));
        satellite[i].OMEGA=float((str[42:61]));
        satellite[i].Cis=float((str[61:80]));
    
    if line%8==4:
        satellite[i].i0=float((str[4:23]));
        satellite[i].Crc=float((str[23:42]));
        satellite[i].omega=float((str[42:61]));
        satellite[i].OMEGADOT=float((str[61:80]));
    
    if line%8==5:
        satellite[i].IDOT=float((str[4:23]));
        satellite[i].CodesOnL2Channel=float((str[23:42]));
        satellite[i].GPSWeek=float((str[42:61]));
        satellite[i].L2_p_data_flag=float((str[61:80]));
    
    if line%8==6:
        satellite[i].SV_accuracy=float((str[4:23]));
        satellite[i].SV_health=float((str[23:42]));
        satellite[i].TGD=float((str[42:61]));
        satellite[i].IODC_Issue_of_Data=float((str[61:80]));
   
    if line%8==7:
        satellite[i].Transmission_time_of_message=float((str[4:23]));
        satellite[i].Fit_interval=float((str[23:42]));
   #第一列为卫星PRN号，2:4列为三维坐标，4:10列为时间
        po=GPS(satellite[i].Year,satellite[i].Month,satellite[i].Day,satellite[i].Hour,satellite[i].Minute ,
        satellite[i].Second,satellite[i].ClockBias,satellite[i].ClockDirft,satellite[i].ClockDirftRate, 
        satellite[i].IodeIssueofData,satellite[i].Crs ,      satellite[i].Delta_n,   satellite[i].M0, 
        satellite[i].Cuc,            satellite[i].e,         satellite[i].Cus,       satellite[i].sqrt_A, 
        satellite[i].toe,            satellite[i].Cic,       satellite[i].OMEGA,     satellite[i].Cis,               
        satellite[i].i0, satellite[i].Crc,  satellite[i].omega,satellite[i].OMEGADOT,satellite[i].IDOT);
        #卫星定轨只需要传到IDOT
        position[i,0]=satellite[i].PRN;position[i,1]=po[0];position[i,2]=po[1]; position[i,3]=po[2];
        position[i,4]=satellite[i].Year;position[i,5]=satellite[i].Month;position[i,6]=satellite[i].Day;
        position[i,7]=satellite[i].Hour;position[i,8]=satellite[i].Minute;position[i,9]=satellite[i].Second;
        i=i+1