Open gregoriomarchesini opened 2 years ago
Hi Greg, we had a long weekend holiday here, hence the late reply. This looks like a nice piece of functionality again :) One question on the 'settings' input. Wouldn't it be easier to have each of the settings entries be an explicit input (with default) for the function itself? That way:
But, as you know, I am not that good a Python programmer, nor am I very familiar with common practices. Let me know what you think about the above suggestion
Absolutely
I will correct it soon! It is indeed a good suggestion
Super I will repost here the updated version
Hi again
I have corrected the function and I believe now it could be more useful. I inserted the parameters as you asked in the function signature
I have also created a minimal working example so that you can test if everything is working as you wish
import numpy as np
from matplotlib import pyplot as plt
from tudatpy.io import save2txt
from tudatpy.util import result2array
from tudatpy.kernel import constants
from tudatpy.kernel.interface import spice
from tudatpy.kernel import numerical_simulation # KC: newly added by me
from tudatpy.kernel.numerical_simulation import environment_setup
from tudatpy.kernel.numerical_simulation import propagation_setup
from scipy.spatial.transform import Rotation
import cartopy
def split_history(state_history:dict,n_bodies:int):
"""
Parameters
-----------
state_history (dict ) : dictionary containing the state_history
of N propagated bodies.
The dictionary is defined as a standard
state_history following the TUdat convention.
the keys of the dictionay are the time stamps
corresponding to each state and the argument is
the state of the body at the given time stamp
n_rotation (int) : number of propagated bodies
Returns
-----------
state_hostory_book (list) : list of state_history dictionaries for each propagated body.
The order of the list is given directly by the order
in which the bodies are propagated.
Descriptiom
-----------
Creates a list of state histories based on the unified state_history
from from the propagation of N satellites
"""
state_history_book = []
time = [key for key in state_history.keys() ]
n_variables = len(state_history[time[0]])
step = int(n_variables/n_bodies)
if n_variables%n_bodies :
print(" ERROR : Not possible to split. Number of variables not divisible by number of bodies")
raise
for n in range(n_bodies) :
current_history = dict()
for t in state_history.keys():
current_history[t]= state_history[t][step*n:step*n+step]
state_history_book.append(current_history)
return state_history_book
def groundtrack_generator(groundtracks_list :list = [],
labels :list = [],
projection :str = "PlateCarree",
colors :list = [],
linewidth :float = 1 ,
grid :bool = True,
grid_linewidth :float = 1,
grid_color :str = "gray",
grid_alpha :float = 0.3) :
"""
Parameters
----------
groundtracks_list (List) : List of numpy.ndarray[[N,2]] defining the latitude and longitude of each groundtrack
col[0] = longitude
col[1] = latitude.
Units are radiants.
labels (List ) : label per groundtrack
projection (Str ) : projection to be used (see available projections below)
colors (List ) : list of colours per groundtrack
linewidth (Float) : linewidth of the groundtracks
grid (Bool ) : activates grid on map, defining latitude and longitude reference
grid_linewidth (Float) : grid width
grid_color (Str ) : grid color (refer to https://matplotlib.org/3.5.0/tutorials/colors/colors.html#:~:text=Single%20character%20shorthand%20notation%20for%20some%20basic%20colors. for a list of colors)
grid_alpha (Float) : grid alpha
Returns
--------
Returns plot of a given set of ground tracks as definied by the groundtrack_list
Description :
-------------
plots a set of ground tracks over a specified map.
have a look at :
https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html
for all maps projections and types
Maps
1 PlateCarree
2 Mollweide
3 AlbersEqualArea
4 AzimuthalEquidistant
5 LambertConformal
6 LambertCylindrical
7 Mercator
8 Miller
9 Mollweide
10 Orthographic
11 Robinson
12 Sinusoidal
13 Stereographic
14 TransverseMercator
15 InterruptedGoodeHomolosine
16 RotatedPole
"""
try :
map_projection = getattr(cartopy.crs,projection)
except AttributeError :
raise ValueError("unknown map type {}.\n Check again the correct name of the map you want to apply".format(projection))
# start defining the figure
figure = plt.figure(figsize=(8,6))
ax = plt.axes(projection=map_projection())
rad2deg = 180/np.pi
if len(colors) == 0 :
colors = [(.0,.0,.0,1.)]*len(groundtracks_list)
elif len(colors) != len(groundtracks_list):
raise ValueError("colors and groundtracks_list have different lenghts")
if len(labels) == 0 :
labels = ["sat{}".format(n) for n in range(len(groundtracks_list))]
elif len(labels) != len(groundtracks_list):
raise ValueError("labels sand groundtracks_list have different lenghts")
for groundtrack,color,label in zip(groundtracks_list,colors,labels) :
# plot start position with bigger width
long = groundtrack[:,0]*rad2deg
lat = groundtrack[:,1]*rad2deg
plt.plot(long,lat,
color = color,
linewidth = linewidth,
label = label,
transform=cartopy.crs.Geodetic())
plt.scatter(long[0],lat[0],
linewidth=linewidth*5,
color = color,
transform=cartopy.crs.Geodetic(),
marker = "*")
ax.set_global()
if grid :
ax.gridlines(color =grid_color,
linestyle ="--",
draw_labels=True,
linewidth =grid_linewidth,
alpha =grid_alpha)
try :
ax.stock_img()
except :
pass
ax.legend()
# define useful constants
julian_day = constants.JULIAN_DAY
hours = 60*60
minutes = 60
rad2deg = 180/np.pi
# start simulation day 20/10/1998
# simulate for 10 days
simulation_start_epoch = 2451106.50000
simulation_end_epoch = simulation_start_epoch + 4*hours
# Load SPICE.
spice.load_standard_kernels()
# Create settings for celestial bodies
bodies_to_create = ['Earth'] # this must have a list of all the planets to create
global_frame_origin = 'Earth' # this is the origin of the refernce system
global_frame_orientation = 'J2000' # orinetation of the reference system
body_settings = environment_setup.get_default_body_settings(
bodies_to_create, global_frame_origin, global_frame_orientation) # body settings taken from SPICE.
# Create environment
bodies = environment_setup.create_system_of_bodies(body_settings)
# define constellation of satellites
number_of_satellites = 4
satellites_names = ["sat{}".format(str(num)) for num in range(number_of_satellites)]
for name in satellites_names :
bodies.create_empty_body(name)
# Set mass of vehicle
bodies.get_body(name).mass = 500.0
###########################################################################
# CREATE ACCELERATIONS ####################################################
###########################################################################
# Define bodies that are propagated, and their central bodies of propagation.
bodies_to_propagate = satellites_names
central_bodies = ['Earth']*len(bodies_to_propagate) # bodies which are contributing with the full gravity and are not only perturbation
# Define accelerations acting on vehicle.
acceleration_settings_satellite = dict(
Earth=[propagation_setup.acceleration.spherical_harmonic_gravity(2,2)]
)
global_acceleration_settings = {name: acceleration_settings_satellite for name in satellites_names}
# Create acceleration models.
acceleration_models = propagation_setup.create_acceleration_models(
bodies, global_acceleration_settings, bodies_to_propagate, central_bodies)
# Define initial state.
standard_initial_state = np.array([5.194681552632270E+06,-1.833738414039049E+06,4.172975147134797E+06,
4.556144583672027E+03,4.992919206236066E+03,-3.471626733390986E+03])[:,np.newaxis]
phi_intervals = np.linspace(0,1.9*np.pi,number_of_satellites)
multidimenasional_state = []
dependent_variables = []
for phi,name in zip(phi_intervals,satellites_names) :
rot_mat = Rotation.from_rotvec(np.array([0,0,1])*phi).as_matrix()
full_state_rot = np.eye(6)
full_state_rot[:3,:3] = rot_mat
full_state_rot[3:,3:] = rot_mat
updated_state = (full_state_rot @standard_initial_state)
multidimenasional_state.append(updated_state)
dependent_variables.append(propagation_setup.dependent_variable.longitude(name,'Earth'))
dependent_variables.append(propagation_setup.dependent_variable.geodetic_latitude(name,'Earth'))
multidimenasional_state = np.concatenate(tuple(multidimenasional_state),axis=0)
# Create propagation settings.
termination_settings = propagation_setup.propagator.time_termination( simulation_end_epoch )
propagator_settings = propagation_setup.propagator.translational(
central_bodies,
acceleration_models,
bodies_to_propagate,
multidimenasional_state,
termination_settings,
output_variables = dependent_variables
)
# Create numerical integrator settings.
fixed_step_size = 10.0
integrator_settings = propagation_setup.integrator.runge_kutta_4(
simulation_start_epoch,
fixed_step_size
)
# Create simulation object and propagate dynamics.
dynamics_simulator = numerical_simulation.SingleArcSimulator(
bodies, integrator_settings, propagator_settings )
# create list of state dictionaries for each satellite
state_history_book = split_history(dynamics_simulator.state_history,number_of_satellites) #list of state histories
state_history_array_book = [result2array(state_history) for state_history in state_history_book]
dependat_variables_array = result2array(dynamics_simulator.dependent_variable_history)
# divide long/lat info from rotation matrix
longlat_coordinates = dependat_variables_array[::20,1:number_of_satellites*2+1]
groundtracks_list = [longlat_coordinates[:,ii:ii+2] for ii in range(0,number_of_satellites*2,2)]
labels = ["sat{}".format(n) for n in range(len(groundtracks_list))]
color_list = [(.0,.0,.0,1.) for ii in range(number_of_satellites)]
map_settings = {"projection" : "Sinusoidal",
"color" : color_list,
"grid" : True,
"grid_alpha" : 0}
test_maps = ["Sinusoidal",
"Stereographic",
"LambertCylindrical",
"Mercator",
"Miller",
"Mollweide",
"Orthographic",
"InterruptedGoodeHomolosine"]
for map in test_maps :
groundtrack_generator(groundtracks_list = groundtracks_list,
colors = color_list,
labels = labels,
projection = map,
linewidth = 1 ,
grid = True,
grid_linewidth = 1,
grid_color = "gray",
grid_alpha = 0.3)
plt.show()
In addition to Tudat, the packages scipy
and cartopy
are required.
But nothing more special than this
I hope the function will be useful in the future, and if not you can just ignore this current issue.
Thank you for your time and attention Have a nice evening
Greg
Hi Greg, somehow we lost track of this code your wrote, our bad! I think we can split this in two:
Dear Dominic
I will do my best to find this piece of code and restore it. It has probably been an year since the last time I opened it. This week I will try to conclude it. Unfortunately my week is super busy until Friday. In the weekend I can try to fix it and see if the code is still working.
Thank you for your time and attention Best Greg
On 14 Dec 2022, at 14:20, DominicDirkx @.***> wrote:
Hi Greg, somehow we lost track of this code your wrote, our bad! I think we can split this in two:
Add your new Python functions into Tudatpy Add the functionality into an example application. Perhaps this one: tudat-team/tudatpy-examples#14 https://github.com/tudat-team/tudatpy-examples/issues/14 would be a good choice — Reply to this email directly, view it on GitHub https://github.com/tudat-team/tudatpy/issues/173, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQ76ADECCSZWSF5CH55IDYTWNHCTTANCNFSM5TP2ZN7A. You are receiving this because you authored the thread.
Hi Greg, no need to hurry! I didn't mean to ask you to do this, I think we can merge your new Python functions into Tudat already (if you want 'credit' in the tudatpy commit history, you can open a pull request, though). The use in the Galileo example will come when we properly write that example.
Hey Dominic
Ok perfect. There is no need for the pull request. I am very happy that my code gets used and I hope it will serve its purpose. In case you need help, I will do my best to help you in time. Super thank you and have a great day
Best Greg
On 15 Dec 2022, at 15:37, DominicDirkx @.***> wrote:
Hi Greg, no need to hurry! I didn't mean to ask you to do this, I think we can merge your new Python functions into Tudat already (if you want 'credit' in the tudatpy commit history, you can open a pull request, though). The use in the Galileo example will come when we properly write that example.
— Reply to this email directly, view it on GitHub https://github.com/tudat-team/tudatpy/issues/173, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQ76ADCK7SCXEUWGGXDU5YLWNMULBANCNFSM5TP2ZN7A. You are receiving this because you authored the thread.
Hi I have developed a small piece of code that can visualize ground tracks using cartopy The available list of map is listed in the documentation of the function.
I tried it and it works for me. You have to install cartopy from conda, but nothing else is required (apart from numpy, but this is pretty normal to have).
If you think it is useful, let me know. I hope it will help
Have a super day!