Large difference of proton fluxes when compared with Spenvis

[leoghizoni]

Wrote a code to propagate an orbit and get the respective fluxes to compare with the results from Spenvis. The electron fluxes match very well with those from Spenvis, but the respective proton fluxes via radbelt are orders of magnitude higher than those from Spenvis.
I can share the code in case interested.

[lpsinger]

Yes please. I'd need a more detailed bug report.

[leoghizoni]

Yes please. I'd need a more detailed bug report.

I just shared a Drive folder with your email. Please let me know if it doesn't work.

[lpsinger]

Would you please post the reproducer in this GitHub issue?

[leoghizoni]

Would you please post the reproducer in this GitHub issue?

Sure. Graphs? Or table with numbers? Or both? Sorry for the ignorance.

[lpsinger]

https://en.wikipedia.org/wiki/Minimal_reproducible_example

[leoghizoni]

Here's a minimal reproducible example. It's not a small script, but it contains all the necessary libraries to run.

In order to run, the script needs the following libraries: numpy, scipy, pandas, matplotlib, beyond, astropy, orekit(10), and naturally radbelt. I'm running it via conda.

For each of the XYZ coordinates, the code retrieves the fluxes per energy (the exact same energies used by Spenvis), which are later averaged for the full orbit (in the same fashion as Spenvis). In order to compare the fluxes with those from Spenvis, I attach here the trapped electron and proton fluxes for two independent cases: 590x590km, 97.8deg orbit; and 1,000x1,000km, 97.8deg orbit.

The script produces three text files:

• orbit.txt: contains the stime stamped XYZ coordinates of 1 orbit, propagated by means of Orekit
• trp_electrons.txt: The averaged electron fluxes for the orbit, both differential and integral
• trp_protons.txt: The averaged proton fluxes for the orbit, both differential and integral

The electron fluxes from radbelt match very well with those from Spenvis, but the proton fluxes from radbelt are orders of magntitude higher than the Spenvis ones.

