Merge pull request #48 from larshinueber/feature/minor-fixes

Minor fixes

.gitignore

@@ -33,6 +33,8 @@ cmake-build-*/
 .pydevproject
 *.sublime-project
 *.sublime-workspace
+.vscode/
+
 # Rever
 rever/
 

estimation/covariance_estimated_parameters.ipynb

@@ -189,12 +189,15 @@
    "source": [
     "### Create the acceleration model\n",
     "Subsequently, all accelerations (and there settings) that act on `Delfi-C3` have to be defined. In particular, we will consider:\n",
+    "\n",
     "* Gravitational acceleration using a spherical harmonic approximation up to 8th degree and order for Earth.\n",
     "* Aerodynamic acceleration for Earth.\n",
     "* Gravitational acceleration using a simple point mass model for:\n",
+    "\n",
     "    - The Sun\n",
     "    - The Moon\n",
     "    - Mars\n",
+    "\n",
     "* Radiation pressure experienced by the spacecraft - shape-wise approximated as a spherical cannonball - due to the Sun.\n",
     "\n",
     "The defined acceleration settings are then applied to `Delfi-C3` by means of a dictionary, which is finally used as input to the propagation setup to create the acceleration models."
@@ -697,13 +700,10 @@
   }
  ],
  "metadata": {
-  "interpreter": {
-   "hash": "4a4d53b53330cd83e1499268313a4bcd5eafe4bf50523883929af79f2dd687b2"
-  },
   "kernelspec": {
-   "display_name": "tudatpy-examples",
+   "display_name": "tudat-space",
    "language": "python",
-   "name": "tudatpy-examples"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -715,7 +715,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.10.14"
   }
  },
  "nbformat": 4,

estimation/covariance_estimated_parameters.py

@@ -40,6 +40,7 @@
 from tudatpy.astro.time_conversion import DateTime
 from tudatpy.astro import element_conversion
 
+
 ## Configuration
 """
 First, NAIF's `SPICE` kernels are loaded, to make the positions of various bodies such as the Earth, the Sun, or the Moon known to `tudatpy`.
@@ -56,6 +57,7 @@
 simulation_start_epoch = DateTime(2000, 1, 1).epoch()
 simulation_end_epoch   = DateTime(2000, 1, 4).epoch()
 
+
 ## Set up the environment
 """
 We will now create and define the settings for the environment of our simulation. In particular, this covers the creation of (celestial) bodies, vehicle(s), and environment interfaces.
@@ -83,6 +85,7 @@
 # Create system of bodies
 bodies = environment_setup.create_system_of_bodies(body_settings)
 
+
 ### Create the vehicle and its environment interface
 """
 We will now create the satellite - called Delfi-C3 - for which an orbit will be simulated. Using an `empty_body` as a blank canvas for the satellite, we define mass of 400kg, a reference area (used both for aerodynamic and radiation pressure) of 4m$^2$, and a aerodynamic drag coefficient of 1.2. Idem for the radiation pressure coefficient. Finally, when setting up the radiation pressure interface, the Earth is set as a body that can occult the radiation emitted by the Sun.
@@ -111,6 +114,7 @@
 # Add the radiation pressure interface to the environment
 environment_setup.add_radiation_pressure_interface(bodies, "Delfi-C3", radiation_pressure_settings)
 
+
 ## Set up the propagation
 """
 Having the environment created, we will define the settings for the propagation of the spacecraft. First, we have to define the body to be propagated - here, the spacecraft - and the central body - here, Earth - with respect to which the state of the propagated body is defined.
@@ -122,15 +126,19 @@
 # Define central bodies of propagation
 central_bodies = ["Earth"]
 
+
 ### Create the acceleration model
 """
 Subsequently, all accelerations (and there settings) that act on `Delfi-C3` have to be defined. In particular, we will consider:
+
 * Gravitational acceleration using a spherical harmonic approximation up to 8th degree and order for Earth.
 * Aerodynamic acceleration for Earth.
 * Gravitational acceleration using a simple point mass model for:
+
     - The Sun
     - The Moon
     - Mars
+
 * Radiation pressure experienced by the spacecraft - shape-wise approximated as a spherical cannonball - due to the Sun.
 
 The defined acceleration settings are then applied to `Delfi-C3` by means of a dictionary, which is finally used as input to the propagation setup to create the acceleration models.
@@ -163,6 +171,7 @@
     bodies_to_propagate,
     central_bodies)
 
+
 ### Define the initial state
 """
 Realise that the initial state of the spacecraft always has to be provided as a cartesian state - i.e. in the form of a list with the first three elements representing the initial position, and the three remaining elements representing the initial velocity.
