Swarm-NG
1.1
|
In this tutorial, we explore the different ways to examine an ensemble before and after integration.
We also demonstrate how to convert the position and velocity of planets into orbital elements and finally, we plot the difference in semi-major axis of planets.
Source code for this tutorial can be found at py/integration_plot_tutorial.py
The set-up is quite similar to approach taken in Resume an integration.
Import the packages and parse the command-line arguments
Next, load the config file and the initial conditions
Integration is the same as the basic tutorial
At this point, the integration is completed. For analysis we have two ensembles: ref
, the reference ensemble that keeps a copy of initial conditions and ens
, the working copy that was simulated to the destination_time
.
In this example, we want to find the semi-major axis of planets and find the evolution of this variable over time.
The basic access methods only give us position and velocity of planets. To calculate the orbital elements, we use the function swarmng.keplerian_for_cartesian.
The function extarct_semi_major
takes an ensemble and calculate the semi-major axis for all planets in all of systems and returns a two-dimentional array of the results. For easier plotting, the first dimension is the planet number and the second dimension is the system number.
For large two-dimensional arrays, we use NumPy. It is much faster than regular Python arrays.
The core function to calculate the semi-major axis is swarmng.keplerian_for_cartesian. It takes two 3-tuples for the information of the planet and the star and returns a tuple of all orbital elemets. For more details refer to documentation for swarmng.keplerian_for_cartesian.
We calculate semi-major axis values for both ensembles: initial and final conditions (final means after the simulation).
The results are basically two-dimensional arrays of same shape. The next line normalizes the changes on the values. Note that the operations are carried out element-wise, so that the line of code here is equivalent to executing change[i,j] = (final[i,j]-initial[i,j])/initial[i,j]
for all valid i
,j
.
For every planet, we plot a scatter plot of change of semi-major axis where each data point represents the change in one of the systems. The plots from different planets are overlapped but distinguished by different markers.
X-axis represents the initial value of the semi-major axis calculated from the initial conditions.
Y-axis represents the normalized change in the semi-major axis for individual planets in different systems.
The generated plot is saved as semi_major_axis_comparison
in the currenty working directory.
In this, tutorial we explored how to extract information from an ensemble in a real application. You can change the tutorial to plot other orbital elements, or change the simulation part to more complicated one.
We used NumPy and MatPlotLib packages, it is not possible to cover those packages in this tutorial. It is best to look for comprehensive tutorials of these packages on their respective websites.
Where to go from here: