starlink卫星轨道预报

SpaceX 公司发射starlink卫星243颗星链卫星，卫星会很大程度影响天文观测。但是有爱好者又想观测到starlink卫星过境，所以我写了以下代码方便大家使用。

需要准备python的轮子包括：urllib，sgp4，spiceypy。

首先从celestrak下载starlink两行报，并将两行报保存在列表嵌套的字典中（不知道卫星两行报的自行百度）。

然后指定卫星名称，读取两行报文件进行预报（或者遍历全部,遍历的话需要并行处理。）。但是由于python模块SGP4输出位置位于TEME参考架，需要转换至J2000进而转换至WGS84坐标系（经度纬度坐标系），就得到卫星的星下点坐标。如果要计算某一位置的starlink观测的方位角和天顶角需要计算卫星与此处的相对位置关系。（这一节先计算到卫星星下点，有需要的可以留言.... frame.py文件为坐标转换文件。见后续blog说明）

# -*- coding: utf-8 -*-
"""
Created on Sun Feb 23 19:54:02 2020

@author: panda
"""
import urllib.request
import os
import platform
from datetime import *
from sgp4.earth_gravity import wgs72,wgs84
from sgp4.io import twoline2rv
import frame
import spiceypy as spice
from math import atan2, cos, pi, sin, sqrt, tan
import numpy as np
def main():

    #星历文件
    f_KERNELS='/home/panda/SPICE/gen/'
    spice.furnsh(f_KERNELS + 'pck/pck00010.tpc')
    spice.furnsh(f_KERNELS+'lsk/naif0011.tls')
    time_intev=60*60#时间间隔,单位：s
    loop=False#True#loop satellite or not（True=会循环所有starlink卫星）
    #star_name='STARLINK-53'#star_link_name
    rad2deg=180./pi
    file_nut80 = 'nut80.dat'地球旋转矩阵
    url='http://www.celestrak.com/NORAD/elements/starlink.txt'
    tle_path='tle'
    result_path='result'
    today=datetime.utcnow().strftime('%Y%m%d')#time utc
    sys_p = platform.system()
    if sys_p == "Windows":
        print("OS is Windows!!!")
        path_tle=os.getcwd()+'\\'+tle_path+'\\'
        path_result=os.getcwd()+'\\'+result_path+'\\'
    elif sys_p == "Linux":
        print("OS is Linux!!!")
        path_tle=os.getcwd()+'/'+tle_path+'/'
        path_result=os.getcwd()+'/'+result_path+'/'+today+'/'
    else:
        print(sys_p)
        path_tle=os.getcwd()+'/'+tle_path+'/'
        path_result=os.getcwd()+'/'+result_path+'/'+today+'/'
    
    print(today)
    data = urllib.request.urlopen(url).read()
    save_tle(data,today,path_tle) 
    tles=split_tle((data))
    try:
        if not os.path.exists(path_result):
            os.makedirs(path_result)
    except IOError as e:
        print(e)
    for tle in tles:
        #print(tle['name'])
        #if tle['name']==star_name:
        new_list = [ x for x in tle['line1'].split(' ') if x != '' ]
        time_y=float(new_list[3])
        year_tle=str(int(time_y/1000)+2000)
        day_lte=(time_y-(int(time_y/1000)*1000))
        utc=year_tle+'-01-01T00:00:00'
        et=spice.spiceypy.utc2et(utc)+(day_lte-1)*24.*60.*60.   
        end_et=et+24.*60.*60.*15
        file_name=path_result+tle['name']+'.txt'
        print(tle['name'])
        with open(file_name, 'w+') as fh:
            while et<= end_et:
                utc=spice.spiceypy.et2utc(et,'ISOC',0,24)

                #得到卫星地面星下点，的经纬度
                Earth2Sat, v_Earth2Sat,lat1,lon1,alt1=satellite_j2000(utc,et,tle['name'],tle['line1'],tle['line2'])
                #print(type(utc),type(lat1),type(lon1))
                fh.write('%s %-8.2f %-8.2f\n'%(str(utc),lat1,lon1))
                et = et + time_intev

    print(len(tles))
   
    spice.kclear()