@@ -178,6 +187,7 @@
 delfi_ephemeris = environment.TleEphemeris( "Earth", "J2000", delfi_tle, False )
 initial_state = delfi_ephemeris.cartesian_state( simulation_start_epoch )
 
+
 ### Create the integrator settings
 """
 For the problem at hand, we will use an RKF78 integrator with a fixed step-size of 60 seconds. This can be achieved by tweaking the implemented RKF78 integrator with variable step-size such that both the minimum and maximum step-size is equal to 60 seconds and a tolerance of 1.0
@@ -188,6 +198,7 @@
     runge_kutta_fixed_step_size(initial_time_step=60.0,
                                 coefficient_set=propagation_setup.integrator.CoefficientSets.rkdp_87)
 
+
 ### Create the propagator settings
 """
 By combining all of the above-defined settings we can define the settings for the propagator to simulate the orbit of `Delfi-C3` around Earth. A termination condition needs to be defined so that the propagation stops as soon as the specified end epoch is reached. Finally, the translational propagator's settings are created.
@@ -207,6 +218,7 @@
     termination_condition
 )
 
+
 ## Set up the observations
 """
 Having set the underlying dynamical model of the simulated orbit, we can define the observational model. Generally, this entails the addition all required ground stations, the definition of the observation links and types, as well as the precise simulation settings.
@@ -232,6 +244,7 @@
     [station_altitude, delft_latitude, delft_longitude],
     element_conversion.geodetic_position_type)
 
+
 ### Define Observation Links and Types
 """
 To establish the links between our ground station and `Delfi-C3`, we will make use of the [observation module](https://py.api.tudat.space/en/latest/observation.html#observation) of tudat. During th link definition, each member is assigned a certain function within the link, for instance as "transmitter", "receiver", or "reflector". Once two (or more) members are connected to a link, they can be used to simulate observations along this particular link. The precise type of observation made along this link - e.g., range, range-rate, angular position, etc. - is then determined by the chosen observable type.
@@ -250,6 +263,7 @@
 link_definition = observation.LinkDefinition(link_ends)
 observation_settings_list = [observation.one_way_doppler_instantaneous(link_definition)]
 
+
 ### Define Observation Simulation Settings
 """
 We now have to define the times at which observations are to be simulated. To this end, we will define the settings for the simulation of the individual observations from the previously defined observation models. Bear in mind that these observation simulation settings are not to be confused with the ones to be used when setting up the estimator object, as done just above.
@@ -282,6 +296,7 @@
     [viability_setting]
 )
 
+
 ## Set up the estimation
 """
 Using the defined models for the environment, the propagator, and the observations, we can finally set the actual presentation up. In particular, this consists of defining all parameter that should be estimated, the creation of the estimator, and the simulation of the observations.
@@ -303,6 +318,7 @@
 # Create the parameters that will be estimated
 parameters_to_estimate = estimation_setup.create_parameter_set(parameter_settings, bodies)
 
+
 ### Creating the Estimator object
 """
 Ultimately, the `Estimator` object consolidates all relevant information required for the estimation of any system parameter:
@@ -321,6 +337,7 @@
     observation_settings_list,
     propagator_settings)
 
+
 ### Perform the observations simulation
 """
 Using the created `Estimator` object, we can perform the simulation of observations by calling its [`simulation_observations()`](https://py.api.tudat.space/en/latest/estimation.html#tudatpy.numerical_simulation.estimation.simulate_observations) function. Note that to know about the time settings for the individual types of observations, this function makes use of the earlier defined observation simulation settings.
@@ -332,6 +349,7 @@
     estimator.observation_simulators,
     bodies)
 
+
 # <a id='covariance_section'></a>
 
 ## Perform the covariance analysis
@@ -357,6 +375,7 @@
 weights_per_observable = {estimation_setup.observation.one_way_instantaneous_doppler_type: noise_level ** -2}
 covariance_input.set_constant_weight_per_observable(weights_per_observable)
 
+
 ### Propagate the covariance matrix
 """
 Using the just defined inputs, we can ultimately run the computation of our covariance matrix. Printing the resulting formal errors will give us the diagonal entries of the matrix - while the first six entries represent the uncertainties in the (cartesian) initial state, the seventh and eighth are the errors associated with the gravitational parameter of Earth and the aerodynamic drag coefficient, respectively.
@@ -365,9 +384,11 @@
 # Perform the covariance analysis
 covariance_output = estimator.compute_covariance(covariance_input)
 
+
 # Print the covariance matrix
 print(covariance_output.formal_errors)
 
+
 ## Results post-processing
 """
 Finally, to further process the obtained data, one can - exemplary - plot the correlation between the individual estimated parameters, or the behaviour of the formal error over time.
@@ -391,6 +412,7 @@
 plt.tight_layout()
 plt.show()
 
+
 ### Propagated Formal Errors
 """
 initial_covariance = covariance_output.covariance
@@ -420,3 +442,4 @@
 
 plt.tight_layout()
 plt.show()
+

estimation/estimation_dynamical_models.ipynb

@@ -300,12 +300,15 @@
    "source": [
     "## Define the Dynamical Model(s)\n",
     "Note that unlike it has usually been the case so far - be it with examples dealing with propagation or the prior estimation ones - we have always defined a mere single dynamical model. The modular structure of tudat, however, enables us to simulate the observations using a dynamical model that is (theoretically entirely) different from the one used to perform the estimation. Hence, we will now first define the model that will be used during the simulation of observations. In particular, we will consider:\n",
+    "\n",
     "* Gravitational acceleration using a spherical harmonic approximation up to 4th degree and order for Mars.\n",
     "* Gravitational acceleration using a simple point mass model for:\n",
+    "\n",
     "    - Mars' two moons Phobos and Deimos\n",
     "    - Earth\n",
     "    - Jupiter\n",
     "    - The Sun\n",
+    "\n",
     "* Radiation pressure experienced by the spacecraft - shape-wise approximated as a spherical cannonball - due to the Sun."
    ]
   },
@@ -668,9 +671,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "tudatpy-examples",
+   "display_name": "tudat-space",
    "language": "python",
-   "name": "tudatpy-examples"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -682,7 +685,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.16"
+   "version": "3.10.14"
   }
  },
  "nbformat": 4,

estimation/estimation_dynamical_models.py

@@ -42,6 +42,7 @@
 # Retrieve current directory
 current_directory = os.getcwd()
 
+
 ## Simulation Settings
 """
 After having defined the general configuration of our simulation (i.e. importing required `SPICE` kernels, defining start and end epoch of the simulation) we will create the main celestial bodies involved in the simulation (mainly Mars, its two moons, the two neighbouring planets, and the Sun), the spacecraft itself, and its environment interface.
@@ -89,6 +90,7 @@
 # Define central bodies of propagation
 central_bodies = ["Mars"]
 
+
 time2plt = np.arange(simulation_start_epoch, simulation_end_epoch, 60)
 mex2plt = list()
 for epoch in time2plt:
@@ -117,6 +119,7 @@
 plt.tight_layout()
 plt.show()
 
+
 ## Set Up the Observations
 """
 Having set the underlying environment model of the simulated orbit, we can define the observational model. This entails the addition all required ground stations, the definition of the observation links and types, as well as the precise simulation settings.
@@ -139,6 +142,7 @@
     [station_altitude, new_norcia_latitude, new_norcia_longitude],
     element_conversion.geodetic_position_type)
 
+
 ### Define Observation Model Settings
 """
 Within this example - as it is common practice when tracking deep-space missions using the ESTRACK system - Mars Express will not be tracked using a set of one-way signal path, but a n-way one (realised as two-way link ends in this example). For our example at hand this means that the signal travels from Earth to the spacecraft where it gets re-transmitted and subsequently has to travel back to Earth where it is recorded and processed. In particular, we will model two-way range and range-rate (Doppler) observables.
@@ -168,6 +172,7 @@
     one_way_nno_mex_link_definition,
     light_time_correction_settings = [light_time_correction_settings]))
 
