Computing ionospheric delay from a IONEX file

[AlekLee]

I’ve been trying to compute the ionospheric delay between a specific GPS satellite and a receiver at a specific time. I’m doing this because I want to compare the uncorrected measured pseudorange values to other measurements.

I thought it’d be the easiest to start with delays computed from a IONEX file, so I’ve used this: GlobalIonosphereMapModel (OREKIT 13.1.2 API)

I’ve been using this IONEX file: EMR0OPSFIN_20252880000_01D_01H_GIM.INX (1.6 MB), which I’ve retrieved from here: Earthdata Login .

I load the IONEX file like this, and this runs w/o issues:

from java.io import File
from org.orekit.data import DataContext, DirectoryCrawler
from org.orekit.models.earth.ionosphere import GlobalIonosphereMapModel

# Point to the IONEX file
ionex_dir = File(r'C:\Users\MyName\Downloads\IONEX')
ionex_fname = 'EMR0OPSFIN_20252880000_01D_01H_GIM.INX'
manager = DataContext.getDefault().getDataProvidersManager()
manager.addProvider(DirectoryCrawler(ionex_dir))

# Create the ionosphere model
iono_model = GlobalIonosphereMapModel(ionex_fname)

I then take my receiver’s position, which I know from another source, and construct the needed TopocentricFrame to get this:

# receiver_pv comes from an OEM, which I don't show here for brevity.
epoch = AbsoluteDate(2025, 10, 15, 0, 0, 6.0, TimeScalesFactory.getUTC())
receiver_geo_point = earth.transform(receiver_pv.getPosition(),
                                     j2000, epoch)
receiver_topo_frame = TopocentricFrame(earth, receiver_geo_point, 'RX')

These are initialised correctly:

In [4]: receiver_geo_point
Out[4]: <GeodeticPoint: {lat: -37.8356164374 deg, lon: -84.0953655701 deg, alt: 576,329.4309704629}>
In [5]: receiver_topo_frame
Out[5]: <TopocentricFrame: RX>

I then get the state of the GPS satellite, here PRN09, from an SP3 file:

In [6]: almanac # From SP3, which I've tested against a reference and it's read correctly.
Out[6]: <SP3Ephemeris: org.orekit.files.sp3.SP3Ephemeris@341814d3>

propagator = almanac.getPropagator()
prn09_pv = propagator.getPVCoordinates(epoch, j2000)
prn09_pt = EphemerisHandlers.EphemerisPoint.from_orekit_j2k_pv(prn09_pv)
prn09_state = SpacecraftState(prn09_pt.get_orekit_cartesian_orbit())

This gives a correct state:

In [7]: prn09_state
Out[7]: <SpacecraftState: SpacecraftState{orbit=Cartesian parameters: {P(3751761.187842153, -2.6251963800329804E7, -163328.23606034936), V(2184.7523157592145, 319.63584050937396, -3190.500677653312)}, attitude=org.orekit.attitudes.Attitude@4397ad89, mass=1000.0, additional={}, additionalDot={}}>

Computing the delay should now be easy:

frequency = 1575.42e6  # L1 frequency in Hz
delay_prn09 = iono_model.pathDelay(prn09_state, receiver_topo_frame,
                                   epoch, frequency, None)

However, the delay I get is zero:

In [8]: delay_prn09
Out[8]: 0.0

How is this possible? What am I doing wrong? The elevation of PRN09 is pretty high, and the IONEX file time span also should be OK (I get the same error even further away from the start of the IONEX file epoch range):

In [12]: receiver_topo_frame.getElevation(prn09_state.getPVCoordinates().getPosition(), prn09_state.getFrame(), epoch)
Out[12]: 0.6177845133215872

I ‘ve noticed this sentence in the docs:

Note that this model pathDelay methods requires the topocentric frame to lie on a OneAxisEllipsoid body shape, because the single layer on which pierce point is computed must be an ellipsoidal shape at some altitude.

So, I’ve also tried it with a receiver located on the Earth surface, and that seems to work:

