Application Programming Interface

This library contains components for PyRayHF software.

PyRayHF.library.azimuth_between_points(lon1_deg, lat1_deg, lon2_deg, lat2_deg)

Compute forward azimuth from point 1 to point 2.

Parameters:

lon1_deg, lat1_degfloat or array-like

Longitude and latitude and of the starting point [degrees].

lon2_deg, lat2_degfloat or array-like

Longitude and Latitude of the destination point [degrees].

Returns:

azimuth_degfloat or ndarray

Azimuth measured clockwise from geographic North [degrees], in the range [0, 360).

PyRayHF.library.build_mup_function(mup_field: ndarray, x_grid: ndarray, z_grid: ndarray, *, geometry: str = 'cartesian', R_E: float = None, bounds_error: bool = False, fill_value: float = nan) → Callable[[ndarray, ndarray], ndarray]

Construct callable for evaluating μ’(x, z) for group delay integration.

Parameters:

mup_fieldndarray

Grid of μ’ (group refractive index) values.

x_gridndarray

Horizontal coordinate grid:

• For “cartesian”: horizontal distance [km].

• For “spherical”: surface arc distance [km].

z_gridndarray

Vertical coordinate grid [km].

geometry{‘cartesian’, ‘spherical’}, optional

Coordinate system type. Default is ‘cartesian’.

R_Efloat, optional

Earth radius [km], used only for spherical geometry.

bounds_errorbool, optional

Passed to RegularGridInterpolator. If True, raise error outside grid.

fill_valuefloat, optional

Fill value used for extrapolation. Default np.nan.

Returns:

mup_funccallable

Function mup_func(x, z) that evaluates μ’ at given coordinates.

Notes

• For Cartesian geometry, (z, x) order is used in the interpolator.

• For Spherical geometry, (r, φ) order is used, where

  r = R_E + z and φ = x / R_E.

PyRayHF.library.build_refractive_index_interpolator_cartesian(z_grid: ndarray, x_grid: ndarray, n_field: ndarray, *, fill_value_n: float = nan, fill_value_grad: float = 0.0, bounds_error: bool = False, edge_order: int = 2) → Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]]

Construct interpolators for refractive index n(z, x) and its gradients.

Parameters:

z_gridndarray, shape (nz,)

Altitude coordinates [km], strictly increasing.

x_gridndarray, shape (nx,)

Horizontal coordinates [km], strictly increasing.

n_fieldndarray, shape (nz, nx)

Refractive index values on (z, x) grid.

fill_value_nfloat

Fill value for n outside grid (default NaN).

fill_value_gradfloat

Fill value for gradients outside grid (default 0.0).

bounds_errorbool

If True, raise error outside grid. If False, use fill values.

edge_orderint

Accuracy order for finite differences (default 2).

Returns:

n_and_gradcallable

Function (x, z) → (n, dndx, dndz).

PyRayHF.library.build_refractive_index_interpolator_spherical(z_grid: ndarray, x_grid: ndarray, n_field: ndarray, *, fill_value_n: float = nan, fill_value_grad: float = 0.0, bounds_error: bool = False, R_E: float = None, edge_order: int = 2) → Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]]

Construct interpolators for refractive index μ(r, φ) and its gradients.

Parameters:

z_gridndarray, shape (nz,)

Altitude coordinates [km], strictly increasing.

x_gridndarray, shape (nx,)

Horizontal coordinates [km], strictly increasing.

n_fieldndarray, shape (nr, nφ)

Refractive index values on (r, φ) grid.

fill_value_nfloat

Fill value for μ outside grid (default NaN).

fill_value_gradfloat

Fill value for gradients outside grid (default 0.0).

bounds_errorbool

If True, raise error outside grid. If False, use fill values.

R_Eflt

Radius of the Earth (km).

edge_orderint

Accuracy order for finite differences (default 2).

Returns:

n_and_grad_rphicallable

Function (φ, r) → (μ, ∂μ/∂r, ∂μ/∂φ).

PyRayHF.library.calculate_gcd(lon0, lat0, lon1, lat1)

Calculate great circle distance.

Parameters:

lon0array-like

value or np.array of longitudes of first position in degrees.

lat0array-like

value or np.array of latitudes of first position in degrees.

lon1array-like

value or np.array of longitudes of second position in degrees.

lat1array-like

value or np.array of latitudes of second position in degrees.

Returns:

gcdarray-like

Great circle distance in degrees.

Notes

This function calculates great circle distance in degrees.

PyRayHF.library.calculate_magnetic_field(year, month, day, lat, lon, aalt)

Get magnetic field strength and angle from vertical.

Parameters:

year, month, dayint

Date at which to evaluate B field

latarray-like

Latitude in degrees

lonarray-like

Longitude in degrees

aaltarray-like

Altitude array in km

Returns:

magarray-like

Magnetic field strength at each altitude [nT]

psiarray-like

Magnetic field angle from vertical at each altitude [degrees]

Notes

This returns the total field strength and angle from vertical

PyRayHF.library.constants()

Define constants for virtual height calculation.

Returns:

cpfloat

Plasma frequency constant in Hz per sqrt(m^-3). f_p [Hz] = cp * sqrt(n_e [m^-3]).

g_pfloat

Electron gyrofrequency constant in Hz/T (f_ce = g_p * B).