`import orekit # type: ignore
vm = orekit.initVM()

from os.path import exists
if not exists('orekit-data.zip'):
from orekit.pyhelpers import download_orekit_data_curdir # type: ignore
download_orekit_data_curdir()

from orekit.pyhelpers import setup_orekit_curdir # type: ignore
setup_orekit_curdir()

from org.orekit.frames import FramesFactory # type: ignore
from org.orekit.bodies import OneAxisEllipsoid # type: ignore
from org.orekit.utils import Constants # type: ignore
from org.orekit.utils import IERSConventions # type: ignore
from org.orekit.time import TimeScalesFactory # type: ignore
from org.orekit.time import AbsoluteDate # type: ignore

import numpy as np
from math import radians
from datetime import datetime

import astropy.constants as const
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import EarthLocation

import matplotlib.pyplot as plt
from scipy.integrate import cumulative_trapezoid

import pandas as pd
from radbelt import get_flux

utc = TimeScalesFactory.getUTC()
gcrf = FramesFactory.getGCRF()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
eci = gcrf # Inertial frame
ecef = itrf # Non-inertial frame

earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ecef)

eEnergies = np.array([4.0e-02,1.0e-01,2.0e-01,3.0e-01,4.0e-01,5.0e-01,6.0e-01,7.0e-01,8.0e-01,
1.00e+00,1.25e+00,1.50e+00,1.75e+00,2.00e+00,2.25e+00,2.50e+00,2.75e+00,
3.00e+00,3.25e+00,3.50e+00,3.75e+00,4.00e+00,4.25e+00,4.50e+00,4.75e+00,
5.00e+00,5.50e+00,6.00e+00,6.50e+00,7.00e+00])

pEnergies = np.array([1.0e-01,1.5e-01,2.0e-01,3.0e-01,4.0e-01,5.0e-01,6.0e-01,7.0e-01,
1.0e+00,1.5e+00,2.0e+00,3.0e+00,4.0e+00,5.0e+00,6.0e+00,7.0e+00,
1.0e+01,1.5e+01,2.0e+01,3.0e+01,4.0e+01,5.0e+01,6.0e+01,7.0e+01,
1.0e+02,1.5e+02,2.0e+02,3.0e+02,4.0e+02])

class KeplerPropagator:
def init(self, wrkDir, spacecraft, orbit, attitude, epoch):
self.wrkDir = wrkDir
self.spacecraft = spacecraft
self.orbit = orbit
self.attitude = attitude
self.epoch = epoch

def analyticalPropagatorBuilder(self, duration,
                                anomaly='TRUE'):
    from org.orekit.orbits import PositionAngle # type: ignore
    from org.orekit.orbits import KeplerianOrbit # type: ignore
    from org.orekit.orbits import EquinoctialOrbit # type: ignore
    from org.orekit.propagation.analytical import KeplerianPropagator # type: ignore
    
    startTime = AbsoluteDate(self.epoch.year, self.epoch.month, self.epoch.day,
                             self.epoch.hour, self.epoch.minute, float(self.epoch.second), utc)
    endTime = startTime.shiftedBy(duration)
    
    # Orbit construction as Keplerian
    if anomaly == 'TRUE':
        initialOrbit = KeplerianOrbit(self.orbit.sma(), self.orbit.eccentricity(),
                                      self.orbit.inclination, self.orbit.aop,
                                      self.orbit.raan, self.orbit.lv,PositionAngle.TRUE,
                                      eci, startTime, Constants.WGS84_EARTH_MU)
    elif anomaly == 'MEAN':
        initialOrbit = KeplerianOrbit(self.orbit.sma(), self.orbit.eccentricity(),
                                      self.orbit.inclination, self.orbit.aop,
                                      self.orbit.raan, self.orbit.lv,PositionAngle.MEAN,
                                      eci, startTime, Constants.WGS84_EARTH_MU)
    else:
        raise ValueError('Unrecognized anomaly type')
    
    initialOrbit = EquinoctialOrbit(initialOrbit)
    propagator = KeplerianPropagator(
        initialOrbit,
        self.attitude,
        Constants.WGS84_EARTH_MU,
        self.spacecraft.mass)
    
    return propagator, startTime, endTime

def propagate(self, propagator, startTime, endTime, step,
              printPV = True, plot=False):  
    from orekit.pyhelpers import absolutedate_to_datetime # type: ignore

    file = open(self.wrkDir+"orbit.txt", "w")
    file.write("Date X[km] Y[km] Z[km]\n")

    while(startTime.compareTo(endTime) <= 0.0):
        
        stateECI = propagator.propagate(startTime)
        pvPosition = stateECI.getPVCoordinates().getPosition()

        date = absolutedate_to_datetime(startTime)
        x = pvPosition.getX()*1e-3
        y = pvPosition.getY()*1e-3
        z = pvPosition.getZ()*1e-3

        file.write(f"{date},{x},{y},{z}\n")
                            
        if printPV:
            print('{} {} {} {}'.format(date, x, y, z))
        startTime = startTime.shiftedBy(step)

    file.close()

class KeplerOrbit:
def init(self, apogee, perigee, inclination, aop, raan, lv):
self.apogee = apogee1000
self.perigee = perigee1000
self.inclination = radians(inclination)
self.aop = radians(aop)
self.raan = radians(raan)
self.lv = radians(lv)

def sma(self):
	return (self.perigee + self.apogee + 2*Constants.WGS84_EARTH_EQUATORIAL_RADIUS)/2.0

def eccentricity(self):
	return abs((self.apogee - self.perigee)/(self.apogee + self.perigee + 2*Constants.WGS84_EARTH_EQUATORIAL_RADIUS))

def meanAnomaly(self):
	return mean_anomaly_from_true(self.eccentricity(),self.lv)

def meanAltitude(self):
	return np.mean([self.perigee, self.apogee])

class Spacecraft:
def init(self, name, mass, cross_section, cd, cr, identification, norad_id):
self.name = name
self.mass = mass
self.cross_section = cross_section
self.cd = cd
self.cr = cr
self.identification = identification
self.norad_id = norad_id

def mod(x, y):
"""Return the modulus after division of x by y.

Python's x % y is best suited for integers, and math.mod returns with the
sign of x.

This function is modelled after Matlab's mod() function.
"""

from math import isinf, isnan, floor

if isnan(x) or isnan(y):
    return float('NaN')

if isinf(x):
    raise ValueError('math domain error')

if isinf(y):
    return x

if y == 0:
    return x

n = floor(x/y)
return x - n*y

def eccentric_anomaly_from_true(e, f):
"""Convert true anomaly to eccentric anomaly."""
E = np.arctan2(np.sqrt(1 - e ** 2) * np.sin(f), e + np.cos(f))
E = mod(E, 2 * np.pi)
return E

def mean_anomaly_from_true(e, f):
"""Convert true anomaly to mean anomaly."""
E = eccentric_anomaly_from_true(e, f)
return E - e * np.sin(E)

def nadirPointing(frame, body):
from org.orekit.attitudes import NadirPointing # type: ignore
return NadirPointing(frame, body)

def orbit_period(ra, rp):
"""Computes the period of 1 orbit"""

sma = (ra*1.e3 + rp*1.e3 + 2*const.R_earth.value)/2.0
period = float(2*np.pi*np.sqrt(sma**3/const.GM_earth.value))

return period

class AEP8:

def __init__(self, wrkDir, orbit):
    self.wkrDir = wrkDir
    self.orbit = orbit

    orbit_file = pd.read_csv(wrkDir+'orbit.txt', delimiter=',', header=1, engine='python')
    orbit_file.columns = ['date', 'x', 'y', 'z']

    self.date = orbit_file['date']
    self.x = orbit_file['x']
    self.y = orbit_file['y']
    self.z = orbit_file['z']

def _ae8(self, case='max', plot=False):

    fluxes = np.zeros((len(self.x), len(eEnergies)))

    i=0
    for x, y, z, date in zip(self.x, self.y, self.z, self.date):

        coords = EarthLocation(x*u.km, y*u.km, z*u.km)

        date = datetime.strptime(date, "%Y-%m-%d %H:%M:%S.%f")
        date = Time(date.strftime("%Y-%m-%d"))

        j=0
        for e in eEnergies:
            fluxes[i][j] = get_flux(coords, date, e*u.MeV, 'e', case).value
            j+=1
        
        i+=1

    differential_flux = np.mean(fluxes, axis=0)
    integral_flux = np.abs(cumulative_trapezoid(differential_flux[::-1],
                                                eEnergies[::-1],
                                                initial=0)[::-1])

    eFile = open(self.wkrDir + "trp_electrons.txt", "w")
    eFile.write("Energy[MeV] IntegralFlux[/cm2/s] DifferentialFlux[/cm2/MeV/s]\n")
    for e, i, f in zip(eEnergies, integral_flux, differential_flux):
        eFile.write(f"{e} {i} {f}\n")
    eFile.close()

    if plot:

        apogee = self.orbit.apogee*1.0e-03
        perigee = self.orbit.perigee*1.0e-03
        inclination = np.rad2deg(self.orbit.inclination)

        fig, ax1 = plt.subplots()
        ax2 = ax1.twinx()
        plot1, = ax1.loglog(eEnergies, differential_flux, '#1f77b4',
                            linewidth=2,
                            label='Differential Flux')
        plot2, = ax2.loglog(eEnergies, integral_flux, '#ff7f0e',
                            linewidth=2,
                            label='Integral Flux')
        ax1.set_xlabel('Energy [MeV]')
        ax1.set_ylabel('Differential Flux [/cm2/MeV/s]')
        ax2.set_ylabel('Integral Flux [/cm2/s]')
        plt.title(f"Apogee = {apogee}km, Perigee = {perigee}km, Inclination = {inclination} deg",
                  fontsize=10)
        plt.suptitle('Trapped Electrons Fluxes')
        plt.grid(color='gray',linestyle='-',linewidth=0.5,which='both')

        # Explicitly enable minor ticks for both axes
        plt.minorticks_on()

        # Optional: Customize the grid and ticks further
        plt.gca().xaxis.set_minor_locator(plt.LogLocator(base=10.0, subs=np.arange(1.0, 10.0) * 0.1, numticks=10))
        plt.gca().yaxis.set_minor_locator(plt.LogLocator(base=10.0, subs=np.arange(1.0, 10.0) * 0.1, numticks=10))

        # Combine legends
        lines = [plot1, plot2]
        labels = [line.get_label() for line in lines]
        ax1.legend(lines, labels)

        plt.show()

def _ap8(self, case='min', plot=False):

    fluxes = np.zeros((len(self.x), len(pEnergies)))

    i=0
    for x, y, z, date in zip(self.x, self.y, self.z, self.date):

        coords = EarthLocation(x*u.km, y*u.km, z*u.km)

        date = datetime.strptime(date, "%Y-%m-%d %H:%M:%S.%f")
        date = Time(date.strftime("%Y-%m-%d"))

        j=0
        for e in pEnergies:
            fluxes[i][j] = get_flux(coords, date, e*u.MeV, 'p', case).value
            j+=1
        
        i+=1

    differential_flux = np.mean(fluxes, axis=0)
    integral_flux = np.abs(cumulative_trapezoid(differential_flux[::-1], pEnergies[::-1], initial=0)[::-1])

    pFile = open(self.wkrDir + "trp_protons.txt", "w")
    pFile.write("Energy[MeV] IntegralFlux[/cm2/s] DifferentialFlux[/cm2/MeV/s]\n")
    for e, i, f in zip(pEnergies, integral_flux, differential_flux):
        pFile.write(f"{e} {i} {f}\n")
    pFile.close()

    if plot:

        apogee = self.orbit.apogee*1.0e-03
        perigee = self.orbit.perigee*1.0e-03
        inclination = np.rad2deg(self.orbit.inclination)

        fig, ax1 = plt.subplots()
        ax2 = ax1.twinx()
        plot1, = ax1.loglog(pEnergies, differential_flux, '#1f77b4', linewidth=2, label='Differential Flux')
        plot2, = ax2.loglog(pEnergies, integral_flux, '#ff7f0e', linewidth=2, label='Integral Flux')
        ax1.set_xlabel('Energy [MeV]')
        ax1.set_ylabel('Differential Flux [/cm2/MeV/s]')
        ax2.set_ylabel('Integral Flux [/cm2/s]')
        plt.title(f"Apogee = {apogee}km, Perigee = {perigee}km, Inclination = {inclination} deg",
                  fontsize=10)
        plt.suptitle('Trapped Protons Fluxes')
        plt.grid(color='gray',linestyle='-',linewidth=0.5,which='both')

        # Explicitly enable minor ticks for both axes
        plt.minorticks_on()

        # Optional: Customize the grid and ticks further
        plt.gca().xaxis.set_minor_locator(plt.LogLocator(base=10.0, subs=np.arange(1.0, 10.0) * 0.1, numticks=10))
        plt.gca().yaxis.set_minor_locator(plt.LogLocator(base=10.0, subs=np.arange(1.0, 10.0) * 0.1, numticks=10))

        # Combine legends
        lines = [plot1, plot2]
        labels = [line.get_label() for line in lines]
        ax1.legend(lines, labels)

        plt.show()

def main():

wrkDir = ''
sat_name = 'Project'

'''
Spacecraft Definition
'''
mass = 161.0 # [kg]
cross_section = 0.6*0.3 + 1.6*1.0 # [m2]
cd = 2.2
cr = 1.8
identification = 1111
norad_id = 41752
spacecraft = Spacecraft(sat_name, mass, cross_section, cd, cr, identification, norad_id)
attitude = nadirPointing(eci, earth)

'''
Keplerian Orbit
'''
epoch = datetime(2026, 1, 1, 0, 0, 0)
ra = 590.0                   # apogee
rp = 590.0                   # perigee
i =  97.8#ssoInclination(ra, rp)  # inclination
omega = 0.0                  # argument of perigee
raan = 0.0#localTime2RAAN(epoch) # right ascension of the ascending node
lv = 0.0#180.0                   # true anomaly
orbit = KeplerOrbit(ra, rp, i, omega, raan, lv)

duration = orbit_period(ra, rp)
step = duration/100.0

'''
Keplerian Propagator
'''
keplerianProp = KeplerPropagator(wrkDir, spacecraft, orbit, attitude, epoch)
propagator, startTime, endTime = keplerianProp.analyticalPropagatorBuilder(duration, anomaly='TRUE')

keplerianProp.propagate(propagator, startTime, endTime, step,
                        printPV=True)

model = AEP8(wrkDir, orbit)
eFluxes = model._ae8(plot=True)
pFluxes = model._ap8(plot=True)

if name == 'main':
main()`

spenvis_590km_sso_trpe.txt
spenvis_590km_sso_trpp.txt
spenvis_1000km_sso_trpe.txt
spenvis_1000km_sso_trpp.txt

[lpsinger]

Would you please turn this into a minimal example? For example, if the details of the orbit are not necessary, can you provide an example that compares the fluxes at a fixed time and location?