In [9]: land_geo_point = GeodeticPoint(35.7, 139.8, 0.0)
In [10]: land_topo_frame = TopocentricFrame(earth, land_geo_point , 'Tokyo')
In [11]: iono_model.pathDelay(prn09_state, land_topo_frame , epoch, frequency, None)
Out[11]: 3.9532514475323284 # Makes sense.

Is there a way to use the IONEX files for non-earthbound (i.e. spaceborne) GPS receivers? If not, what other options exist to compute the ionospheric delay?

I’ve seen the GRAPHIC method was implemented in Orekit, and I could try to use this (I only have single-frequency measurements): GRAPHICCombination (OREKIT 13.1.2 API) , but I can’t see any examples of how to use it. The documenation doesn’t even mention units that the methods take… Are there any examples of how to use this and similar measurement combinations? I’ve only found this, but it’s a bit sant in details: Ionosphere free LC - #7 by bcazabonne . And anyway, this doesn’t help me compute the ionospheric delay, but the corrected measurement, it seems.

[luc]

It is possible to compute ionospheric delays for spaceborne receivers, but for now you still have to create a temporary topocentric frame from a geodetic point. This is not a huge overhead, though, so it is not a big deal. Note that @baubin is changing the measurement hierarchy anf API, and it does involve dealing with ionospheric delays for spaceborne receivers.

Also note that for a spaceborne receiver, the altitude along the signal path is not always structly descending, you may have a lowest point and a decreasing part from emitter to lowest point then an increasing part from lowest point to receiver. Two ionospheric delays should be computed setting the topocentric frame at this lowest point (one towards emitter and one towards receiver) and then the two delays should be added.

[AlekLee]

@luc

ntric frame from a geodetic point. This is not a huge overhead, though, so it is not a big deal. Note that @baubin is changing the measurement hierarchy anf API, and it does involve dealing with ionospheric delays for spaceborne receivers.

Also note that for a spaceborne receiver, the altitude along the signal path is not always structly descending, you may have a lowest point and a decreasing pa

Thanks! It’s good to hear there’s some work being done on this.

That’s exactly what I did (crearted a topocentric frame at the spaceborne receiver location), but the delay I got was zero. Did I do something wrong? Does it have to do with the explanation about the adding of the two delays? I’m sorry, but I didn’t quite get that. I tried to draw it to understand it better, but I sadly still don’t:

At the end of the day, the receiver records the pseudorange at the reception epoch, which is light-time delay after the transmission epoch. In the example above, I took the toopocentric position of the receiver at the reception epoch, and computed the delay then. Where else should I place the other topocentric frame?

Also, from your response I take it there’s no guide/tutorial/more docs on how to use the GRAPHIC algorithm? Just to be sure.

[baubin]

1. Yes I am going to rewrite the ionospheric delay routines so they automatically handle spaceborne receivers, but I’m afraid that will take some time. The MR is not small. Just a fair warning.
2. Just a check, but do you know if the sender of the signal was above the horizon when you ran the example check? If the signal sender is below the horizon of the receiver, then the code automatically gives a delay of zero. What might be tripping you up here is that if your receiver is in space, then a satellite that is actually in view of your receiver, but has a lower altitude than your receiver, would have a < 0 deg elevation from the POV of your artificial ground point.

[baubin]

@AlekLee Here are a couple of diagrams to illustrate the points. First one shows what Luc means about the lowest-altitude point between sender and receiver not necessarily being either sender OR receiver:

what Luc is saying here is that in order to properly calculate the delay between sender and receiver using the model, we’d have to calculate two legs: the first leg between the GPS sat and the lowest-altitude point, and the second between the lowest-altitude point and the receiver. Once my MR is done, this will happen automatically however and no one will need to worry about it anymore.

Second one shows what I mean about the GPS receiver possibly being “below your elevation”

[baubin]

I also just noticed that the altitude of your geodetic point was 576 km. That might just be too high for some atmospheric models to deal with. That’s a maybe though. I can’t be sure without actually digging into which ionospheric model you’re using.

[AlekLee]

@baubin

at might be tripping you up here is that if your receiver is in space, then a satellite that is actually in view of your receiver, but has a lower altitude than your receiver, would have a < 0 deg el

Yes, I actually know the receiver started tracking PRN09 at this time. Here’s a figure showing the geometry. I computed the elevation as 12 degrees (beginning of tracking). Thanks for the suggestion, though. I did have an issue with this earlier

[AlekLee]

Thanks, this was the issue here, apparently. I’ve checked, and the model works for the sub-satellite point but then breaks down at about 455 km, even though the elevations is still comfortably above 0 deg:

[AlekLee]

I’ve been learning about pseudorange and GNSS in Orekit lately, and haven’t seen many comprehensive examples out there. Here’s my script that I used with this question, hoping it’ll help someone. It works with the IONEX and SP3 files that I attached to the original question. I haven’t yet implemented the “piecewise” delay calculation that @luc and @baubin mentioned above, but hopefully it’s a start.

import os
import math
import numpy
import logging
import orekit_jpype as orekit
import matplotlib.pyplot as plt

frequency = 1575.42e6  # L1 frequency in Hz.

psr_measured = 23972370.14 # m

#%% Configure the application logger.
logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)
handler = logging.StreamHandler()
handler.setLevel(logging.INFO)
formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
handler.setFormatter(formatter)
logger.addHandler(handler)

#%% Set-up Orekit, import its modules and prepare the needed coordinate frames.
orekit.initVM()

from orekit_jpype.pyhelpers import absolutedate_to_datetime, download_orekit_data_curdir, setup_orekit_curdir

if not os.path.isfile('orekit-data.zip'):
    download_orekit_data_curdir()
setup_orekit_curdir()

from org.orekit.utils import Constants, IERSConventions, PVCoordinates, AbsolutePVCoordinates
from org.orekit.models.earth.ionosphere import GlobalIonosphereMapModel
from org.orekit.data import DataSource, DataContext, DirectoryCrawler
from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint
from org.orekit.frames import FramesFactory, TopocentricFrame
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.propagation import SpacecraftState
from org.orekit.files import sp3
from java.io import File
from org.hipparchus.geometry.euclidean.threed import Vector3D

j2000 = FramesFactory.getEME2000()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                         Constants.WGS84_EARTH_FLATTENING, itrf)

#%% Read the SP3 file.
# 5-minute final, fitted after applying a Helmert transformation from Natural Resources Canada.
# Source: https://cddis.nasa.gov/archive/gnss/products/2388/
SP3_FNAME = 'EMR0OPSFIN_20252880000_01D_05M_ORB.SP3'

sp3Parser = sp3.SP3Parser()
sp3DataSource = DataSource(os.path.join(os.getcwd(), SP3_FNAME))
almanacsSP3 = sp3Parser.parse(sp3DataSource).getSatellites()

almanacs_used_sp3 = [almanacsSP3.get(key) for key in almanacsSP3.keySet()]

# Same for all PRNs, actually, but just to be safe.
assert all([absolutedate_to_datetime(al.getStart()) ==
            absolutedate_to_datetime(almanacs_used_sp3[0].getStart())
            for al in almanacs_used_sp3])

#%% Find GPS satellite's position.
# Reception epoch of interest.
epoch = AbsoluteDate(2025, 10, 15, 0, 0, 4.05, TimeScalesFactory.getUTC())

almanac = almanacs_used_sp3[8]
assert almanac.getId() == 'G09' # Just because.
propagator = almanac.getPropagator()

# Find transmission epoch, i.e. account for light-time delay.
light_time_delay = psr_measured / Constants.SPEED_OF_LIGHT
epoch_transmission = epoch.shiftedBy(-1*light_time_delay)

# This is where PRN09 was when it sent the signal.
prn09_pv = propagator.getPVCoordinates(epoch_transmission, j2000)
prn09_apv = AbsolutePVCoordinates(j2000, epoch, prn09_pv.getPosition(), prn09_pv.getVelocity())
prn09_state = SpacecraftState(prn09_apv)

#%% Read the IONEX file.
# Source: https://cddis.nasa.gov/archive/gnss/products/ionex/2025/288/
# Point to the IONEX file
ionex_dir = File(os.getcwd())
ionex_fname = 'EMR0OPSFIN_20252880000_01D_01H_GIM.INX'
manager = DataContext.getDefault().getDataProvidersManager()
manager.addProvider(DirectoryCrawler(ionex_dir))

# Create the ionosphere model
iono_model = GlobalIonosphereMapModel(ionex_fname)

#%% Compute ionospheric delay at a particular terrestrial location at the same epoch.
# More or less below the spaceborne receiver.
ground_geo_point = GeodeticPoint(math.radians(-59.4), math.radians(-76.4), 0.0)
ground_topo_frame = TopocentricFrame(earth, ground_geo_point, 'Ground')
delay_ground = iono_model.pathDelay(prn09_state, ground_topo_frame, epoch, frequency, None)

print(f"Ionospheric delay to ground receiver is {delay_ground} metres")

# Check elevation of the PRN09.
elevation_from_ground = ground_topo_frame.getElevation(
    prn09_state.getPVCoordinates().getPosition(), prn09_state.getFrame(), epoch)

print(f"PRN09 elevation from ground receiver is {math.degrees(elevation_from_ground)} deg")

#%% Compute the ionospheric delay at a particular spaceborne location at a particular epoch.
# Receiver position - taken from an OEM.
receiver_pv = PVCoordinates(Vector3D(2123944.7, -2834143.8, -5968585.3),
                            Vector3D(1059.1, -4650.3, 5882.4))
receiver_geo_point = earth.transform(receiver_pv.getPosition(), j2000, epoch)
receiver_topo_frame = TopocentricFrame(earth, receiver_geo_point, 'RX')

#TODO Implement the method suggested by Luc and Baubin to compute the delay along the "two legs" https://forum.orekit.org/t/computing-ionospheric-delay-from-a-ionex-file/5137/5

# Compute ionospheric delay.
delay_receiver = iono_model.pathDelay(prn09_state, receiver_topo_frame,
                                   epoch, frequency, None)

# The spaceborne receiver is at 576 km altitude, which is too high for this
# model, apparently. Run the next cell to convince yourself about it.
print(f"Ionospheric delay to spaceborne receiver at {receiver_geo_point.getAltitude()/1000:.1f} km is {delay_receiver} metres")

# Check elevation.
elevation_from_receiver = receiver_topo_frame.getElevation(
    prn09_state.getPVCoordinates().getPosition(), prn09_state.getFrame(), epoch)

print(f"PRN09 elevation from spaceborne receiver at {receiver_geo_point.getAltitude()/1000:.1f} km is {math.degrees(elevation_from_receiver)} deg")

#%% Compare measured to simulated PSR.
# Find the geometric range between the receiver and the GPS satellite.
psr_sim = float(receiver_pv.getPosition().distance(prn09_pv.getPosition()))

# Get the PRN09 clock offset form the SP3 file and correct for it.
clock_model = almanac.extractClockModel()
offset = clock_model.getOffset(epoch).getOffset()
psr_sim = psr_sim - Constants.SPEED_OF_LIGHT * offset

# Here, apply further corrections - ionosphere, troposphere etc.

print(f'Measured PSR={psr_measured:.1f} m, simulated={psr_sim:.1f} m. Difference is {psr_measured-psr_sim:.1f} m.')

#%% Plot all of this in 3D to make sure it all makes sense.
fig = plt.figure()
ax3d = fig.add_subplot(projection='3d')

pos = ground_topo_frame.getCartesianPoint()
ax3d.scatter(pos.getX(), pos.getY(), pos.getZ(), c='k', label='Ground RX')

pos = j2000.getTransformTo(itrf, epoch).transformPVCoordinates(receiver_pv).getPosition()
ax3d.scatter(pos.getX(), pos.getY(), pos.getZ(), c='tab:blue', label='Space RX')

pos = j2000.getTransformTo(itrf, epoch).transformPVCoordinates(prn09_pv).getPosition()
ax3d.scatter(pos.getX(), pos.getY(), pos.getZ(), c='tab:red', label='PRN09 TX')

noEarthPoints = 100
RE = Constants.EGM96_EARTH_EQUATORIAL_RADIUS
uEarth = numpy.linspace(0.0, 2.0 * numpy.pi, noEarthPoints)
vEarth = numpy.linspace(0.0, numpy.pi, noEarthPoints)
xEarth = RE * numpy.outer(numpy.cos(uEarth), numpy.sin(vEarth))
yEarth = RE * numpy.outer(numpy.sin(uEarth), numpy.sin(vEarth))
zEarth = RE * numpy.outer(numpy.ones_like(uEarth), numpy.cos(vEarth))
ax3d.plot_wireframe(xEarth, yEarth, zEarth, rstride=4, cstride=4, color='tab:grey',
                    alpha=0.5, linewidths=0.1)