R_Efloat

Earth radius [km].

c_km_sfloat

Speed of light [km/s].

Notes

This function provides constants used in virtual height calculations.

PyRayHF.library.den2freq(density)

Convert plasma density to plasma frequency.

Parameters:

densityfloat or array-like

Plasma density in m^-3.

Returns:

frequencyfloat or array-like

Plasma frequency in Hz.

PyRayHF.library.earth_radius_at_latitude(latitude)

Calculate the radius of the Earth (in kilometers) for a given latitude.

Parameters:

latitudefloat

Geo Latitude in degrees.

Returns:

radiusfloat

Radius of the Earth in km.

PyRayHF.library.eval_refractive_index_and_grad(x: ndarray, z: ndarray, n_interp: RegularGridInterpolator, dn_dx_interp: RegularGridInterpolator, dn_dz_interp: RegularGridInterpolator) → Tuple[ndarray, ndarray, ndarray]

Evaluate refractive index and its gradients at (x, z).

Parameters:

xarray_like

Horizontal positions [km]. Scalar or array.

zarray_like

Altitudes [km]. Scalar or array.

n_interpRegularGridInterpolator

Interpolator for n(z, x).

dn_dx_interpRegularGridInterpolator

Interpolator for ∂n/∂x(z, x).

dn_dz_interpRegularGridInterpolator

Interpolator for ∂n/∂z(z, x).

Returns:

nnp.ndarray

Refractive index values at (x, z).

dndxnp.ndarray

∂n/∂x [1/km] values at (x, z).

dndznp.ndarray

∂n/∂z [1/km] values at (x, z).

Notes

Inputs x and z are broadcast to a common shape. Points are stacked as (z, x) before passing to interpolators. Output arrays match the broadcasted shape of inputs.

PyRayHF.library.event_ground(s: float, y: ndarray, z_ground_km: float) → float

Stop when ray hits or goes below the ground.

PyRayHF.library.event_x_left(s: float, y: ndarray, x_min_km: float) → float

Stop when ray exits left boundary.

PyRayHF.library.event_x_right(s: float, y: ndarray, x_max_km: float) → float

Stop when ray exits right boundary.

PyRayHF.library.event_z_bottom(s: float, y: ndarray, z_min_km: float) → float

Stop when ray leaves the bottom of the domain.

PyRayHF.library.event_z_top(s: float, y: ndarray, z_max_km: float) → float

Stop when ray leaves the top of the domain.

PyRayHF.library.find_X(n_e, f)

Calculate the square of the plasma freq over the square of the ion freq.

Parameters:

n_earray-like

Electron density in m^-3.

farray-like

Ionosonde frequency in Hz.

Returns:

Xfloat or array-like

Ratio (f_N / f)^2.

PyRayHF.library.find_Y(f, b)

Calculate Y: the electron gyrofrequency to ionosonde frequency ratio.

Parameters:

farray-like

Ionosonde frequency in Hz.

barray-like

Magnetic field magnitude in Tesla.

Returns:

Yarray-like

Electron gyrofrequency / ionosonde frequency.

PyRayHF.library.find_mean_gradient_error(atlon, atlat, arlon, arlat, year, month, day, UT, F107)

Estimate the mean horiz-gradient-induced foF2 error along T and R GCD.

Parameters:

atlonndarray

Transmitter longitudes [deg].

atlatndarray

Transmitter latitudes [deg].

arlonndarray

Receiver longitudes [deg].

arlatndarray

Receiver latitudes [deg].

yearint

Year of the PyIRI specification.

monthint

Month of the PyIRI specification.

dayint

Day of the PyIRI specification.

UTfloat

Universal Time [hours].

F107float

F10.7 solar flux index [sfu].

Returns:

a_mean_errorndarray

Mean percent foF2 deviation along each T-R path relative to the midpoint value [%]. Positive values indicate higher foF2 away from the midpoint.

F2_middict

PyIRI F2-region output at the midpoint locations, including foF2 and associated parameters.

Notes

For each transmitter-receiver (T-R) pair, this function samples the ionosphere along the great-circle path between the two locations, computes the foF2 variation relative to the midpoint value, and returns the mean percent deviation. The midpoint foF2 is treated as the reference, consistent with the midpoint-assumption framework.

Sampling is performed uniformly along the great-circle distance between each transmitter and receiver. foF2 (rather than NmF2) is used to quantify horizontal gradients, avoiding logarithmic scaling effects. The result provides a compact scalar metric describing how strongly horizontal structure violates the midpoint assumption for each link.

PyRayHF.library.find_mu_mup(X, Y, bpsi, mode, *, y_tol: float = 1e-12) → tuple[ndarray, ndarray]

Calculate the phase and group refractive indices (μ and μ’).

Parameters:

Xarray-like

Ratio of plasma and transmission frequencies squared.

Yarray-like

Ratio of electron gyrofrequency and transmission frequency.

bpsiarray-like

Angle ψ between wave vector and magnetic field in degrees.

modestr

‘O’ for ordinary or ‘X’ for extraordinary wave mode.

y_tolflt

Tollerance to determine if plasma is magnetized or not.

Returns:

muarray-like

Phase refractive index μ.

muparray-like

Group refractive index μ’.

Notes

When Y ≈ 0 (unmagnetized plasma), use isotropic formulas: mu = sqrt(1 - X) for X < 1, else NaN mup = 1 / mu Otherwise, fall through to the existing magnetized Appleton-Hartree logic.

PyRayHF.library.find_turning_point(z: ndarray, mu: ndarray, p: float) → float

Locate altitude where μ crosses the Snell invariant p.

Uses linear interpolation between bracketing nodes.

Parameters:

znp.ndarray

1D array of altitudes [km], monotonically increasing.

munp.ndarray

1D array of refractive index values corresponding to each altitude.

pfloat

Snell invariant (μ sinθ), constant along the ray path.

Returns:

z_turnfloat

Altitude of turning point [km].

PyRayHF.library.find_vh(X, Y, bpsi, dh, alt_min, mode)

Calculate virtual height for given mode.

Parameters:

Xarray-like

Plasma to ionosonde frequency ratio squared.

Yarray-like

Electron gyrofrequency to ionosonde frequency ratio.

bpsiarray-like

Angle between wave vector and magnetic field (degrees).

dharray-like

Vertical layer thickness in km.

alt_minfloat

Minimum altitude in km.

modestr

“O” or “X” mode.

Returns:

vharray-like

Virtual height in km.

PyRayHF.library.freq2den(frequency)

Convert plasma frequency to plasma density.

Parameters:

frequencyfloat or array-like

Plasma frequency in Hz.

Returns:

densityfloat or array-like

Plasma density in m^-3.

PyRayHF.library.generate_input_1D(year, month, day, UT, tlat, tlon, aalt, F107, save_path='')

Compute 1D PyIRI/IGRF input data for raytracing in PyRayHF.

Parameters:

yearint

Four digit year in C.E.

monthint

Integer month

dayint

Integer day of month

UTint

Universal time (UT) in hours

tlatfloat

Latitude of transmitter in degrees

tlonfloat

Longitude of transmitter in degrees

aaltarray-like

Array of altitude grid points in km (max 1000 km)

F107float

User provided F10.7 solar flux index in SFU.

save_pathstr

Full file path for saving output data (must include .p file extension) If left blank, the output will not be saved toa file. (default=’’)

Returns:

out_datadict

‘alt’ : z coordinate array of vertical grid in km ‘den’ : Array of electron density in m^-3 ‘bmag’ : Array of magnetic field strenth in T ‘bpsi’ : Array of magnetic field angle to vertical in degrees ‘F2’ : PyIRI output dictionary for F2 region at (tlon, xlat) ‘F1’ : PyIRI output dictionary for F1 region at (tlon, tlat) ‘E’ : PyIRI output dictionary for E region at (tlon, tlat) ‘Es’ : PyIRI output dictionary for Es region at (tlon, tlat) ‘year’ : Year used to run PyIRI for this data ‘month’ : Month used to run PyIRI for this data ‘day’ : Day of month used to run PyIRI for this data ‘UT’ : Unniversal time used to run PyIRI for this data ‘F107’ : F10.7 in SFU used to run PyIRI for this data ‘tlat’ : Latitude in degrees of transmitter ‘tlon’ : Longitude in degrees of transmitter

Notes

Given date/time, solar conditions, HF transmitter location, and desired grid parameters, this function uses PyIRI and IGRF to generate electron density and magnetic field data in a 1D (vertical) grid that can be used for vertical or 2D Snell’s-law ray tracing input to PyRayHF.

PyRayHF.library.generate_input_2D(year, month, day, UT, tlat, tlon, dx, aalt, gcd, az, F107, save_path='')

Compute 2D PyIRI/IGRF input data for raytracing in PyRayHF.

Parameters:

yearint

Four digit year in C.E.

monthint

Integer month

dayint

Integer day of month

UTint

Universal time (UT) in hours

tlatfloat

Latitude of transmitter in degrees

tlonfloat

Longitude of transmitter in degrees

dxfloat

Horizontal resolution of output grid in km

aaltarray-like

Array of altitude grid points in km (max 1000 km)

gcdfloat

Great circle distance in km from transmitter to end of domain (The horizontal span of the output domain)

azfloat

Azimuth from transmitter in degrees that defines the orientation of the output domain (measured clockwise from geographic North).

F107float

User provided F10.7 solar flux index in SFU.

save_pathstr

Full file path for saving output data (must include .p file extension) If left blank, the output will not be saved toa file. (default=’’)

Returns:

out_datadict

‘xgrid’ : x coordinate array of horizontal grid in km ‘zgrid’ : z coordinate array of vertical grid in km ‘xlat’ : Latitude in degrees of each point in xgrid ‘xlon’ : Longitude in degrees of each point in xgrid ‘den’ : 2D array of electron density in m^-3 ‘bmag’ : 2D array of magnetic field strenth in T ‘bpsi’ : 2D array of magnetic field angle to vertical in degrees ‘F2’ : PyIRI output dictionary for F2 region at all (xlon, xlat) ‘F1’ : PyIRI output dictionary for F1 region at all (xlon, xlat) ‘E’ : PyIRI output dictionary for E region at all (xlon, xlat) ‘Es’ : PyIRI output dictionary for Es region at all (xlon, xlat) ‘year’ : Year used to run PyIRI for this data ‘month’ : Month used to run PyIRI for this data ‘day’ : Day of month used to run PyIRI for this data ‘UT’ : Unniversal time used to run PyIRI for this data ‘F107’ : F10.7 in SFU used to run PyIRI for this data ‘tlat’ : Latitude in degrees of transmitter ‘tlon’ : Longitude in degrees of transmitter ‘az’ : Azimuth from transmitter in degrees (measured clockwise from

geographic north) that defines the orientation of the output domain

