Working with space weather data

Added in version 7.3.0.

heyoka.py provides the ability to access space weather (SW) data and use it in the expression system to represent SW indices as time-varying quantities. Space weather support enables accurate modeling of atmospheric density variations, capturing effects such as solar activity and geomagnetic storms that influence atmospheric drag on Earth-orbiting spacecraft (especially in low Earth orbit).

Accessing the raw data

The raw SW data is encapsulated in the sw_data class. Let us explore it a bit.

import heyoka as hy

# Default-construct an sw_data instance.
data = hy.sw_data()

The raw data is available via the table property, which returns a structured NumPy array containing a time series of the SW data:

data.table

array([(36113., 16.5, 264.63333333, 267.  ),
       (36114., 15.5, 257.09166667, 267.75),
       (36115., 15.5, 258.10416667, 268.45), ...,
       (61271., 13.5, 147.21666667, 133.  ),
       (61272., 10. , 151.31666667, 133.15),
       (61273.,  8. , 153.03333333, 133.25)],
      shape=(25161,), dtype={'names': ['mjd', 'Ap_avg', 'f107', 'f107a_center81'], 'formats': ['<f8', '<f8', '<f8', '<f8'], 'offsets': [0, 8, 16, 24], 'itemsize': 32, 'aligned': True})

The dtype of the array is sw_data_row, which contains:

• the UTC modified Julian date (mjd),

• the 24-hour running average of the Ap index centred on the date (Ap_avg),

• the observed 10.7-cm solar radio flux F10.7 (f107),

• the 81-day running average of the observed F10.7 centred on the date (f107a_center81).

Updating the SW data

A default-constructed sw_data instance uses a builtin SW dataset from CelesTrak. This dataset is comprehensive, including both historical data dating back to the 50s and predictions for the near future.

Although the builtin dataset is updated at every new heyoka.py release, it is likely to be outdated, at least for operational uses. For this reason, the sw_data class provides static factory methods to construct sw_data instances from up-to-date datasets downloaded from the internet.

As an example, we can use fetch_latest_celestrak() to download the dataset covering the last 5 years before present from celestrak:

updated_data = hy.sw_data.fetch_latest_celestrak()

Please see the documentation of sw_data for more detailed information.

Using the SW data in the expression system

heyoka.py’s expression system includes several SW indices implemented as time-dependent functions. These SW functions can then be employed, e.g., to formulate the dynamics of Earth-orbiting spacecraft.

The SW functions are all implemented as piecewise linear functions of the input time coordinate, and currently they include:

• the 24-hour running average of the Ap index Ap_avg() and its derivative Ap_avgp(),

• the observed 10.7-cm solar radio flux f107() and its derivative f107p(),

• the 81-day running average of the observed 10.7-cm solar radio flux f107a_center81() and its derivative f107a_center81p().

All these functions take in input a time expression assumed to represent the number of Julian centuries elapsed since the epoch of J2000 in the terrestrial time (TT) timescale. They also accept as a second (optional) argument the sw_data dataset to be used in the computation - if not provided, the builtin SW dataset will be used.

Here is an example of computation and visualisation of the time evolution of the observed 10.7-cm solar radio flux and its running average using heyoka.py’s compiled functions:

import numpy as np

# Define the computation timespan: 22 years following J2000.
tspan = np.linspace(0, 0.22, 10000)

# Introduce a symbolic variable for the representation of time.
tm = hy.make_vars("tm")

# Construct a compiled function for the computation of the 10.7-cm solar radio flux.
cf = hy.cfunc(
    [hy.model.f107(time_expr=tm), hy.model.f107a_center81(time_expr=tm)], [tm]
)

# Compute the values of f107 and its average over tspan.
f107, f107a = cf(inputs=[tspan])

%matplotlib inline
from matplotlib.pylab import plt

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot()

# Convert the timespan to days for visualisation.
cy_to_days = 36525

ax.plot(tspan * cy_to_days, f107, label="F10.7 flux")
ax.plot(tspan * cy_to_days, f107a, label="81-days average of F10.7")
ax.set_xlabel("Days since J2000")
ax.set_ylabel("Solar radio flux (sfu)")
ax.legend();

Using custom SW datasets

Added in version 7.12.0.

In addition to the builtin and downloadable data sets, it is also possible to manually construct SW datasets from user-supplied data.

First, we generate the raw data table as a structured NumPy array (using some invented values for the MJDs and the SW indices):

tbl = np.array(
    [
        (36113.0, 16.5, 264.63333333, 267.0),
        (36114.0, 15.5, 257.09166667, 267.75),
        (36115.0, 15.5, 258.10416667, 268.45),
    ],
    dtype=hy.sw_data_row,
)

Then, we need to pick an identifier and a timestamp for our custom SW dataset. heyoka.py internally uses timestamp and identifier to discriminate between SW datasets: two datasets with the same timestamp, identifier and number of data rows are considered identical.

Both timestamp and identifier can be arbitrary strings, as long as they are non-empty and they do not contain the - or . characters:

timestamp = "20260611"
identifier = "my_custom_data"

Now we can proceed to the construction of our custom SW dataset:

custom_data = hy.sw_data(data=tbl, timestamp=timestamp, identifier=identifier)

The custom dataset is now ready to be used in the expression system:

# Construct an expression for the solar radio flux
# using the custom dataset.
hy.model.f107(sw_data=custom_data)

sw_f107_3_20260611-my_custom_data(t)

As you can see, the timestamp, the identifier and the number of rows have all been mangled into the name of the function.