I created a function to generate and propagate a satellite orbit. Now, I want to save all my data in a .dat file. I'm not sure how many for loops I need or quite how to do it.
I want time, latitude, longitude, and altitude all on one line for each propagation step.
Code for data:
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]
Desired Output in .dat file:
J2000 Time, Lat, Lon, Alt
1085.0, 48.6, -144.2, 264779.5
2170.0, 26.5, -50.4, 262446.1
3255.0, -28.7, -1.6, 319661.8
4340.1, -48.0, 91.5, 355717.2
5425.1, 0.1, 152.8, 06129.0
Full Code:
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)
To answer your general question, instead of defining a bunch of Python lists (which are slow to append to and work with, especially when dealing with a large number of values as you scale up your resolution), you can start by creating Numpy arrays instead of lists. When initializing Numpy arrays it's usually necessary to specify the size of the array, but in this case that's easy since you know how many propagations you want. For example:
>>> orbitX = np.empty(propNum, dtype=float)
creates an empty Numpy array of propNum floating-point values (here "empty" means the array is just not initialized with any values--this is the quickest way to create a new array as we're going to fill all the values in it later anyways.
Then in your loop instead of using orbitX.append(x) you would assign to the value in the array corresponding to the current tick: orbitX[i] = x. Same for the other cases.
Then there are several possibilities for how to output your data, but using the Astropy Table affords a lot of flexibility. You can create a table containing the columns you want easily like:
>>> from astropy.table import Table
>>> table = Table([J2000, lat, lon, alt], names=('J2000', 'lat', 'lon', 'alt'))
The nice thing about having the Table object is that there are a ton of options for output format. E.g. at the Python prompt:
>>> table
J2000 lat lon alt
float64 float64 float64 float64
------------- -------------- -------------- -------------
1085.01128806 48.5487129749 -144.175054697 264779.500823
2170.02257613 26.5068122883 -50.3805485685 262446.085716
3255.03386419 -28.7915478311 -1.6090285674 319661.817451
4340.04515225 -48.0536526356 91.5416828221 355717.274021
5425.05644032 0.084422392655 152.811717713 306129.02576
To output to a file first you have to consider how you want the data formatted. There are already many common data formats that you can consider using, but it depends on what the data is for and who will be consuming it (".dat file" doesn't mean anything in terms of file formats; or rather, it could mean anything). But in your question you indicated that what you want is what is called "comma-separated values" (CSV) where the data is formatted with columns going down, with each value within a row separated by a comma. The Table class can output CSV (and any variant) quite easily:
>>> import sys
>>> table.write(sys.stdout, format='ascii.csv') # Note: I'm just using sys.stdout for demonstration purposes; normally you would give a filename
J2000,lat,lon,alt
1085.011288063746,48.54871297493748,-144.17505469691633,264779.5008225624
2170.022576127492,26.506812288280788,-50.38054856853237,262446.0857159357
3255.0338641912376,-28.79154783108773,-1.6090285674024463,319661.81745081506
4340.045152254984,-48.05365263557444,91.54168282208444,355717.2740210131
5425.05644031873,0.08442239265500713,152.81171771323176,306129.0257600865
There are lots of other options though. For example, if you want the data nicely formatted in justified columns you could do that too. You can read more about it here. (Also, I would suggest that if you want a plain-text file format you try ascii.ecsv which has the advantage that it outputs additional metadata that allows it to be easily read back into an Astropy Table):
>>> table.write(sys.stdout, format='ascii.ecsv')
# %ECSV 0.9
# ---
# datatype:
# - {name: J2000, datatype: float64}
# - {name: lat, datatype: float64}
# - {name: lon, datatype: float64}
# - {name: alt, datatype: float64}
# schema: astropy-2.0
J2000 lat lon alt
1085.01128806 48.5487129749 -144.175054697 264779.500823
2170.02257613 26.5068122883 -50.3805485685 262446.085716
3255.03386419 -28.7915478311 -1.6090285674 319661.817451
4340.04515225 -48.0536526356 91.5416828221 355717.274021
5425.05644032 0.084422392655 152.811717713 306129.02576
Another unrelated thing I'll note is that many objects in Astropy can operate over whole arrays, in addition to single values, and this can often be much more efficient, especially for many values. In particular, you have this Python loop:
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)
This can be rewritten entirely without an explicit loop by treating them as arrays of coordinates instead. For example:
>>> times = time.Time(myT) # All the times, not just a single one
>>> cartrep = coord.CartesianRepresentation(orbitX, orbitY, orbitZ, unit=u.m) # Again, an array of coordinates
>>> gcrs = coord.GCRS(cartrep, obstime=times)
>>> itrs = gcrs.transform_to(coord.ITRS(obstime=ts))
>>> loc = coord.EarthLocation(*itrs.cartesian.xyz) # I'm not sure this is the most efficient way to do this but I'm not an expert on the coordinates package
With this we can get all the coordinates as arrays as well. E.g.:
>>> loc.lat
<Latitude [ 48.54871297, 26.50681229,-28.79154783,-48.05365264,
0.08442239] deg>
So this is in general a more efficient way to work with large numbers of coordinate values. Similarly when converting the myT in your original code, instead of looping over all time offsets you can create a single TimeDelta array and add that to your initial time.
Unfortunately I'm not an expert on the orbital package, but it doesn't seem like it makes it as easy as one would like to get the coordinates at different points on an orbit. ISTM like there should be.