Notes

Given date/time, solar conditions, HF transmitter location, ray azimuth, and desired grid parameters, this function uses PyIRI and IGRF to generate electron density and magnetic field data in a 2D (horizontal VS vertical) grid that can be used for oblique ray tracing input to PyRayHF.

PyRayHF.library.great_circle_point(tlat, tlon, gcd, az)

Get lat/lon of a GCD point from an origin point.

Parameters:

tlat, tlonfloat

Lattitude and longitude of origin [degrees]

gcdarray-like

Great circle distance from origin to destination point [km]

azfloat

Azimuth (clockwise from north) from origin to destination point [degrees]

Returns:

rlat, rlonarray-like

Lattitude and longitude of destination point [degrees]

Notes

Assumes spherical earth (not ellipsoid)

PyRayHF.library.make_n_and_grad(n_interp: RegularGridInterpolator, dn_dx_interp: RegularGridInterpolator, dn_dz_interp: RegularGridInterpolator) → Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]]

Return a function (x, z) -> (n, dndx, dndz).

PyRayHF.library.minimize_parameters(F2, F1, E, f_in0, vh_obs0, alt, b_mag, b_psi, method='brute', percent_sigma=20.0, step=1.0, mode='O', n_points=200, bottom_type='B_bot')

Minimize F2 layer parameters (hmF2 and B_bot) to fit observed VH.

Parameters:

F2dict

Initial F2 layer parameters. Must include ‘Nm’, ‘hm’, and ‘B_bot’.

F1dict

Initial F1 layer parameters.

Edict

Initial E layer parameters.

f_in0ndarray

Input frequencies [MHz].

vh_obs0ndarray

Observed virtual heights [km].

altndarray

Altitude array [km].

b_magndarray

Magnetic field magnitude array [nT].

b_psindarray

Magnetic field dip angle array [degrees].

methodstr

Method of minimization in lmfit: “brute”: (default) A grid search method for finding a global minimum. “levenberg-marquardt”: Generally fast and effective for many curve-fitting needs. “powell”: Another derivative-free method.

percent_sigmaflt

How far off from the background value to deviate. Default is 20%. If the speed needs to be increase, decrease this parameter.

stepflt

Step size in km for brute minimization. If the speed needs to be increase, increase this parameter.

modestr

‘O’ or ‘X’ mode. The default is ‘O’.

n_pointsint

Number of desired vertical grid points. Default is 200.

bottom_typestr

Type of the F2 bottom formalizm. Default is ‘B_bot’, which construncts bottom side of F2 using a single thickness parameter. Other option is ‘B0_B1’ which uses an IRI formalizm.

Returns:

vh_resultndarray

Virtual height after parameter fitting [km].

EDP_resultndarray

Reconstructed electron density profile after fitting [m^-3].

F2_fitdict

F2 dictionary with updated parameters.

PyRayHF.library.model_VH(F2, F1, E, f_in, alt, b_mag, b_psi, mode='O', n_points=200, bottom_type='B_bot')

Compute vertical virtual height using a modeled EDP and raytrace.

Parameters:

F2dict

Dictionary of F2 layer parameters. Must include: - ‘Nm’: peak electron density (NmF2) - ‘hm’: peak height (hmF2) - ‘B_bot’: thickness of the bottomside of the F2 layer

F1dict

Dictionary of F1 layer parameters. Must include: - ‘P’: shape factor or profile parameter

Edict

Dictionary of E layer parameters. Must include: - ‘hm’: peak height of the E layer

f_inndarray

Input frequency [MHz].

altndarray

1D array of altitudes [km].

b_magndarray

1D array of magnetic field magnitudes [nT].

b_psindarray

1D array of magnetic field dip angles [rad].

modestr

‘O’ or ‘X’ mode. Default is ‘O’ mode.

n_pointsint

Number of vertical grid points. Default is 200.

bottom_typestr

Type of the F2 bottom formalizm. Default is ‘B_bot’, which construncts bottom side of F2 using a single thickness parameter. Other option is ‘B0_B1’ which uses an IRI formalizm.

Returns:

vh_Ondarray

Virtual height trace (O-mode) [km].

EDPndarray

Reconstructed electron density profile [m^-3].

PyRayHF.library.n_and_grad(x: ndarray, z: ndarray, n_interp: RegularGridInterpolator, dn_dx_interp: RegularGridInterpolator, dn_dz_interp: RegularGridInterpolator) → Tuple[ndarray, ndarray, ndarray]

Evaluate n, ∂n/∂x, and ∂n/∂z at points (x, z).

Parameters:

xarray_like

Horizontal positions [km]; scalar or array.

zarray_like

Altitudes [km]; scalar or array.

Returns:

nnp.ndarray

Refractive index at (x, z). Shape equals the broadcasted shape of inputs.

dndxnp.ndarray

Partial derivative ∂n/∂x [1/km] at (x, z). Same shape as n.

dndznp.ndarray

Partial derivative ∂n/∂z [1/km] at (x, z). Same shape as n.

Notes

Inputs are broadcast to a common shape, flattened for one batched interpolation, then reshaped back—this minimizes overhead during ODE stepping. Outside the grid hull, values are fill_value_n and fill_value_grad (unless bounds_error=True).

PyRayHF.library.n_and_grad_rphi(phi: ndarray, r: ndarray, n_interp: RegularGridInterpolator, dn_dr_interp: RegularGridInterpolator, dn_dphi_interp: RegularGridInterpolator) → Tuple[ndarray, ndarray, ndarray]

Evaluate μ(r, φ) and its gradients at given coordinates.

Parameters:

phiarray_like

Azimuthal coordinate [rad]. Scalar or array. Will broadcast with r.

rarray_like

Radial coordinate [km] (Earth radius + altitude). Scalar or array.

n_interpRegularGridInterpolator

Interpolator for μ(r, φ).

dn_dr_interpRegularGridInterpolator

Interpolator for ∂μ/∂r.

dn_dphi_interpRegularGridInterpolator

Interpolator for ∂μ/∂φ.

Returns:

nnp.ndarray

Refractive index μ at (r, φ). Shape matches broadcasted inputs.

dn_drnp.ndarray

Partial derivative ∂μ/∂r at (r, φ).

dn_dphinp.ndarray

Partial derivative ∂μ/∂φ at (r, φ).

PyRayHF.library.oblique_to_vertical(range_km, group_path_km, freq_oblique_mhz, R_E=6371.0)

Convert oblique ionogram to vertical equivalent using spherical Earth.

Parameters:

range_kmfloat

Ground distance between transmitter and receiver [km]

group_path_kmarray-like

Oblique group path (total propagation distance) [km]

freq_oblique_mhzarray-like

Oblique frequency [MHz]

R_Eflt

Earth’s radius in [km]. Default is 6371 km.

Returns:

freq_vertical_mhzarray-like

Equivalent vertical frequency [MHz]

height_virtual_kmarray-like

Virtual height at midpoint [km]

PyRayHF.library.ray_rhs_cartesian(s: float, y: ndarray, n_and_grad: Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]], renormalize_every: int, eval_counter: Dict[str, int]) → ndarray

Right-hand side of the 2D ray equations in Cartesian coordinates.

Parameters:

sfloat

Arc length [km].

yndarray

State vector [x, z, vx, vz].

n_and_gradcallable

Function returning (n, dndx, dndz) at (x, z).

renormalize_everyint

Frequency (in calls) to re-normalize direction vector.

eval_counterdict

Mutable counter used for tracking RHS evaluations.

Returns:

dydsndarray

Derivatives [dx/ds, dz/ds, dvx/ds, dvz/ds].

PyRayHF.library.regrid_to_nonuniform_grid(f, n_e, b, bpsi, aalt, mode='O', n_points=200, dh=1e-06)

Regrid profile to smooth non-uniform vertical grid.

Parameters:

farray-like

Ionosonde frequency in Hz.

n_earray-like

Electron density in m^-3.

barray-like

Magnetic field magnitude.

bpsiarray-like

Angle to magnetic field vector in degrees.

aaltarray-like

Altitude profile in km.

modestr

‘O’ or ‘X’ propagation mode. Default ‘O’.

n_pointsint

Points in new vertical grid.

dhflt

How close to the reflection height do we want to get. Default is 1e-6 km.

Returns:

regriddeddict

Dictionary with re-gridded arrays

PyRayHF.library.residual_VH(params, F2_init, F1_init, E_init, f_in, vh_obs, alt, b_mag, b_psi, mode='O', n_points=200, bottom_type='B_bot')

Compute the residuals between observed and modeled VHs.

Parameters:

paramslmfit.Parameters

Parameters to be optimized, containing: - ‘NmF2’: peak electron density of F2 layer - ‘hmF2’: peak height of F2 layer - ‘B_bot’: thickness of F2 bottomside

F2_initdict

Initial F2 layer parameters.

F1_initdict

Initial F1 layer parameters.

E_initdict

Initial E layer parameters.

f_infloat

Input frequency [MHz].

vh_obsndarray

Observed virtual heights [km].

altndarray

Altitude array [km].

b_magndarray

Magnetic field magnitude array [nT].

b_psindarray

Magnetic field dip angle array [rad].

modestr

‘O’ or ‘X’ mode. Default is ‘O’ mode.

n_pointsint

Number of vertical grid points. Default is 200.

bottom_typestr

Type of the F2 bottom formalizm. Default is ‘B_bot’, which construncts bottom side of F2 using a single thickness parameter. Other option is ‘B0_B1’ which uses an IRI formalizm.

Returns:

residualndarray

Flattened array of residuals between observed and modeled virtual heights [km].

PyRayHF.library.rhs_spherical(s: float, y: ndarray, n_and_grad_rphi: Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]], renormalize_every: int, eval_counter: Dict[str, int]) → ndarray

Calculate the right-hand side of the spherical ray equations.

Parameters:

sfloat

Current arc length along the ray [km].

yndarray, shape (4,)

State vector [r, φ, v_r, v_φ]: - r : radial coordinate [km] - φ : azimuthal angle [rad] - v_r : radial direction cosine - v_φ : tangential direction cosine

