Add updated NRLMSIS Models

[cfleslie]

The latest version of the empirical atmosphere model NRLMSIS (2.0 and 2.1) has made significant changes based of the satellite data available since 2000. It would be great to have this added as an option like the currently implemented nrlmsise00.

Link to model:
https://ccmc.gsfc.nasa.gov/models/NRLMSIS~2.0/

Publications for versions since 00:
https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2020EA001321
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2022JA030896

[DominicDirkx]

Hi, I think that's a great idea! We'd need a wrapper for the Fortran code, but I think this package does exactly this for Python:
https://swxtrec.github.io/pymsis/
Using the custom atmoshere model in Tudat, you could link pymsis to the Tudat encironment

[cfleslie]

Great! Yeah my plan was to use custom_four_dimensional_constant_temperature with pymsis. Unless you can recommend a better way? From the documentation on the atmosphere acceleration model constant temperature doesnt seem to have any effect on calculating drag, is this correct?

Thanks for your help!

[DominicDirkx]

For now, I think this is the best way to address it, without having to add stuff to the C++ layer of tudat. Whether the temperature impacts the aerodynamic force will depend on your definition of the coefficients. For instance, if your coefficients depend on Mach number (or temperature directly), there is a dependency of the force on the temperature. To remedy this, we should add a custom atmosphere model where you can also define a custom temperature (and composition).

When you have the pymsis-Tudat link working, could you post a message? It would be great to add this as a ready-made option in Tudat

[DominicDirkx]

By the way, we're busy working a bunch of extensions to the code and documentation. I'd love to know more about what you're using Tudat for, and what your experience has been with it. If you'd need it, we could also add the more comprehensive atmosphere model. It's been on our todo list, and I don't think it'd be all that time-consuming. Drop me a message at d.dirkx at tudelft.nl if you'd like to get in touch.

[cfleslie]

Here is a snippet for anybody who wants to get it working with pymsis:

import numpy as np
from pymsis import msis
from tudatpy import numerical_simulation
from tudatpy.numerical_simulation import environment_setup
from tudatpy.interface import spice
from tudatpy.astro.time_conversion import datetime_to_python
from tudatpy.astro.time_conversion import date_time_from_epoch

# Define density function with 4 dimensions as input [alt, lat, long, time]
def density_function(altitude, latitude, longitude, time):
    
    # Convert to python datetime from datetime module
    date_time = datetime_to_python(date_time_from_epoch(time))
    
    # Run msis model
    data = msis.run(
        date_time, 
        np.rad2deg(longitude), # Tudat gives rad, msis expects degrees
        np.rad2deg(latitude),  # Tudat gives rad, msis expects degrees
        altitude/1000.0, # Tudat gives m, msis expects km
        version=2.1) # Version of the model (2.1 is latest)
    return data[0][0] # Return only mass density (kg/m3)

# Standard tudat body definition
spice.load_standard_kernels()

bodies_to_create = ["Earth"]
global_frame_origin = "Earth"
global_frame_orientation = "J2000"

body_settings = environment_setup.get_default_body_settings(
    bodies_to_create, global_frame_origin, global_frame_orientation
)

# Apply custom atmosphere model using density function
body_settings.get("Earth").atmosphere_settings = environment_setup.atmosphere.custom_four_dimensional_constant_temperature(
    density_function,
    constant_temperature=1000.0,
)

I tested the density output from tudat against using pymsis manually with the same altitudes, lats, longs and times and they are in agreement. Let me know if you spot any mistakes I have overlooked!

[DominicDirkx]

Awesome, thanks! One small point: the input time when called from the propagation will be in TDB. Which timescale does msis expect? If it's UTC or TAI, doing an approximate conversion would be good.