ax3d.set_xlim(-4*RE, 4*RE)
ax3d.set_ylim(-4*RE, 4*RE)
ax3d.set_zlim(-4*RE, 4*RE)
ax3d.set_xlabel('X ECEF (m)')
ax3d.set_ylabel('Y ECEF (m)')
ax3d.set_zlabel('Z ECEF (m)')
ax3d.view_init(elev=33, azim=-51, roll=0)
ax3d.legend()
fig.show()

#%% Compute ionospheric delay at several intermediate altitudes and plot.
alts = list(range(0, 550, 5))
delays = []
elevations = []

for alt in alts:
    receiver2_geo_point = GeodeticPoint(math.radians(-59.4), math.radians(-76.4),
                                        alt*1000.0)
    receiver2_topo_frame = TopocentricFrame(earth, receiver2_geo_point, 'Space')
    delay_receiver2 = iono_model.pathDelay(prn09_state, receiver2_topo_frame, epoch, frequency, None)
    
    # print(f"Ionospheric delay to spaceborne receiver 2 at {alt} km is {delay_receiver2} metres")
    
    # Check elevation.
    elevation_from_receiver2 = receiver2_topo_frame.getElevation(
        prn09_state.getPVCoordinates().getPosition(), prn09_state.getFrame(), epoch)
    
    el = math.degrees(elevation_from_receiver2)
    # print(f"PRN09 elevation from spaceborne receiver 2 at {alt} km is {el} deg")
    
    delays.append(delay_receiver2)
    elevations.append(el)

# See how the delay and elevation vary with altitude.
fig, ax1 = plt.subplots()
ax1.set_xlabel('Altitude (km)')
ax1.set_ylabel('Delay (m)', color='tab:red')
ax1.plot(alts, delays, color='tab:red')
ax1.tick_params(axis='y', labelcolor='tab:red')

ax2 = ax1.twinx()
ax2.set_ylabel('PRN09 elevation (deg)', color='tab:blue')
ax2.plot(alts, elevations, color='tab:blue')
ax2.tick_params(axis='y', labelcolor='tab:blue')

fig.tight_layout()
fig.show()

[baubin]

So just an FYI, a couple of years ago I compared the model delay values for a spaceborne receiver to the L1/L2 correction calculations, and they were… very different. Most ionospheric models assume your receiver is on the ground, so a lot of them will not work great for spaceborne receivers. If you’re getting two frequencies in your actual data, you’re much better off using the correction equation than a model.

Ionospheric delay calculation

[AlekLee]

Haha, I wish I had L2… Thanks for the suggestion, though. I’ll try out a few other models implemented in Orekit, and then maybe try the GRAPHIC algorithm: https://www.dlr.de/de/rb/medien/publikationen/publikationen/gsoc-dokumente/ast_03396.pdf/@@download/file (original, behind a paywall: https://www.sciencedirect.com/science/article/pii/S0273117708006340).

[luc]

@baubin

Most ionospheric models assume your receiver is on the ground,

Yes, this is the case in particular for models that are based on a single layer for electronic content or that already assume everything has been integrated along the path, taking a mapping function into account.

There are models that do not make this assumption: NeQuick ITU and NeQuick-G (the second one is a variation of the first one, specialized for Galileo using a slightly different solar flux input with a quadratic polynomial and a slightly different integration method). These models integrate the electronic content along the path. The drawback of these models is that they are computationally intensive. One often has to wrap them with some caching to trade accuracy against computing cost. This is what I had to do recently for ionospheric scintillation computation, which uses NeQuick-ITU underneath.

[AlekLee]

@luc

is is the case in particular for models that are based on a single layer for electronic content or that already assume everything has been integrated along the path, taking a mapping function into account.