n_and_grad_rphicallable

Function (φ, r) → (μ, ∂μ/∂r, ∂μ/∂φ), typically constructed using build_refractive_index_interpolator_rphi.

renormalize_everyint

Frequency (in RHS evaluations) to re-normalize the tangent vector to ensure |v| = 1. Use 0 or None to disable.

eval_counterdict

Dictionary used to track number of RHS evaluations: e.g., {‘n’: 0} will be incremented in-place.

Returns:

dydsndarray

Derivatives with respect to arc length s: [dr/ds, dφ/ds, dv_r/ds, dv_φ/ds].

Notes

• Implements 2D spherical geometry (flat-Earth limit not assumed).

• The equations conserve |v| ≈ 1 under small step sizes.

• NaN or non-positive μ values return zero derivatives (halts ray).

Computes derivatives for the ODE system governing 2D spherical ray propagation (r, φ) in a refractive index field μ(r, φ). The equations describe the evolution of position and tangent components along the raypath with respect to arc length s.

The system is defined as:

\[ \begin{align}\begin{aligned}\begin{split}\\frac{dr}{ds} = v_r \\qquad \\frac{dφ}{ds} = \\frac{v_φ}{r}\end{split}\\\begin{split}\\frac{dv_r}{ds} = \\frac{1}{μ} \\left[ \\frac{∂μ}{∂r} - (∇μ · v)v_r \\right] + \\frac{v_φ^2}{r}\end{split}\\\begin{split}\\frac{dv_φ}{ds} = \\frac{1}{μ} \\left[ \\frac{1}{r} \\frac{∂μ}{∂φ} - (∇μ · v)v_φ \\right] - \\frac{v_r v_φ}{r}\end{split}\end{aligned}\end{align} \]

where the gradient projection is given by: .. math:

∇μ · v = \\frac{∂μ}{∂r} v_r + \\frac{1}{r} \\frac{∂μ}{∂φ} v_φ

References

• Haselgrove, J., “The Hamiltonian Ray Path Equations”, Proc. IEE (1955)

• Budden, K. G., The Propagation of Radio Waves, Cambridge Univ. Press,

1985

PyRayHF.library.save_to_file(output, file_path)

Save dictionary to a pickle file.

Parameters:

outputdict

Dictionary to save to pickle file

file_pathstring

Full filepath (including .p extension) of output file

PyRayHF.library.smooth_nonuniform_grid(start, end, n_points, sharpness)

Generate smooth non-uniform grid from start to end.

Parameters:

startfloat

endfloat

n_pointsint

sharpnessfloat

Controls how sharply resolution increases near end.

Returns:

xndarray

Non-uniformly spaced grid.

PyRayHF.library.tan_from_mu_scalar(mu_val: float, p: float) → float

Compute tanθ safely for Snell’s law in plasma.

Parameters:

mu_valfloat

Phase refractive index μ at altitude.

pfloat

Snell invariant (μ0 sinθ0).

Returns:

tanθfloat

Tangent of propagation angle relative to vertical.

Notes

This function finds tah_theta from μ * sin(θ) = p. tan(θ) = sin(θ) / cos(θ)​, and cos(θ) = sqrt(1−sin^2(θ)) tan(θ) = ​(p/μ​) / sqrt(1−(p/μ)^2) = ​p​ / sqrt(μ^2 − p^2)

PyRayHF.library.trace_ray_cartesian_gradient(n_and_grad: Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]], mup_func: Callable[[ndarray, ndarray], ndarray], x0_km: float, z0_km: float, elevation_deg: float, s_max_km: float = 5000.0, *, rtol: float = 1e-07, atol: float = 1e-09, max_step_km: float | None = None, z_ground_km: float = 0.0, z_min_km: float = -1.0, z_max_km: float = 1000.0, x_min_km: float = -1000000.0, x_max_km: float = 1000000.0, renormalize_every: int = 50) → Dict[str, Any]

Trace a 2D ray in a horizontally varying refractive index field μ(x, z).

Geometry is integrated via:

dr/ds = v, ||v|| = 1 dv/ds = (1/μ) [∇μ - (∇μ·v) v]

Group delay is integrated using μ’(x, z):

τ = ∫ (μ’/c) ds

Parameters:

n_and_gradcallable

(x, z) → (μ, ∂μ/∂x, ∂μ/∂z). Usually built with build_refractive_index_interpolator(…).

mup_funccallable

(x, z) → μ’(x, z). Must be built with build_mup_function(…).

x0_km, z0_kmfloat

Launch point [km].

elevation_degfloat

Launch elevation above horizontal [deg].

s_max_kmfloat

Max path length [km].

rtol, atolfloat

ODE solver tolerances.

max_step_kmfloat or None

Max solver step size [km].

z_ground_km, z_min_km, z_max_kmfloat

Vertical domain limits [km].

x_min_km, x_max_kmfloat

Horizontal domain limits [km].

renormalize_everyint

Re-normalize v every N evaluations.

Returns:

resultdict

Notes

• This model assumes a 2D Cartesian (flat-Earth) geometry.

• μ controls bending; μ’ controls group delay.

• NaNs or invalid μ terminate integration.