+
 ### Define Observation Simulation Settings
 """
 Finally, for each above-defined observation model, we will define the noise of the observation-type and several, general viability criteria. We impose the spacecraft to be at a certain minimum angle (15 degrees) of elevation above the horizon as seen from the ground station, as well as introduce Mars as a body that can potentially occult the line-of-sight between Mars Express and New Norcia (i.e. when the spacecraft dives 'behind' Mars, we will not simulate any observations).
@@ -210,15 +215,19 @@
     viability_settings
 )
 
+
 ## Define the Dynamical Model(s)
 """
 Note that unlike it has usually been the case so far - be it with examples dealing with propagation or the prior estimation ones - we have always defined a mere single dynamical model. The modular structure of tudat, however, enables us to simulate the observations using a dynamical model that is (theoretically entirely) different from the one used to perform the estimation. Hence, we will now first define the model that will be used during the simulation of observations. In particular, we will consider:
+
 * Gravitational acceleration using a spherical harmonic approximation up to 4th degree and order for Mars.
 * Gravitational acceleration using a simple point mass model for:
+
     - Mars' two moons Phobos and Deimos
     - Earth
     - Jupiter
     - The Sun
+
 * Radiation pressure experienced by the spacecraft - shape-wise approximated as a spherical cannonball - due to the Sun.
 """
 
@@ -244,6 +253,7 @@
         propagation_setup.acceleration.cannonball_radiation_pressure()
     ])
 
+
 ### Perform the observations simulation
 """
 However, following the known - trivial - estimation pipeline, the observations are simulated using the `simulation_observations()` function of the respective `Estimator` object. However, to avoid having to create two distinct estimators, we will manually implement a set of observation simulators upfront, before altering the dynamical model and creating the actual estimator.
