Hi. I have tried to plot the coastlines as a globe as you did in this video but it is not working as expected. Actually I converted the coastlines csv lats and longs into rectangular coords (x,y and z) and plotted it in 3d. It does not show as expected. I am sure I am wrongly doing the conversion. Could you please advise and share how you did it? Thanks.
Yes so this process is different than doing it in 2D. It doesn't look right when you try to plot 60k points in 3D, and its really slow. So instead I reduced it down to ~1-2k points, and separated each large land mass (africa, americas, australia, etc.) So I wrote a quick program to open the image in matplotlib and manually click on the coastlines to get the data that way. I just posted this new data in the same repository I showed before in the coastlines_3d directory. Here is the link: github.com/alfonsogonzalez/orbital_mechanics_with_python And then to plot them, I plot them as lines, not points. Here is a snippet of the software to load in the lat long coordinates and convert to ECEF vectors. This is within the OrbitalAnimator class I have so there are some extra names but I think you'll be able to figure out the concept. The first two lines are just global variables and the indented part is a method of the class: LAND_MASSES = [ 'africa', 'americas', 'antarctica', 'australia', 'eurasia', 'greenland', 'new_zealand' ] COASTLINE_DATA_DIR = '/home/alfonso/AWP/earth_coastline_data' def load_body_fixed_land_masses( self ): for name in LAND_MASSES: print( 'Loading %s..' % name ) latlons = np.genfromtxt( os.path.join( COASTLINE_DATA_DIR, name + '.csv' ), delimiter = ',' ) * nt.d2r self.land_masses[ name ] = np.zeros( ( latlons.shape[ 0 ], 3 ) ) for n in range( latlons.shape[ 0 ] ): self.land_masses[ name ][ n ] = spice.latrec( self.config[ 'cb_radius' ], latlons[ n, 0 ], latlons[ n, 1 ] ) And then in the plotting function: if self.config[ 'show_coastlines' ]: for name in LAND_MASSES: if self.config[ 'frame' ] == 'inertial': rs = self.land_masses[ name ][ step ] * 1.01 else: rs = self.land_masses[ name ] * 1.01 plt.plot( rs[ :, 0 ], rs[ :, 1 ], rs[ :, 2 ], self.config[ 'land_mass_color' ], markersize = 0.5, zorder = 10 ) So I implemented the choice of whether im doing the animation in the ECI or ECEF frames, because if its in ECEF I only have to load one frame of the coastlines since they are stationary. But for ECI, the earth is rotating so I have to load in the data at each time step. Also, I only have this in the animations, not for stationary 3D plots. Theres a lot in there so let me know if you have more questions
@@alfonsogonzalez-astrodynam2207 Thanks man, this works great. I did it for the ecef frame only. I have 2 questions: 1) why do you multiply the coordinates by a factor 1.01. 2) I skipped the inertial plot, i.e. the line "rs = self.land_masses[ name ][ step ] * 1.01". I do understand that the data for each step need to be loaded as earth is rotating. However I do not understand that single line with [step]. If step is a number, then rs will be a single x,y,z point.
For point 2, alternatively I have used the spice axisar function to obtain a rotation matrix about the z-axis from the angle calculated for the tstep. Using a for loop i kept rotating it to the last step. Then I took the final xs and ys. The reason for the for loop is to be able to use each rotation step for an animation when required later. Not sure if this is the right approach but it seems to work.
1. I added the 1.01 while I was experimenting with how the plots actually looked, and the zorder. I don't think its necessary to have that. Instead, I ended up making Earth more transparent. Here is the part where I do that, where the "alpha" parameter is a number between 0 and 1, and the closer it is to 0, the more transparent Earth is (by default I have it set at 0.2): _u, _v = np.mgrid[ 0:2*np.pi:20j, 0:np.pi:10j ] _x = self.config[ 'cb_radius' ] * np.cos( _u ) * np.sin( _v ) _y = self.config[ 'cb_radius' ] * np.sin( _u ) * np.sin( _v ) _z = self.config[ 'cb_radius' ] * np.cos( _v ) ax.plot_surface( _x, _y, _z, cmap = self.config[ 'cb_cmap' ], zorder = self.config[ 'cb_zorder' ], alpha = self.config[ 'cb_alpha' ] ) 2. That extra [step] is because I load in the inertial land masses differently than the ECEF. Here is the method inside OrbitalAnimator that loads in inertial land masses: def load_inertial_land_masses( self ): for name in LAND_MASSES: print( 'Loading %s..' % name ) latlons = np.genfromtxt( os.path.join( COASTLINE_DATA_DIR, name + '.csv' ), delimiter = ',' ) * nt.d2r self.land_masses[ name ] = np.zeros( ( self.steps, latlons.shape[ 0 ], 3 ) ) for step in range( self.steps ): for n in range( latlons.shape[ 0 ] ): vector = spice.latrec( self.config[ 'cb_radius' ], latlons[ n, 0 ], latlons[ n, 1 ] ) matrix = spice.pxform( self.config[ 'bf_frame' ], self.config[ 'inert_frame' ], self.tspan [ step ] ) vector = np.dot( matrix, vector ) self.land_masses[ name ][ step, n ] = vector
For this case, this will work, since the earth is spinning around its ECEF z-axis. So actually, that might be faster than finding the ECEF frame using spice.pxform, so that was a good observation! I've covered the principal axis rotation matrices in the numerical methods with python series, so you may be interested in checking those out to learn more about those rotations. But here is how I have a z-axis rotation matrix in the numerical_tools.py file: def Cz( a ): ''' Principal Z axis rotation matrix by an angle ''' return np.array( [ [ math.cos( a ), -math.sin( a ), 0 ], [ math.sin( a ), math.cos( a ), 0 ], [ 0, 0, 1 ] ] )
Is a numerical integrator same as a propogator? I feel like all propogators should arrive at same solution if same forces are included and the steps are small enough and numerical precision is not an issue.
It just depends how precise the application needs to be. But my understanding is that even with relative and absolute tolerances the same, different solvers will give different solutions (even if they are really small, say 1E-15). One could do a test with 2 body dynamics of the different SciPy ODE solvers to confirm or disprove that: docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
Hi. I have tried to plot the coastlines as a globe as you did in this video but it is not working as expected. Actually I converted the coastlines csv lats and longs into rectangular coords (x,y and z) and plotted it in 3d. It does not show as expected. I am sure I am wrongly doing the conversion. Could you please advise and share how you did it? Thanks.
Yes so this process is different than doing it in 2D. It doesn't look right when you try to plot 60k points in 3D, and its really slow. So instead I reduced it down to ~1-2k points, and separated each large land mass (africa, americas, australia, etc.)
So I wrote a quick program to open the image in matplotlib and manually click on the coastlines to get the data that way. I just posted this new data in the same repository I showed before in the coastlines_3d directory. Here is the link: github.com/alfonsogonzalez/orbital_mechanics_with_python
And then to plot them, I plot them as lines, not points.
Here is a snippet of the software to load in the lat long coordinates and convert to ECEF vectors. This is within the OrbitalAnimator class I have so there are some extra names but I think you'll be able to figure out the concept. The first two lines are just global variables and the indented part is a method of the class:
LAND_MASSES = [ 'africa', 'americas', 'antarctica', 'australia', 'eurasia',
'greenland', 'new_zealand' ]
COASTLINE_DATA_DIR = '/home/alfonso/AWP/earth_coastline_data'
def load_body_fixed_land_masses( self ):
for name in LAND_MASSES:
print( 'Loading %s..' % name )
latlons = np.genfromtxt( os.path.join( COASTLINE_DATA_DIR, name + '.csv' ),
delimiter = ',' ) * nt.d2r
self.land_masses[ name ] = np.zeros( ( latlons.shape[ 0 ], 3 ) )
for n in range( latlons.shape[ 0 ] ):
self.land_masses[ name ][ n ] = spice.latrec(
self.config[ 'cb_radius' ],
latlons[ n, 0 ], latlons[ n, 1 ] )
And then in the plotting function:
if self.config[ 'show_coastlines' ]:
for name in LAND_MASSES:
if self.config[ 'frame' ] == 'inertial':
rs = self.land_masses[ name ][ step ] * 1.01
else:
rs = self.land_masses[ name ] * 1.01
plt.plot( rs[ :, 0 ], rs[ :, 1 ], rs[ :, 2 ],
self.config[ 'land_mass_color' ], markersize = 0.5, zorder = 10 )
So I implemented the choice of whether im doing the animation in the ECI or ECEF frames, because if its in ECEF I only have to load one frame of the coastlines since they are stationary. But for ECI, the earth is rotating so I have to load in the data at each time step.
Also, I only have this in the animations, not for stationary 3D plots.
Theres a lot in there so let me know if you have more questions
@@alfonsogonzalez-astrodynam2207 Thanks man, this works great. I did it for the ecef frame only. I have 2 questions: 1) why do you multiply the coordinates by a factor 1.01. 2) I skipped the inertial plot, i.e. the line "rs = self.land_masses[ name ][ step ] * 1.01". I do understand that the data for each step need to be loaded as earth is rotating. However I do not understand that single line with [step]. If step is a number, then rs will be a single x,y,z point.
For point 2, alternatively I have used the spice axisar function to obtain a rotation matrix about the z-axis from the angle calculated for the tstep. Using a for loop i kept rotating it to the last step. Then I took the final xs and ys.
The reason for the for loop is to be able to use each rotation step for an animation when required later.
Not sure if this is the right approach but it seems to work.
1. I added the 1.01 while I was experimenting with how the plots actually looked, and the zorder. I don't think its necessary to have that. Instead, I ended up making Earth more transparent. Here is the part where I do that, where the "alpha" parameter is a number between 0 and 1, and the closer it is to 0, the more transparent Earth is (by default I have it set at 0.2):
_u, _v = np.mgrid[ 0:2*np.pi:20j, 0:np.pi:10j ]
_x = self.config[ 'cb_radius' ] * np.cos( _u ) * np.sin( _v )
_y = self.config[ 'cb_radius' ] * np.sin( _u ) * np.sin( _v )
_z = self.config[ 'cb_radius' ] * np.cos( _v )
ax.plot_surface( _x, _y, _z, cmap = self.config[ 'cb_cmap' ],
zorder = self.config[ 'cb_zorder' ], alpha = self.config[ 'cb_alpha' ] )
2. That extra [step] is because I load in the inertial land masses differently than the ECEF. Here is the method inside OrbitalAnimator that loads in inertial land masses:
def load_inertial_land_masses( self ):
for name in LAND_MASSES:
print( 'Loading %s..' % name )
latlons = np.genfromtxt( os.path.join( COASTLINE_DATA_DIR, name + '.csv' ),
delimiter = ',' ) * nt.d2r
self.land_masses[ name ] = np.zeros(
( self.steps, latlons.shape[ 0 ], 3 ) )
for step in range( self.steps ):
for n in range( latlons.shape[ 0 ] ):
vector = spice.latrec( self.config[ 'cb_radius' ],
latlons[ n, 0 ], latlons[ n, 1 ] )
matrix = spice.pxform(
self.config[ 'bf_frame' ],
self.config[ 'inert_frame' ],
self.tspan [ step ] )
vector = np.dot( matrix, vector )
self.land_masses[ name ][ step, n ] = vector
For this case, this will work, since the earth is spinning around its ECEF z-axis. So actually, that might be faster than finding the ECEF frame using spice.pxform, so that was a good observation! I've covered the principal axis rotation matrices in the numerical methods with python series, so you may be interested in checking those out to learn more about those rotations. But here is how I have a z-axis rotation matrix in the numerical_tools.py file:
def Cz( a ):
'''
Principal Z axis rotation matrix by an angle
'''
return np.array( [
[ math.cos( a ), -math.sin( a ), 0 ],
[ math.sin( a ), math.cos( a ), 0 ],
[ 0, 0, 1 ]
] )
Is a numerical integrator same as a propogator? I feel like all propogators should arrive at same solution if same forces are included and the steps are small enough and numerical precision is not an issue.
It just depends how precise the application needs to be. But my understanding is that even with relative and absolute tolerances the same, different solvers will give different solutions (even if they are really small, say 1E-15). One could do a test with 2 body dynamics of the different SciPy ODE solvers to confirm or disprove that: docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
nice vid! btw, can you make a video about instantaneous launch windows and reentry?
I can move it up in my to-do list, but there are a lot of videos that need to be made before that / a bit complicated.
@@alfonsogonzalez-astrodynam2207 ill look forward to it then, godspeed!
The sound is low ...