Return dictionary keys: ‘sol’ : OdeResult, # full solve_ivp solution object ‘t’ : ndarray, # integration parameter (path length) [km] ‘x’ : ndarray, # horizontal coordinate [km] ‘z’ : ndarray, # altitude [km] ‘vx’ : ndarray, # horizontal unit-velocity component ‘vz’ : ndarray, # vertical unit-velocity component ‘status’ : str, # ‘ground’, ‘domain’, ‘length’, etc. ‘group_path_km’ : float, # geometric path length [km] ‘group_delay_sec’ : float, # group delay τ = ∫ μ’/c ds [s] ‘x_midpoint’ : float, # midpoint x-coordinate [km] ‘z_midpoint’ : float, # midpoint altitude [km] ‘ground_range_km’ : float, # landing x-coordinate or NaN [km] ‘x_apex_km’ : float, # apex horizontal coordinate [km] ‘z_apex_km’ : float # apex altitude [km]

PyRayHF.library.trace_ray_cartesian_snells(f0_Hz: float, elevation_deg: float, alt_km: ndarray, Ne: ndarray, Babs: ndarray, bpsi: ndarray, mode: str) → Dict[str, float]

Stratified Snell’s law ray tracing (flat Earth, 2D Cartesian).

Parameters:

f0_Hzfloat

Frequency [Hz].

elevation_degfloat

Launch elevation above horizontal [deg].

alt_kmndarray

Altitude grid [km].

Nendarray

Electron density [el/m^3].

Babsndarray

Magnetic field strength [T].

bpsindarray

Angle between B and k-vector [degrees].

modestr

Wave mode: “O” or “X”.

Returns:

resultdict

Notes

This function models how a high-frequency radio wave propagates through the ionosphere, using Snell’s law adapted for a plasma medium. It calculates the trajectory of the ray as it leaves the ground, bends through the ionized atmosphere, reaches a turning point, and (if conditions allow) returns back toward Earth.

Background: Snell’s Law in a Plasma. In a uniform dielectric, Snell’s law states that nsinθ=constant, where n is the refractive index and θ is the propagation angle relative to the vertical. In a plasma, the refractive index isn’t constant but depends on: electron density (affects plasma frequency), magnetic field (splits O and X modes), wave frequency, and angle between wave vector and magnetic field. This gives two possible wave modes: the ordinary (O) and extraordinary (X) mode, each with a different effective refractive index. The function uses auxiliary functions find_X, find_Y, and find_mu_mup to compute these refractive indices as functions of altitude. Thus, the plasma-modified Snell’s law is applied: μ’sinθ=constant, where μ’ is the transverse refractive index for the chosen wave mode.

Specifics: Geometry (bending) uses phase index μ. Group delay integrates group index μ’ (mup). Down-leg is a perfect mirror of the up-leg about the apex.

Return dictionary has the following keys:

‘x’ : ndarray # surface distance along the ray [km] ‘z’ : ndarray # altitude along the ray [km] ‘group_path_km’ : float # total geometric path length [km] ‘group_delay_sec’: float # group delay [s] ‘x_midpoint’ : float # midpoint horizontal coordinate [km] ‘z_midpoint’ : float # midpoint altitude [km] ‘ground_range_km’: float # surface distance to landing point [km] ‘x_apex_km’ : float # apex horizontal coordinate [km] ‘z_apex_km’ : float # apex altitude [km]

PyRayHF.library.trace_ray_spherical_gradient(n_and_grad_rphi: Callable[[ndarray, ndarray], Tuple[ndarray, ndarray, ndarray]], mup_func: Callable[[ndarray, ndarray], ndarray], x0_km: float, z0_km: float, elevation_deg: float, s_max_km: float = 6000.0, *, R_E: float | None = None, z_ground_km: float = 0.0, r_max_km: float | None = None, phi_min: float = -3.141592653589793, phi_max: float = 3.141592653589793, rtol: float = 1e-07, atol: float = 1e-09, max_step_km: float | None = 2.0, renormalize_every: int = 50) → Dict[str, Any]

Raytrace in 2-D spherical using μ for geometry and μ’ for delay.

Parameters:

n_and_grad_rphicallable

(φ, r) → (μ, ∂μ/∂r, ∂μ/∂φ).

mup_funccallable (REQUIRED)

(x, z) → μ’(x, z). Build with build_mup_function(…, geometry=”spherical”).

x0_km, z0_kmfloat

Launch point [km]; x0 is surface arc distance, z0 altitude.

elevation_degfloat

Launch elevation above local horizontal [deg].

s_max_kmfloat

Max arclength to integrate [km].

R_Efloat, optional

Earth radius [km]. Defaults to constants()[2] if None.

z_ground_kmfloat

Ground altitude [km]; stop when r ≤ R_E + z_ground_km.

r_max_kmfloat, optional

Max r; default R_E + 1200 km if None.

phi_min, phi_maxfloat

Azimuth bounds [rad].

rtol, atolfloat

ODE tolerances.

max_step_kmfloat or None

Max solver step [km].

renormalize_everyint

Re-normalize v every N RHS calls.

Returns:

dict

Notes