There are models that do not make this assumption: NeQuick ITU and NeQuick-G (the second one is a variation of the first one, specialized for Galileo using a slightly different solar flux input with a quadratic polynomial and a slightly different integration method).

Thanks a lot for pointing this out! Incidentally, I already implemented the NeQuickItu model in my MWE script above last night (I’ll probably update it with a better, more comprehensive version later). For now, I’ve been taking the F10.7 flux from the CSSI weather file like so:

from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
from org.orekit.models.earth.ionosphere.nequick import NeQuickItu
space_weather = CssiSpaceWeatherData('SpaceWeather-All-v1.2.txt')
f107 = space_weather.getDailyFlux(epoch)
iono_model_neq = NeQuickItu(f107, TimeScalesFactory.getUTC())
# ...
delay_ground_neq = iono_model_neq.pathDelay(gps_tx_state, ground_topo_frame,
                                            epoch, frequency, None)

Do you know if using the daily flux is the correct way to go here? On the navipedia page (link), they mention “daily flux”, but it isn’t verry clear if that’s just a coincidence or an actual recommendation. The actual model is available from ITU here: https://www.itu.int/dms_pubrec/itu-r/rec/p/R-REC-P.531-16-202509-I!!ZIP-E.zip. In the README, they say that:

The sub-models contained in NeQuick2 use monthly average values of solar activity in two forms: twelve month running mean sunspot number R12 and 10.7 cm wavelength solar radio flux F107. The latter is considered to be the primary input parameter. The following relation between R12 and F107 is used.

This would point to not using the daily flux, but the mean flux instead: space_weather.getMeanFlux(epoch). In the Galileo model spec (link), they also mention the 12-month average. Other ITU recommendations (link1, link2) don’t even mention “flux” or “F10.7”, so I’m not sure. The documentation in the source (link) is no more verbose than what’s in the documenation. I’d greatly appreciate some advice here, @luc .

I noticed that the NeQuickGalileo ( NeQuickGalileo (OREKIT 13.1.2 API) ) and NeQuickModel ( NeQuickModel (OREKIT 13.1.2 API) ) use the effective ionisation level az, as defined in the Galileo spec (link). Is there a way to retrive this metric in a similar way to how one can get the Klobuchar parameters?

Here’s how the two models compare for my tow example, in case it helps anyone:

[luc]

Ionospheric effects vary a lot daily, so using a current flux rather than some mean flux seems better to me, but be aware I am really not an expert in ionosphere! I just happen to have implemented the models blindly from the specs without understanding what I did…

In order to get the Az for Galileo, I think it is broadcast by the satellites. You may find it in Rinex navigation messages, as they basically dump the broadcast messages.

[AlekLee]

In case you’re interested or it helps anyone else, it’s pretty easy to get the Galileo broadcast coefficients - they are in the combined multi-system RINEX files, available e.g. here: Earthdata Login

The same file also has GPS coefficients needed for the Klobuchar model. Corresponding legacy RINEX files with only GPS and Glonas are also available here: Earthdata Login

It turns out there is a ~5.6 m diference betwee the Klobuchar and the NeQuick models. And even within the NeQuick models, you can get a spread of almost a metre, depending on the source of the flux data (CSSI daily, mean or the broadcast total ionisation level coeffieicnts). Not easy when trying to do precise OD!

Anyway, thanks a lot for all your help!

[AlekLee]

I’ve looked into this more and built the following algorithm to check if the signal gets closer to Earth than the RX altitude after it’s transmitted from the TX:

# Understand whether the receiver is the closest point to the Earth (lowest altitude)
# that the signal ever reaches. There are many ways to parametrise this. Here,
# we check the Eath-RX-TX angle - if it's more than 90 deg, the ray is moving
# away from the Earth. Another way to look at it is: the height of the triangle
# perpendicular to the TX-RX side doesn't intersect this side between RX and TX.
rx2tx_unit_vec = gps_tx_pos.subtract(receiver_pos).normalize()

if Vector3D.angle(rx2tx_unit_vec, receiver_pos.negate().normalize()) >= math.pi/2:
    # RX location is the lowest altitude the ray ever reaches.
    min_alt_pos = receiver_pos

