我创建了一个函数来生成和传播卫星轨道。现在,我想把我所有的数据保存在一个.dat文件中。我不知道我需要多少个for循环,也不知道该怎么做。在
我希望时间,纬度,经度,和高度都在一条线上,为每个传播步骤。在
数据代码:
myOrbitJ2000Time = [1085.0, 2170.0, 3255.0, 4340.1, 5425.1]
lat = [48.5, 26.5, -28.8, -48.1, 0.1]
lon = [-144.1, -50.4, -1.6, 91.5, 152.8]
alt = [264779.5, 262446.1, 319661.8, 355717.3, 306129.0]
所需的.dat文件输出:
J2000时间,纬度,经度,高度
^{pr2}$完整代码:
import orbital
from orbital import earth, KeplerianElements, plot
import matplotlib.pyplot as plt
import numpy as np
from astropy import time
from astropy.time import TimeDelta, Time
from astropy import units as u
from astropy import coordinates as coord
#def orbitPropandcoordTrans(orbitPineapple_J2000time, _ecc, _inc, _raan, _arg_pe, _meananom, meanMotion):
def orbitPropandcoordTrans(propNum, orbitPineapple_J2000time, _ecc, _inc, _raan, _arg_pe, _meananom, meanMotion):
'''
Create original orbit and run for 100 propagations (in total one whole orbit)
in order to get xyz and time for each propagation step.
The end goal is to plot the lat, lon, & alt data to see if it matches ISS groundtrace.
'''
import orbital
from orbital import earth, KeplerianElements, plot
import matplotlib.pyplot as plt
import numpy as np
from astropy import time
from astropy.time import TimeDelta, Time
from astropy import units as u
from astropy import coordinates as coord
'Calculate Avg. Period from Mean Motion'
_avgPeriod = 86400 / meanMotion
#print('_avgPeriod', _avgPeriod)
######
#propNum = int(_avgPeriod)
'Generate Orbit'
#orbitPineapple = KeplerianElements.with_period(5560, body=earth, e=0.0004525, i=(np.deg2rad(51.6414)), raan=(np.deg2rad(247.1662)))
orbitPineapple = KeplerianElements.with_period(_avgPeriod, body=earth, e=_ecc, i=(np.deg2rad(_inc)), raan=(np.deg2rad(_raan)), arg_pe=(np.deg2rad(_arg_pe)), M0=(np.deg2rad(_meananom))) #ref_epoch=
#plot(orbitPineapple)
#plt.show()
#print('')
#print('')
'Propagate Orbit and retrieve xyz'
myOrbitX = [] #X Coordinate for propagated orbit step
myOrbitY = [] #Y Coordinate for propagated orbit step
myOrbitZ = [] #Z Coordinate for propagated orbit step
myOrbitTime = [] #Time for each propagated orbit step
myOrbitJ2000Time = [] #J2000 times
#propNum = 100 #Number of propagations and Mean Anomaly size (one orbit 2pi/propNum)
for i in range(propNum):
orbitPineapple.propagate_anomaly_by(M=(2.0*np.pi/propNum)) #Propagate the orbit by the Mean Anomaly
myOrbitX.append(orbitPineapple.r.x) #x vals
myOrbitY.append(orbitPineapple.r.y) #y vals
myOrbitZ.append(orbitPineapple.r.z) #z vals
myOrbitTime.append(orbitPineapple_J2000time) #J2000 time vals
myOrbitJ2000Time.append(orbitPineapple.t)
#plot(orbitPineapple)
#print('X:', 'Length:', len(myOrbitX))
#print(myOrbitX)
#print('Y:', 'Length:',len(myOrbitY))
#print(myOrbitY)
#print('Z:', 'Length:', len(myOrbitZ))
#print(myOrbitZ)
#print('J2000 Time:', 'Length:',len(myOrbitTime))
#print(myOrbitTime)
'''Because the myOrbitTime is only the time between each step to be the sum of itself plus
all the previous times. And then I need to convert that time from seconds after J2000 to UTC.'''
myT = [] #UTC time list
for i in range(propNum):
myT.append((Time(2000, format='jyear') + TimeDelta(myOrbitTime[i]*u.s)).iso) #Convert time from J2000 to UTC
#print('UTC Time List Length:', len(myT))
#print('UTC Times:', myT)
'''Now I have xyz and time for each propagation step and need to convert the coordinates from
ECI to Lat, Lon, & Alt'''
now = [] #UTC time at each propagation step
xyz =[] #Xyz coordinates from OrbitalPy initial orbit propagation
cartrep = [] #Cartesian Representation
gcrs = [] #Geocentric Celestial Reference System/Geocentric Equatorial Inertial, the default coord system of OrbitalPy
itrs =[] #International Terrestrial Reference System coordinates
lat = [] #Longitude of the location, for the default ellipsoid
lon = [] #Longitude of the location, for the default ellipsoid
alt = [] #Height of the location, for the default ellipsoid
for i in range(propNum):
xyz = (myOrbitX[i], myOrbitY[i], myOrbitZ[i]) #Xyz coord for each prop. step
now = time.Time(myT[i]) #UTC time at each propagation step
cartrep = coord.CartesianRepresentation(*xyz, unit=u.m) #Add units of [m] to xyz
gcrs = coord.GCRS(cartrep, obstime=time.Time(myT[i])) #Let AstroPy know xyz is in GCRS
itrs = gcrs.transform_to(coord.ITRS(obstime=time.Time(myT[i]))) #Convert GCRS to ITRS
loc = coord.EarthLocation(*itrs.cartesian.xyz) #Get lat/lon/height from ITRS
lat.append(loc.lat.value) #Create latitude list
lon.append(loc.lon.value) #Create longitude list
alt.append(loc.height.value)
#print('Current Time:', now)
#print('')
#print('UTC Time:')
#print(myT)
print('myOrbitJ2000Time', myOrbitJ2000Time)
print('')
print('Lat:')
print('')
print(lat)
print('')
print('Lon:')
print(lon)
print('')
print('Alt:')
print(alt)
orbitPropandcoordTrans(5,-34963095,0.0073662,51.5946,154.8079,103.6176,257.3038,15.92610159)
最简单的方法可能是使用zip():
data = zip(myOrbitJ2000Time, lat, lon, alt)
所以“数据”将具有您想要序列化的格式。别忘了序列化您想要的头。在
为了回答您的一般性问题,您可以先创建Numpy数组,而不是定义一堆Python列表(这些列表的附加和处理速度都很慢,尤其是在您扩展分辨率时处理大量的值时)。初始化Numpy数组时,通常需要指定数组的大小,但在这种情况下,这很简单,因为您知道需要多少次传播。例如:
创建一个由
propNum
浮点值组成的空Numpy数组(这里的“empty”表示数组没有用任何值初始化这是创建一个新数组的最快方法,因为我们稍后将填充其中的所有值。在然后在循环中,而不是使用
orbitX.append(x)
,而是在数组中分配与当前记号相对应的值:orbitX[i] = x
。其他情况也一样。在还有几种方法可以输出数据,但是使用Astropy^{} 提供了很大的灵活性。可以创建包含所需列的表,例如:
^{pr2}$拥有Table对象的好处是有大量的输出格式选项。E、 g.在Python提示符下:
要先输出到文件,您必须考虑如何格式化数据。已经有很多常见的数据格式可以考虑使用,但这取决于数据的用途以及谁将使用它(“.dat file”在文件格式方面没有任何含义;或者更确切地说,它可能意味着任何)。但在您的问题中,您指出您想要的是所谓的“逗号分隔值”(CSV),其中数据的格式是按列向下排列,行中的每个值都用逗号分隔。
Table
类可以很容易地输出CSV(以及任何变体):不过,还有很多其他的选择。例如,如果希望数据在对齐列中格式良好,也可以这样做。你可以阅读更多关于它的here。(另外,我建议,如果您想要纯文本文件格式,可以尝试
ascii.ecsv
,它的优点是它输出额外的元数据,从而可以轻松地将其读回Astropy表中):另一件不相关的事我会注意到,除了单个值之外,Astropy中的许多对象可以对整个数组进行操作,这通常会更有效,尤其是对于许多值。特别是,您有一个Python循环:
这可以完全重写,无需显式循环,而是将它们视为坐标的数组。例如:
用这个我们也可以得到所有的坐标。E、 g.:
一般来说,这是处理大量坐标值的一种更有效的方法。类似地,在转换原始代码中的
myT
时,可以创建一个单独的TimeDelta
数组并将其添加到初始时间中,而不是在所有时间偏移上循环。在不幸的是,我不是
orbital
软件包的专家,但它似乎不像人们希望的那样容易获得轨道上不同点的坐标。就像应该有的。在相关问题 更多 >
编程相关推荐