Return dictionary has the following keys: ‘t’ : ndarray, # integration parameter (path length) [km] ‘r’ : ndarray, # radial coordinate [km] ‘phi’ : ndarray, # angular coordinate [rad] ‘v_r’ : ndarray, # radial component of unit tangent vector ‘v_phi’ : ndarray, # angular component of unit tangent vector ‘x’ : ndarray, # surface distance [km] ‘z’ : ndarray, # altitude above Earth’s surface [km] ‘status’ : str, # termination condition (ground, domain, ..) ‘group_path_km’ : float, # geometric path length [km] ‘group_delay_sec’ : float, # group delay [s] ‘x_midpoint’ : float, # midpoint horizontal coordinate [km] ‘z_midpoint’ : float, # midpoint altitude [km] ‘ground_range_km’ : float, # landing surface distance [km] ‘x_apex_km’ : float, # apex horizontal coordinate [km] ‘z_apex_km’ : float # apex altitude [km]

Geometry: Integrates in spherical coordinates (r, φ) with tangent v = (v_r, v_φ):

dr/ds = v_r dφ/ds = v_φ / r dv_r/ds = (1/μ)[∂μ/∂r - (∇μ·v)v_r] + (v_φ²)/r dv_φ/ds = (1/μ)[(∂μ/∂φ)/r - (∇μ·v)v_φ] - (v_r v_φ)/r

where ∇μ·v = (∂μ/∂r)v_r + (∂μ/∂φ)(v_φ / r).

Delay Group delay integrates μ’ along the path:

τ = ∫ ( μ’(x, z) / c ) ds

with x = R_E φ (surface arc) and z = r − R_E.

PyRayHF.library.trace_ray_spherical_snells(f0_Hz: float, elevation_deg: float, alt_km: ndarray, Ne: ndarray, Babs: ndarray, bpsi: ndarray, mode: str = 'O', *, dz_target_km: float = 1.0, apex_boost: float = 200.0, max_substeps: int = 400, R_E: float | None = None) → Dict[str, float]

Stratified Snell’s law ray tracing (spherical Earth, 2D geometry).

Parameters:

f0_Hzfloat

Radio frequency [Hz].

elevation_degfloat

Launch elevation above local horizontal [degrees].

alt_kmndarray

Altitude grid [km].

Nendarray

Electron density profile [el/m³].

Babsndarray

Magnetic field strength [T].

bpsindarray

Magnetic field inclination [degrees].

modestr, default ‘O’

Wave mode: ‘O’ (ordinary) or ‘X’ (extraordinary).

dz_target_kmfloat, default 1.0

Target altitude increment for apex refinement [km].

apex_boostfloat, default 200.0

Sharpness multiplier that increases substeps near turning points.

max_substepsint, default 400

Hard limit on number of adaptive substeps per altitude interval.

R_Efloat or None, optional

Earth radius [km]. If None, defaults to value from constants().

Returns:

resultdict

Notes

Physical formulation:

In spherical geometry, Snell’s invariant becomes:

μ * r * sin(θ) = constant

where μ is the phase refractive index, r = R_E + z, and θ is the angle between the ray and the local vertical.

The ray is launched from ground level, bent according to μ(r), reflected at the turning point where μ * r = p, and mirrored down.

Group delay:

The propagation delay is integrated along the ray using the group refractive index μ’ (mup):

τ = ∫ (μ’ / c) ds

where c is the speed of light in vacuum.

Apex refinement:

Near the turning point, the derivative

dφ/dz = p / [ r * sqrt((μ r)² - p²) ]

becomes sharply peaked, making coarse integration unstable. To resolve this, the algorithm adaptively subdivides altitude steps where |μ r - p| is small, increasing resolution toward the apex while capping the total substeps with max_substeps.

Comparison to Cartesian model:

This function extends the flat-Earth Snell’s law formulation (trace_ray_cartesian_snells) to spherical geometry. For large R_E, both solutions converge within numerical precision.

Return dictionary has the following keys: ‘x’ : ndarray, # horizontal coordinate [km] ‘z’ : ndarray, # altitude along the ray [km] ‘group_path_km’ : float, # total geometric path length [km] ‘group_delay_sec’ : float, # group delay [s] ‘x_midpoint’ : float, # midpoint horizontal coordinate [km] ‘z_midpoint’ : float, # midpoint altitude [km] ‘ground_range_km’ : float, # landing surface distance [km] ‘x_apex_km’ : float, # apex horizontal coordinate [km] ‘z_apex_km’ : float # apex altitude [km]

PyRayHF.library.vertical_forward_operator(freq, den, bmag, bpsi, alt, mode='O', n_points=200)

Calculate virtual height from ionosonde freq and ion profile.

Parameters:

freqarray-like

Frequency in MHz.

denarray-like

Electron density in m^-3.

bmagarray-like

Magnetic field magnitude in Tesla.

bpsiarray-like

Angle to magnetic field vector.

altarray-like

Altitude profile in km.

modestr

‘O’ or ‘X’ propagation mode. Default ‘O’.

n_pointsint

Number of vertical grid points. Default is 200.

Returns:

vhndarray

Virtual height in km.

PyRayHF.library.vertical_to_magnetic_angle(inclination_deg)

Calculate angle between vertical and magnetic field vector.

Parameters:

inclination_degfloat or ndarray

Magnetic inclination in degrees (positive = downward).

Returns:

vertical_anglefloat or ndarray

Angle between vertical and magnetic field in degrees.