else:
    # After transmission from TX, the ray reaches the lowest altitude before
    # it reaches the RX. Need to take that into account when computing the delays.
    
    # Find the height of the triangle perpendicular to TX-RX side.
    rx_tx_distance = gps_tx_pos.distance(receiver_pos)
    s = (receiver_pos.getNorm() + gps_tx_pos.getNorm() +
         rx_tx_distance)/2 # Semi-perimeter.
    A = math.sqrt(s*(s-receiver_pos.getNorm()
                     )*(s-gps_tx_pos.getNorm())*(s-rx_tx_distance)) # Heron's formula.
    h = 2*A/rx_tx_distance
    
    # Find the distance from RX to where the height intercepts the TX-RX side
    # (minimum altitude point).
    d = math.sqrt(receiver_pos.getNorm()*receiver_pos.getNorm() - h*h)
    min_alt_pos = receiver_pos.add(rx2tx_unit_vec.scalarMultiply(d))

It relies on the Earth-RX-TX angle (pink on the diagram), which is also a proxy of elevation at the RX location. If it’s larger than 90 degrees, the RX-TX side gets to the lowest altitude (black line) behind the reicever, so it doesn’t matter to signal propagation:

I guess the above is the typical scenario if your GNSS antenna is pointing up towards space, and you aren’t likely to track anything close to the horizon. It is possible, though, which gives rise to this case that @luc and @baubin mentioned:

Here’s an excerpt from my script (extension of the above) with two hard-coded reciever_pos es that you can use to reproduce the two above figures:

# The lowest-altitude point along the RX-TX line is at:
# 709298.5682995064, 6770452.134745813, 2988265.404524084
# Use the following hard-coded receiver_pos to see what this looks like if the
# ray crosses the lowest altitude between the RX and the TX:
# receiver_pos = Vector3D(3701731.923046343, 5951127.721921541, -1416451.1254569488)
# Use this position to see the ray never getting closer to the Earth thant he RX:
# receiver_pos = Vector3D(1564279.5267986006, 6536359.44536745, 1729774.967386649)
# Use the following to just visualise what the script is running:
receiver_pos = j2000.getTransformTo(itrf, epoch).transformPVCoordinates(receiver_pv).getPosition()
gps_tx_pos = j2000.getTransformTo(itrf, epoch).transformPVCoordinates(gps_tx_pv).getPosition()

fig = plt.figure()
ax3d = fig.add_subplot(projection='3d', aspect='equal')

# Plot the ground and spaceborne receivers, and the GPS TX.
ground_rx_pos = ground_topo_frame.getCartesianPoint()
ax3d.scatter(ground_rx_pos.getX(), ground_rx_pos.getY(), ground_rx_pos.getZ(),
             c='k', label='Ground RX')

ax3d.scatter(receiver_pos.getX(), receiver_pos.getY(), receiver_pos.getZ(),
             c='tab:blue', label='Space RX')

ax3d.scatter(gps_tx_pos.getX(), gps_tx_pos.getY(), gps_tx_pos.getZ(),
             c='tab:red', label='GPS TX')

# Plot the RX - Earth - TX triangle - use it to find the point along the signal
# path with the lowest altitude.
ax3d.plot([0, receiver_pos.getX()], [0, receiver_pos.getY()], [0, receiver_pos.getZ()],
          c='tab:blue')
ax3d.plot([0, gps_tx_pos.getX()], [0, gps_tx_pos.getY()], [0, gps_tx_pos.getZ()],
          c='tab:red')
ax3d.plot([gps_tx_pos.getX(), receiver_pos.getX()], [gps_tx_pos.getY(), receiver_pos.getY()],
          [gps_tx_pos.getZ(), receiver_pos.getZ()], c='tab:orange')

# Understand whether the receiver is the closest point to the Earth (lowest altitude)
# that the signal ever reaches. There are many ways to parametrise this. Here,
# we check the Eath-RX-TX angle - if it's more than 90 deg, the ray is moving
# away from the Earth. Another way to look at it is: the height of the triangle
# perpendicular to the TX-RX side doesn't intersect this side between RX and TX.
rx2tx_unit_vec = gps_tx_pos.subtract(receiver_pos).normalize()