def save_tle(data,today,path):
    
    #today=datetime.utcnow().strftime('%Y%m%d')#time utc
    file_name=path+today+'.txt'# filename
    try:
        if not os.path.exists(path):
            os.makedirs(path)
        if not os.path.exists(file_name):#if file not exist load file from url
            #url='http://www.celestrak.com/NORAD/elements/starlink.txt'
            #data = urllib.request.urlopen(url).read()
            with open(file_name, 'w+b') as fh:
                fh.write(data)
        else:
            fsize = os.path.getsize(file_name)
            if fsize==0:
                print('reloading:',file_name)
                #url='http://www.celestrak.com/NORAD/elements/starlink.txt'
                #data = urllib.request.urlopen(url).read()
                with open(file_name, 'w+b') as fh:
                    fh.write(data)      
    except IOError as e:
        print(e)
def split_tle(data):
    tle={}
    list_data=data.decode('utf-8').split('\n')
    tles=[]
    for i in range(len(list_data)-1):
        if list_data[i][0]=='1' or list_data[i][0]=='2':
            x=1
        else:
            #print(len(list_data[i]))
            tle['name']=list_data[i].strip() 
            tle['line1']=list_data[i+1].strip() 
            tle['line2']=list_data[i+2].strip() 
            tles.append(tle)
            tle={}
    return(tles)

def load_tle(url):
    mode = re.compile(r'\d+')
    data = urllib.request.urlopen(url).read()#
    data = data.decode('UTF-8')
    tle_lines=data.splitlines()
    for i in range(len(tle_lines)):
        line=tle_lines[i]
        if op.eq(line.strip(),celestrak_flag):
            TlE_name=(line)
            lines1=(tle_lines[i+1])
            lines2=(tle_lines[i+2])
    return(TlE_name,lines1,lines2)
def satellite_j2000(utc,et,TlE_name,lines1,lines2):
    yr=int(utc[0:4])          
    mon=int(utc[5:7])           
    day=int(utc[8:10])           
    hr=int(utc[11:13])         
    mins=int(utc[14:16])           
    sec=int(utc[17:19])  
    satellite = twoline2rv(lines1, lines2, wgs84)
    reci1,veci1= satellite.propagate(yr, mon, day, hr, mins, sec)
    jdtt, jdttfrac = jday(yr, mon, day, hr, mins, sec)
    ttt = (jdtt - 2451545.0) / 36525.0
    opt = '80'
    order = 106
    eqeterms = 2
    ateme = [0, 0, 0]
    file_nut80 = 'nut80.dat'
    reci, veci, aeci = frame.teme2eci(reci1, veci1, ateme, ttt, order, eqeterms, opt, file_nut80)
    lat1,lon1,alt1 = calculate(reci,et)
    return(reci,veci,lat1,lon1,alt1)
def calculate(options,et):
    x = options[0]#*1000.0
    y = options[1]#*1000.0
    z = options[2]#*1000.0
    xtrn=spice.spiceypy.pxform('J2000','IAU_EARTH',et)
    jstate=np.dot(options,np.transpose(xtrn))
    raddii=spice.spiceypy.bodvar(399,"RADII", 3)
    re=raddii[0]
    rp=raddii[2]
    flat=(re-rp)/re
    lon,lat,alt=spice.spiceypy.recpgr('EARTH',jstate,re,flat)
    lat = ((lat*180.)/pi)
    lon = ((lon*180.)/pi)
    return (lat, lon, alt)
def jday(yr, mon, day, hr, mins, sec):
    jd = 367.0 * yr - int( (7 * (yr +int( (mon + 9) / 12.0) ) ) * 0.25 ) + int( 275 * mon / 9.0 )+ day + 1721013.5
    jdfrac = (sec + mins * 60.0 + hr *3600.0) / 86400.0
    if jdfrac > 1.0:
        jd = jd + int(jdfrac)
        jdfrac = jdfrac - int(jdfrac)
    return(jd,jdfrac)
if __name__ == "__main__":
    main()