@@ -303,6 +313,7 @@
     observation_simulators,
     bodies)
 
+
 ### Alter the Dynamical Model for Mars
 """
 We will now re-purpose the previously defined dynamical model of accelerations acting on `MEX` by altering the perceived gravitational acceleration of Mars onto the spacecraft. In particular, we will remove the gravitational pull of both of Mars' moons - Phobos and Deimos - from the acceleration settings of the spacecraft. All remaining settings, however, remain untouched.
@@ -333,6 +344,7 @@
                   integrator_settings=integrator_settings,
                   termination_settings=termination_settings)
 
+
 ## Perform the estimation
 """
 Having altered the dynamical model as well as the propagator settings, we create the `Estimator` object and subsequently set up the inversion of the problem - in particular, one has to define which parameters are to be estimated, could potentially include any a-priori information in the form of an a-priori covariance matrix, and define the weights associated with the individual types of observations.
@@ -371,6 +383,7 @@
                           estimation_setup.observation.one_way_range_type: noise_level_range ** -2}
 estimation_input.set_constant_weight_per_observable(weights_per_observable)
 
+
 ### Estimate the individual parameters
 """
 Finally, the actual estimation can be performed - ideally having reached a sufficient level of convergence, the least squares estimator will have found the most suitable parameters for the problem at hand.
@@ -381,10 +394,12 @@
 # Perform the covariance analysis
 estimation_output = estimator.perform_estimation(estimation_input)
 
+
 # Print the covariance matrix
 print(estimation_output.formal_errors)
 print(truth_parameters - parameters_to_estimate.parameter_vector)
 
+
 ## Post-processing
 """
 Finally, to further illustrate the impact certain differences between the applied dynamical models have, we will first plot the behaviour of the simulated observations over time, as well as show how the discrepancy between our estimated solution and the 'ground truth' builds up over time.
@@ -404,6 +419,7 @@
 plt.tight_layout()
 plt.show()
 
+
 simulator_object = estimation_output.simulation_results_per_iteration[-1]
 state_history = simulator_object.dynamics_results.state_history
 
@@ -427,3 +443,4 @@
 
 plt.tight_layout()
 plt.show()
+