if Vector3D.angle(rx2tx_unit_vec, receiver_pos.negate().normalize()) >= math.pi/2:
    # RX location is the lowest altitude the ray ever reaches.
    min_alt_pos = receiver_pos

else:
    # After transmission from TX, the ray reaches the lowest altitude before
    # it reaches the RX. Need to take that into account when computing the delays.
    
    # Find the height of the triangle perpendicular to TX-RX side.
    rx_tx_distance = gps_tx_pos.distance(receiver_pos)
    s = (receiver_pos.getNorm() + gps_tx_pos.getNorm() +
         rx_tx_distance)/2 # Semi-perimeter.
    A = math.sqrt(s*(s-receiver_pos.getNorm()
                     )*(s-gps_tx_pos.getNorm())*(s-rx_tx_distance)) # Heron's formula.
    h = 2*A/rx_tx_distance
    
    # Find the distance from RX to where the height intercepts the TX-RX side
    # (minimum altitude point).
    d = math.sqrt(receiver_pos.getNorm()*receiver_pos.getNorm() - h*h)
    min_alt_pos = receiver_pos.add(rx2tx_unit_vec.scalarMultiply(d))

ax3d.scatter(min_alt_pos.getX(), min_alt_pos.getY(), min_alt_pos.getZ(),
             c='tab:green', label='Min. alt.')

# Plot the Earth.
noEarthPoints = 100
RE = Constants.EGM96_EARTH_EQUATORIAL_RADIUS
uEarth = numpy.linspace(0.0, 2.0 * numpy.pi, noEarthPoints)
vEarth = numpy.linspace(0.0, numpy.pi, noEarthPoints)
xEarth = RE * numpy.outer(numpy.cos(uEarth), numpy.sin(vEarth))
yEarth = RE * numpy.outer(numpy.sin(uEarth), numpy.sin(vEarth))
zEarth = RE * numpy.outer(numpy.ones_like(uEarth), numpy.cos(vEarth))
ax3d.plot_wireframe(xEarth, yEarth, zEarth, rstride=4, cstride=4, color='tab:grey',
                    alpha=0.5, linewidths=0.1)

# Misc. plot formatting.
ax3d.set_xlim(-4*RE, 4*RE)
ax3d.set_ylim(-4*RE, 4*RE)
ax3d.set_zlim(-4*RE, 4*RE)
ax3d.set_xlabel('X ECEF (m)')
ax3d.set_ylabel('Y ECEF (m)')
ax3d.set_zlabel('Z ECEF (m)')
ax3d.view_init(elev=27, azim=-17, roll=0)
ax3d.legend()
fig.show()

[luc]

Note that there is a dedicated method to compute the lowest point: OneAxisEllipsoid.java, line 846.

There is however a catch. Ionospheric models usually consider spherical shells and not ellipsoidal ones, but some also introduce corrections to compensate that. Some standards or ITU recommendation are unclear about that. The method referenced above, in OneAxisEllipsoid, takes Earth oblateness into acount.

[baubin]

I would have just assumed that the ionospheric models all took oblateness into account. Thanks for the tip Luc.

[AlekLee]

Just a word of caution: the approach of finding the point with minimum altitude I proposed will break in situations like this, when the RX and TX are on the opposite sites of the Earth:

More specifically, this line computing the delay from the lowest altitude point to the RX:

delay2 = iono_model.pathDelay(min_alt_state, rx_topo_frame, epoch, frequency, None)

will fall over with a somewhat cryptic, to me, error:

*** java.lang.java.lang.ArrayIndexOutOfBoundsException: java.lang.ArrayIndexOutOfBoundsException: Index -1 out of bounds for length 183

This condition is fairly easy to check because:

IPdb [40]: min_alt_state.getPosition().getNorm() < Constants.EGM96_EARTH_EQUATORIAL_RADIUS
Out  [40]: True

Interestingly, computing the ionospheric path delay from TX to RX directly (bypassing the lowest altitude point) still works in cases like this:

IPdb [41]: iono_model.pathDelay(tx_state, rx_topo_frame, epoch, frequency, None)
Out  [41]: 0.8559012369666256

Honestly, I don’t know if it should because it makes little sense at first glance. However, I guess this is something for you to consider during your pending update, @baubin .