Skip to content

hfpytrace.pharlap

Package

MATLAB Engine wrapper for PHaRLAP 2D and 3D calls.

Key Classes

Class Engine

Key Methods

Function get_matlab_pharlap_lib()
Method Engine.run_pharlap()
Method Engine.run_pharlap_3d()
Method Engine.run_pharlap_3d_sp()
Method Engine.close()

API

hfpytrace.pharlap

PHaRLAP MATLAB engine wrappers.

Provides Python wrappers around the PHaRLAP (Provision of High-frequency Raytracing Laboratory for Propagation studies) MATLAB library, enabling 2D and 3D HF ray-tracing from Python via the MATLAB Engine API.

Two 3D formulations are supported:

  • raytrace_3d — WGS84 ellipsoidal Earth model.
  • raytrace_3d_sp — spherical Earth model (used for the SCurve study).

Requires

matlab.engine : MATLAB Engine API for Python (ships with MATLAB installation). PHaRLAP : MATLAB HF ray-tracing library installed at PHARLAP_LIB_PATH.

Classes / Functions

get_matlab_pharlap_lib Resolve the PHaRLAP library path for a given version tag. Engine MATLAB engine lifecycle manager with PHaRLAP ray-tracing methods.

Engine

MATLAB engine lifecycle manager with PHaRLAP ray-tracing methods.

Starts a MATLAB engine session, adds the PHaRLAP library to the MATLAB path, and exposes Python-callable wrappers for:

  • :meth:run_pharlap — 2D spherical ray-tracing (raytrace_2d_sp).
  • :meth:run_pharlap_3d — 3D WGS84 ray-tracing (raytrace_3d).
  • :meth:run_pharlap_3d_sp — 3D spherical ray-tracing (raytrace_3d_sp).
Parameters
str or None

Explicit path to the PHaRLAP MATLAB library directory. None resolves the path via :func:get_matlab_pharlap_lib.

Notes

Call :meth:close (or use as a context manager) to cleanly quit the MATLAB session and release memory.

Source code in hfpytrace/pharlap.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
class Engine:
    """MATLAB engine lifecycle manager with PHaRLAP ray-tracing methods.

    Starts a MATLAB engine session, adds the PHaRLAP library to the MATLAB
    path, and exposes Python-callable wrappers for:

    * :meth:`run_pharlap` — 2D spherical ray-tracing (``raytrace_2d_sp``).
    * :meth:`run_pharlap_3d` — 3D WGS84 ray-tracing (``raytrace_3d``).
    * :meth:`run_pharlap_3d_sp` — 3D spherical ray-tracing (``raytrace_3d_sp``).

    Parameters
    ----------
    lib_path : str or None
        Explicit path to the PHaRLAP MATLAB library directory.  ``None``
        resolves the path via :func:`get_matlab_pharlap_lib`.

    Notes
    -----
    Call :meth:`close` (or use as a context manager) to cleanly quit the
    MATLAB session and release memory.
    """

    def __init__(self, lib_path: str = None):
        self.eng = matlab.engine.start_matlab()
        env_path = get_matlab_pharlap_lib() if lib_path is None else lib_path
        self.eng.addpath(self.eng.genpath(env_path), nargout=0)
        logger.info("Matlab engine started and library path added.")
        return

    def close(self):
        logger.info("Closing Matlab engine.")
        self.eng.quit()
        return

    def run_pharlap(
        self,
        ne_grid: np.ndarray,
        collision_freq: np.ndarray,
        elevs: np.ndarray,
        rb: float,
        freqs: np.ndarray,
        irreg: np.ndarray,
        nhops: int = 1,
        tol: float = 1e-7,
        radius_earth: float = 6371,
        irregs_flag: int = 0,
        start_height: int = 50,
        height_inc: int = 1,
        range_inc: int = 1,
    ):
        """Run PHaRLAP 2D spherical ray-tracing (``raytrace_2d_sp``).

        Parameters
        ----------
        ne_grid : np.ndarray
            2D electron density grid, shape (nalt, nrange), in cm⁻³.
        collision_freq : np.ndarray
            2D collision frequency grid, same shape as ``ne_grid``, in Hz.
        elevs : np.ndarray, 1D
            Launch elevation angles [°].
        rb : float
            Ground range of the ionospheric profile endpoint [km].
        freqs : np.ndarray, 1D
            Operating frequencies [MHz].
        irreg : np.ndarray
            Irregularity array (passed to PHaRLAP; shape depends on version).
        nhops : int
            Number of ionospheric hops.  Default ``1``.
        tol : float
            ODE solver tolerance.  Default ``1e-7``.
        radius_earth : float
            Earth radius [km].  Default ``6371``.
        irregs_flag : int
            Irregularity flag (0 = no irregularities).  Default ``0``.
        start_height : int
            Starting height for the ray [km].  Default ``50``.
        height_inc, range_inc : int
            Grid increment in height and range [km].  Defaults ``1``.

        Returns
        -------
        ray_data : SimpleNamespace
            Per-ray scalar parameters (group_range, phase_path, …).
        ray_path_data : SimpleNamespace
            Per-ray path arrays (group_range, height, …).
        """
        logger.info("Running Pharlap...")
        self.eng.eval("close all; clear all; clc;", nargout=0)

        self.eng.workspace["ne_grid"] = matlab.double(ne_grid.tolist())
        self.eng.workspace["elevs"] = matlab.double(elevs.tolist())
        self.eng.workspace["rb"] = rb
        self.eng.workspace["collision_freq"] = matlab.double(collision_freq.tolist())
        self.eng.workspace["freqs"] = matlab.double(freqs.tolist())
        self.eng.workspace["nhops"] = nhops
        self.eng.workspace["tol"] = tol
        self.eng.workspace["radius_earth"] = radius_earth
        self.eng.workspace["irregs_flag"] = irregs_flag
        self.eng.workspace["start_height"] = start_height
        self.eng.workspace["height_inc"] = height_inc
        self.eng.workspace["range_inc"] = range_inc
        self.eng.workspace["irreg"] = matlab.double(irreg.tolist())

        self.eng.eval(
            """
            [ray_data, ray_path_data] = ...
                raytrace_2d_sp(elevs, rb, freqs, nhops, tol, ...
                radius_earth, irregs_flag, ne_grid, ne_grid, ...
                collision_freq, start_height, height_inc, range_inc, irreg ...
                );
            """,
            nargout=0,
        )
        try:
            ray_data, ray_path_data = (
                self.eng.workspace["ray_data"],
                self.eng.workspace["ray_path_data"],
            )
        except ValueError:
            # MATLAB Engine cannot directly return non-scalar struct arrays.
            self.eng.eval(
                "ray_data_json = jsonencode(ray_data); "
                "ray_path_data_json = jsonencode(ray_path_data);",
                nargout=0,
            )
            ray_data = json.loads(self.eng.workspace["ray_data_json"])
            ray_path_data = json.loads(self.eng.workspace["ray_path_data_json"])
        logger.info("Pharlap run completed.")
        return to_namespace(ray_data), to_namespace(ray_path_data)

    @staticmethod
    def _as_matlab_double(arr, ensure_row: bool = False):
        a = np.asarray(arr, dtype=float)
        if a.ndim == 0:
            return float(a)
        if ensure_row and a.ndim == 1:
            a = a.reshape(1, -1)
        return matlab.double(a.tolist())

    def _fetch_struct_array(self, keys: list[str]):
        try:
            out = [self.eng.workspace[k] for k in keys]
        except ValueError:
            # MATLAB Engine cannot directly return non-scalar struct arrays.
            enc = "; ".join([f"{k}_json = jsonencode({k})" for k in keys]) + ";"
            self.eng.eval(enc, nargout=0)
            out = [json.loads(self.eng.workspace[f"{k}_json"]) for k in keys]
        return [to_namespace(v) for v in out]

    def run_pharlap_3d(
        self,
        origin_lat: float,
        origin_lon: float,
        origin_ht: float,
        elevs: np.ndarray,
        ray_bearings: np.ndarray,
        freqs: np.ndarray,
        iono_en_grid: np.ndarray,
        iono_en_grid_5: np.ndarray,
        collision_freq: np.ndarray,
        iono_grid_parms: np.ndarray,
        Bx: np.ndarray,
        By: np.ndarray,
        Bz: np.ndarray,
        geomag_grid_parms: np.ndarray,
        OX_mode: int = 1,
        nhops: int = 1,
        tol: float | np.ndarray = 1e-7,
        ray_state_vec_in=None,
    ):
        """
        Wrapper for PHaRLAP `raytrace_3d` (WGS84 formulation).
        """
        logger.info("Running PHaRLAP raytrace_3d...")
        self.eng.eval("close all; clear all; clc;", nargout=0)

        self.eng.workspace["origin_lat"] = float(origin_lat)
        self.eng.workspace["origin_lon"] = float(origin_lon)
        self.eng.workspace["origin_ht"] = float(origin_ht)
        self.eng.workspace["elevs"] = self._as_matlab_double(elevs, ensure_row=True)
        self.eng.workspace["ray_bearings"] = self._as_matlab_double(
            ray_bearings, ensure_row=True
        )
        self.eng.workspace["freqs"] = self._as_matlab_double(freqs, ensure_row=True)
        self.eng.workspace["OX_mode"] = int(OX_mode)
        self.eng.workspace["nhops"] = int(nhops)
        self.eng.workspace["tol"] = self._as_matlab_double(tol, ensure_row=True)

        self.eng.workspace["iono_en_grid"] = self._as_matlab_double(iono_en_grid)
        self.eng.workspace["iono_en_grid_5"] = self._as_matlab_double(iono_en_grid_5)
        self.eng.workspace["collision_freq"] = self._as_matlab_double(collision_freq)
        self.eng.workspace["iono_grid_parms"] = self._as_matlab_double(
            iono_grid_parms, ensure_row=False
        )

        self.eng.workspace["Bx"] = self._as_matlab_double(Bx)
        self.eng.workspace["By"] = self._as_matlab_double(By)
        self.eng.workspace["Bz"] = self._as_matlab_double(Bz)
        self.eng.workspace["geomag_grid_parms"] = self._as_matlab_double(
            geomag_grid_parms, ensure_row=False
        )

        if ray_state_vec_in is None:
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, iono_en_grid, iono_en_grid_5, ...
                    collision_freq, iono_grid_parms, Bx, By, Bz, geomag_grid_parms);
                """,
                nargout=0,
            )
        else:
            self.eng.workspace["ray_state_vec_in"] = ray_state_vec_in
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, iono_en_grid, iono_en_grid_5, ...
                    collision_freq, iono_grid_parms, Bx, By, Bz, ...
                    geomag_grid_parms, ray_state_vec_in);
                """,
                nargout=0,
            )

        ray_data, ray_path_data, ray_state_vec = self._fetch_struct_array(
            ["ray_data", "ray_path_data", "ray_state_vec"]
        )
        logger.info("PHaRLAP raytrace_3d run completed.")
        return ray_data, ray_path_data, ray_state_vec

    def run_pharlap_3d_sp(
        self,
        origin_lat: float,
        origin_lon: float,
        origin_ht: float,
        elevs: np.ndarray,
        ray_bearings: np.ndarray,
        freqs: np.ndarray,
        rad_earth_m: float,
        iono_en_grid: np.ndarray,
        iono_en_grid_5: np.ndarray,
        collision_freq: np.ndarray,
        iono_grid_parms: np.ndarray,
        Bx: np.ndarray,
        By: np.ndarray,
        Bz: np.ndarray,
        geomag_grid_parms: np.ndarray,
        OX_mode: int = 1,
        nhops: int = 1,
        tol: float | np.ndarray = 1e-7,
        ray_state_vec_in=None,
    ):
        """
        Wrapper for PHaRLAP `raytrace_3d_sp` (spherical Earth formulation).
        """
        logger.info("Running PHaRLAP raytrace_3d_sp...")
        self.eng.eval("close all; clear all; clc;", nargout=0)

        self.eng.workspace["origin_lat"] = float(origin_lat)
        self.eng.workspace["origin_lon"] = float(origin_lon)
        self.eng.workspace["origin_ht"] = float(origin_ht)
        self.eng.workspace["elevs"] = self._as_matlab_double(elevs, ensure_row=True)
        self.eng.workspace["ray_bearings"] = self._as_matlab_double(
            ray_bearings, ensure_row=True
        )
        self.eng.workspace["freqs"] = self._as_matlab_double(freqs, ensure_row=True)
        self.eng.workspace["OX_mode"] = int(OX_mode)
        self.eng.workspace["nhops"] = int(nhops)
        self.eng.workspace["tol"] = self._as_matlab_double(tol, ensure_row=True)
        self.eng.workspace["rad_earth_m"] = float(rad_earth_m)

        self.eng.workspace["iono_en_grid"] = self._as_matlab_double(iono_en_grid)
        self.eng.workspace["iono_en_grid_5"] = self._as_matlab_double(iono_en_grid_5)
        self.eng.workspace["collision_freq"] = self._as_matlab_double(collision_freq)
        self.eng.workspace["iono_grid_parms"] = self._as_matlab_double(
            iono_grid_parms, ensure_row=False
        )

        self.eng.workspace["Bx"] = self._as_matlab_double(Bx)
        self.eng.workspace["By"] = self._as_matlab_double(By)
        self.eng.workspace["Bz"] = self._as_matlab_double(Bz)
        self.eng.workspace["geomag_grid_parms"] = self._as_matlab_double(
            geomag_grid_parms, ensure_row=False
        )

        if ray_state_vec_in is None:
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d_sp( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, rad_earth_m, iono_en_grid, ...
                    iono_en_grid_5, collision_freq, iono_grid_parms, ...
                    Bx, By, Bz, geomag_grid_parms);
                """,
                nargout=0,
            )
        else:
            self.eng.workspace["ray_state_vec_in"] = ray_state_vec_in
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d_sp( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, rad_earth_m, iono_en_grid, ...
                    iono_en_grid_5, collision_freq, iono_grid_parms, ...
                    Bx, By, Bz, geomag_grid_parms, ray_state_vec_in);
                """,
                nargout=0,
            )

        ray_data, ray_path_data, ray_state_vec = self._fetch_struct_array(
            ["ray_data", "ray_path_data", "ray_state_vec"]
        )
        logger.info("PHaRLAP raytrace_3d_sp run completed.")
        return ray_data, ray_path_data, ray_state_vec

run_pharlap(ne_grid, collision_freq, elevs, rb, freqs, irreg, nhops=1, tol=1e-07, radius_earth=6371, irregs_flag=0, start_height=50, height_inc=1, range_inc=1)

Run PHaRLAP 2D spherical ray-tracing (raytrace_2d_sp).

Parameters
np.ndarray

2D electron density grid, shape (nalt, nrange), in cm⁻³.

np.ndarray

2D collision frequency grid, same shape as ne_grid, in Hz.

np.ndarray, 1D

Launch elevation angles [°].

float

Ground range of the ionospheric profile endpoint [km].

np.ndarray, 1D

Operating frequencies [MHz].

np.ndarray

Irregularity array (passed to PHaRLAP; shape depends on version).

int

Number of ionospheric hops. Default 1.

float

ODE solver tolerance. Default 1e-7.

float

Earth radius [km]. Default 6371.

int

Irregularity flag (0 = no irregularities). Default 0.

int

Starting height for the ray [km]. Default 50.

height_inc, range_inc : int Grid increment in height and range [km]. Defaults 1.

Returns
SimpleNamespace

Per-ray scalar parameters (group_range, phase_path, …).

SimpleNamespace

Per-ray path arrays (group_range, height, …).

Source code in hfpytrace/pharlap.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def run_pharlap(
    self,
    ne_grid: np.ndarray,
    collision_freq: np.ndarray,
    elevs: np.ndarray,
    rb: float,
    freqs: np.ndarray,
    irreg: np.ndarray,
    nhops: int = 1,
    tol: float = 1e-7,
    radius_earth: float = 6371,
    irregs_flag: int = 0,
    start_height: int = 50,
    height_inc: int = 1,
    range_inc: int = 1,
):
    """Run PHaRLAP 2D spherical ray-tracing (``raytrace_2d_sp``).

    Parameters
    ----------
    ne_grid : np.ndarray
        2D electron density grid, shape (nalt, nrange), in cm⁻³.
    collision_freq : np.ndarray
        2D collision frequency grid, same shape as ``ne_grid``, in Hz.
    elevs : np.ndarray, 1D
        Launch elevation angles [°].
    rb : float
        Ground range of the ionospheric profile endpoint [km].
    freqs : np.ndarray, 1D
        Operating frequencies [MHz].
    irreg : np.ndarray
        Irregularity array (passed to PHaRLAP; shape depends on version).
    nhops : int
        Number of ionospheric hops.  Default ``1``.
    tol : float
        ODE solver tolerance.  Default ``1e-7``.
    radius_earth : float
        Earth radius [km].  Default ``6371``.
    irregs_flag : int
        Irregularity flag (0 = no irregularities).  Default ``0``.
    start_height : int
        Starting height for the ray [km].  Default ``50``.
    height_inc, range_inc : int
        Grid increment in height and range [km].  Defaults ``1``.

    Returns
    -------
    ray_data : SimpleNamespace
        Per-ray scalar parameters (group_range, phase_path, …).
    ray_path_data : SimpleNamespace
        Per-ray path arrays (group_range, height, …).
    """
    logger.info("Running Pharlap...")
    self.eng.eval("close all; clear all; clc;", nargout=0)

    self.eng.workspace["ne_grid"] = matlab.double(ne_grid.tolist())
    self.eng.workspace["elevs"] = matlab.double(elevs.tolist())
    self.eng.workspace["rb"] = rb
    self.eng.workspace["collision_freq"] = matlab.double(collision_freq.tolist())
    self.eng.workspace["freqs"] = matlab.double(freqs.tolist())
    self.eng.workspace["nhops"] = nhops
    self.eng.workspace["tol"] = tol
    self.eng.workspace["radius_earth"] = radius_earth
    self.eng.workspace["irregs_flag"] = irregs_flag
    self.eng.workspace["start_height"] = start_height
    self.eng.workspace["height_inc"] = height_inc
    self.eng.workspace["range_inc"] = range_inc
    self.eng.workspace["irreg"] = matlab.double(irreg.tolist())

    self.eng.eval(
        """
        [ray_data, ray_path_data] = ...
            raytrace_2d_sp(elevs, rb, freqs, nhops, tol, ...
            radius_earth, irregs_flag, ne_grid, ne_grid, ...
            collision_freq, start_height, height_inc, range_inc, irreg ...
            );
        """,
        nargout=0,
    )
    try:
        ray_data, ray_path_data = (
            self.eng.workspace["ray_data"],
            self.eng.workspace["ray_path_data"],
        )
    except ValueError:
        # MATLAB Engine cannot directly return non-scalar struct arrays.
        self.eng.eval(
            "ray_data_json = jsonencode(ray_data); "
            "ray_path_data_json = jsonencode(ray_path_data);",
            nargout=0,
        )
        ray_data = json.loads(self.eng.workspace["ray_data_json"])
        ray_path_data = json.loads(self.eng.workspace["ray_path_data_json"])
    logger.info("Pharlap run completed.")
    return to_namespace(ray_data), to_namespace(ray_path_data)

run_pharlap_3d(origin_lat, origin_lon, origin_ht, elevs, ray_bearings, freqs, iono_en_grid, iono_en_grid_5, collision_freq, iono_grid_parms, Bx, By, Bz, geomag_grid_parms, OX_mode=1, nhops=1, tol=1e-07, ray_state_vec_in=None)

Wrapper for PHaRLAP raytrace_3d (WGS84 formulation).

Source code in hfpytrace/pharlap.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
def run_pharlap_3d(
    self,
    origin_lat: float,
    origin_lon: float,
    origin_ht: float,
    elevs: np.ndarray,
    ray_bearings: np.ndarray,
    freqs: np.ndarray,
    iono_en_grid: np.ndarray,
    iono_en_grid_5: np.ndarray,
    collision_freq: np.ndarray,
    iono_grid_parms: np.ndarray,
    Bx: np.ndarray,
    By: np.ndarray,
    Bz: np.ndarray,
    geomag_grid_parms: np.ndarray,
    OX_mode: int = 1,
    nhops: int = 1,
    tol: float | np.ndarray = 1e-7,
    ray_state_vec_in=None,
):
    """
    Wrapper for PHaRLAP `raytrace_3d` (WGS84 formulation).
    """
    logger.info("Running PHaRLAP raytrace_3d...")
    self.eng.eval("close all; clear all; clc;", nargout=0)

    self.eng.workspace["origin_lat"] = float(origin_lat)
    self.eng.workspace["origin_lon"] = float(origin_lon)
    self.eng.workspace["origin_ht"] = float(origin_ht)
    self.eng.workspace["elevs"] = self._as_matlab_double(elevs, ensure_row=True)
    self.eng.workspace["ray_bearings"] = self._as_matlab_double(
        ray_bearings, ensure_row=True
    )
    self.eng.workspace["freqs"] = self._as_matlab_double(freqs, ensure_row=True)
    self.eng.workspace["OX_mode"] = int(OX_mode)
    self.eng.workspace["nhops"] = int(nhops)
    self.eng.workspace["tol"] = self._as_matlab_double(tol, ensure_row=True)

    self.eng.workspace["iono_en_grid"] = self._as_matlab_double(iono_en_grid)
    self.eng.workspace["iono_en_grid_5"] = self._as_matlab_double(iono_en_grid_5)
    self.eng.workspace["collision_freq"] = self._as_matlab_double(collision_freq)
    self.eng.workspace["iono_grid_parms"] = self._as_matlab_double(
        iono_grid_parms, ensure_row=False
    )

    self.eng.workspace["Bx"] = self._as_matlab_double(Bx)
    self.eng.workspace["By"] = self._as_matlab_double(By)
    self.eng.workspace["Bz"] = self._as_matlab_double(Bz)
    self.eng.workspace["geomag_grid_parms"] = self._as_matlab_double(
        geomag_grid_parms, ensure_row=False
    )

    if ray_state_vec_in is None:
        self.eng.eval(
            """
            [ray_data, ray_path_data, ray_state_vec] = raytrace_3d( ...
                origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                freqs, OX_mode, nhops, tol, iono_en_grid, iono_en_grid_5, ...
                collision_freq, iono_grid_parms, Bx, By, Bz, geomag_grid_parms);
            """,
            nargout=0,
        )
    else:
        self.eng.workspace["ray_state_vec_in"] = ray_state_vec_in
        self.eng.eval(
            """
            [ray_data, ray_path_data, ray_state_vec] = raytrace_3d( ...
                origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                freqs, OX_mode, nhops, tol, iono_en_grid, iono_en_grid_5, ...
                collision_freq, iono_grid_parms, Bx, By, Bz, ...
                geomag_grid_parms, ray_state_vec_in);
            """,
            nargout=0,
        )

    ray_data, ray_path_data, ray_state_vec = self._fetch_struct_array(
        ["ray_data", "ray_path_data", "ray_state_vec"]
    )
    logger.info("PHaRLAP raytrace_3d run completed.")
    return ray_data, ray_path_data, ray_state_vec

run_pharlap_3d_sp(origin_lat, origin_lon, origin_ht, elevs, ray_bearings, freqs, rad_earth_m, iono_en_grid, iono_en_grid_5, collision_freq, iono_grid_parms, Bx, By, Bz, geomag_grid_parms, OX_mode=1, nhops=1, tol=1e-07, ray_state_vec_in=None)

Wrapper for PHaRLAP raytrace_3d_sp (spherical Earth formulation).

Source code in hfpytrace/pharlap.py 
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
def run_pharlap_3d_sp(
    self,
    origin_lat: float,
    origin_lon: float,
    origin_ht: float,
    elevs: np.ndarray,
    ray_bearings: np.ndarray,
    freqs: np.ndarray,
    rad_earth_m: float,
    iono_en_grid: np.ndarray,
    iono_en_grid_5: np.ndarray,
    collision_freq: np.ndarray,
    iono_grid_parms: np.ndarray,
    Bx: np.ndarray,
    By: np.ndarray,
    Bz: np.ndarray,
    geomag_grid_parms: np.ndarray,
    OX_mode: int = 1,
    nhops: int = 1,
    tol: float | np.ndarray = 1e-7,
    ray_state_vec_in=None,
):
    """
    Wrapper for PHaRLAP `raytrace_3d_sp` (spherical Earth formulation).
    """
    logger.info("Running PHaRLAP raytrace_3d_sp...")
    self.eng.eval("close all; clear all; clc;", nargout=0)

    self.eng.workspace["origin_lat"] = float(origin_lat)
    self.eng.workspace["origin_lon"] = float(origin_lon)
    self.eng.workspace["origin_ht"] = float(origin_ht)
    self.eng.workspace["elevs"] = self._as_matlab_double(elevs, ensure_row=True)
    self.eng.workspace["ray_bearings"] = self._as_matlab_double(
        ray_bearings, ensure_row=True
    )
    self.eng.workspace["freqs"] = self._as_matlab_double(freqs, ensure_row=True)
    self.eng.workspace["OX_mode"] = int(OX_mode)
    self.eng.workspace["nhops"] = int(nhops)
    self.eng.workspace["tol"] = self._as_matlab_double(tol, ensure_row=True)
    self.eng.workspace["rad_earth_m"] = float(rad_earth_m)

    self.eng.workspace["iono_en_grid"] = self._as_matlab_double(iono_en_grid)
    self.eng.workspace["iono_en_grid_5"] = self._as_matlab_double(iono_en_grid_5)
    self.eng.workspace["collision_freq"] = self._as_matlab_double(collision_freq)
    self.eng.workspace["iono_grid_parms"] = self._as_matlab_double(
        iono_grid_parms, ensure_row=False
    )

    self.eng.workspace["Bx"] = self._as_matlab_double(Bx)
    self.eng.workspace["By"] = self._as_matlab_double(By)
    self.eng.workspace["Bz"] = self._as_matlab_double(Bz)
    self.eng.workspace["geomag_grid_parms"] = self._as_matlab_double(
        geomag_grid_parms, ensure_row=False
    )

    if ray_state_vec_in is None:
        self.eng.eval(
            """
            [ray_data, ray_path_data, ray_state_vec] = raytrace_3d_sp( ...
                origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                freqs, OX_mode, nhops, tol, rad_earth_m, iono_en_grid, ...
                iono_en_grid_5, collision_freq, iono_grid_parms, ...
                Bx, By, Bz, geomag_grid_parms);
            """,
            nargout=0,
        )
    else:
        self.eng.workspace["ray_state_vec_in"] = ray_state_vec_in
        self.eng.eval(
            """
            [ray_data, ray_path_data, ray_state_vec] = raytrace_3d_sp( ...
                origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                freqs, OX_mode, nhops, tol, rad_earth_m, iono_en_grid, ...
                iono_en_grid_5, collision_freq, iono_grid_parms, ...
                Bx, By, Bz, geomag_grid_parms, ray_state_vec_in);
            """,
            nargout=0,
        )

    ray_data, ray_path_data, ray_state_vec = self._fetch_struct_array(
        ["ray_data", "ray_path_data", "ray_state_vec"]
    )
    logger.info("PHaRLAP raytrace_3d_sp run completed.")
    return ray_data, ray_path_data, ray_state_vec

get_matlab_pharlap_lib(trace_spec=None, version='4.5.3')

Resolve the PHaRLAP MATLAB library path for a given version tag.

Parameters
Path or None

Override the root search path. None uses :data:PHARLAP_LIB_PATH.

str

PHaRLAP version string, e.g. "4.5.3". The directory pharlap_lib/pharlap_<version> must exist under the resolved root.

Returns

str Absolute path to the PHaRLAP sub-directory.

Raises

FileNotFoundError If the expected PHaRLAP directory does not exist.

Source code in hfpytrace/pharlap.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def get_matlab_pharlap_lib(trace_spec: Path | None = None, version: str = "4.5.3"):
    """Resolve the PHaRLAP MATLAB library path for a given version tag.

    Parameters
    ----------
    trace_spec : Path or None
        Override the root search path.  ``None`` uses :data:`PHARLAP_LIB_PATH`.
    version : str
        PHaRLAP version string, e.g. ``"4.5.3"``.  The directory
        ``pharlap_lib/pharlap_<version>`` must exist under the resolved root.

    Returns
    -------
    str
        Absolute path to the PHaRLAP sub-directory.

    Raises
    ------
    FileNotFoundError
        If the expected PHaRLAP directory does not exist.
    """
    base_path = PHARLAP_LIB_PATH if trace_spec is None else trace_spec
    lib_path = (
        base_path / f"pharlap_{version}"
        if base_path.name == "pharlap_lib"
        else base_path / "pharlap_lib" / f"pharlap_{version}"
    )
    if lib_path.exists():
        logger.info(f"Matlab library path found: {lib_path}")
        return str(lib_path)
    raise FileNotFoundError(f"Matlab library path not found: {lib_path}")

Source Code

hfpytrace/pharlap.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
"""PHaRLAP MATLAB engine wrappers.

Provides Python wrappers around the PHaRLAP (Provision of High-frequency
Raytracing Laboratory for Propagation studies) MATLAB library, enabling
2D and 3D HF ray-tracing from Python via the MATLAB Engine API.

Two 3D formulations are supported:

* ``raytrace_3d`` — WGS84 ellipsoidal Earth model.
* ``raytrace_3d_sp`` — spherical Earth model (used for the SCurve study).

Requires
--------
matlab.engine : MATLAB Engine API for Python (ships with MATLAB installation).
PHaRLAP : MATLAB HF ray-tracing library installed at ``PHARLAP_LIB_PATH``.

Classes / Functions
-------------------
get_matlab_pharlap_lib
    Resolve the PHaRLAP library path for a given version tag.
Engine
    MATLAB engine lifecycle manager with PHaRLAP ray-tracing methods.
"""

import json
from pathlib import Path

import matlab
import matlab.engine
import numpy as np
from loguru import logger

from hfpytrace import PHARLAP_LIB_PATH
from hfpytrace.utils import to_namespace


def get_matlab_pharlap_lib(trace_spec: Path | None = None, version: str = "4.5.3"):
    """Resolve the PHaRLAP MATLAB library path for a given version tag.

    Parameters
    ----------
    trace_spec : Path or None
        Override the root search path.  ``None`` uses :data:`PHARLAP_LIB_PATH`.
    version : str
        PHaRLAP version string, e.g. ``"4.5.3"``.  The directory
        ``pharlap_lib/pharlap_<version>`` must exist under the resolved root.

    Returns
    -------
    str
        Absolute path to the PHaRLAP sub-directory.

    Raises
    ------
    FileNotFoundError
        If the expected PHaRLAP directory does not exist.
    """
    base_path = PHARLAP_LIB_PATH if trace_spec is None else trace_spec
    lib_path = (
        base_path / f"pharlap_{version}"
        if base_path.name == "pharlap_lib"
        else base_path / "pharlap_lib" / f"pharlap_{version}"
    )
    if lib_path.exists():
        logger.info(f"Matlab library path found: {lib_path}")
        return str(lib_path)
    raise FileNotFoundError(f"Matlab library path not found: {lib_path}")


class Engine:
    """MATLAB engine lifecycle manager with PHaRLAP ray-tracing methods.

    Starts a MATLAB engine session, adds the PHaRLAP library to the MATLAB
    path, and exposes Python-callable wrappers for:

    * :meth:`run_pharlap` — 2D spherical ray-tracing (``raytrace_2d_sp``).
    * :meth:`run_pharlap_3d` — 3D WGS84 ray-tracing (``raytrace_3d``).
    * :meth:`run_pharlap_3d_sp` — 3D spherical ray-tracing (``raytrace_3d_sp``).

    Parameters
    ----------
    lib_path : str or None
        Explicit path to the PHaRLAP MATLAB library directory.  ``None``
        resolves the path via :func:`get_matlab_pharlap_lib`.

    Notes
    -----
    Call :meth:`close` (or use as a context manager) to cleanly quit the
    MATLAB session and release memory.
    """

    def __init__(self, lib_path: str = None):
        self.eng = matlab.engine.start_matlab()
        env_path = get_matlab_pharlap_lib() if lib_path is None else lib_path
        self.eng.addpath(self.eng.genpath(env_path), nargout=0)
        logger.info("Matlab engine started and library path added.")
        return

    def close(self):
        logger.info("Closing Matlab engine.")
        self.eng.quit()
        return

    def run_pharlap(
        self,
        ne_grid: np.ndarray,
        collision_freq: np.ndarray,
        elevs: np.ndarray,
        rb: float,
        freqs: np.ndarray,
        irreg: np.ndarray,
        nhops: int = 1,
        tol: float = 1e-7,
        radius_earth: float = 6371,
        irregs_flag: int = 0,
        start_height: int = 50,
        height_inc: int = 1,
        range_inc: int = 1,
    ):
        """Run PHaRLAP 2D spherical ray-tracing (``raytrace_2d_sp``).

        Parameters
        ----------
        ne_grid : np.ndarray
            2D electron density grid, shape (nalt, nrange), in cm⁻³.
        collision_freq : np.ndarray
            2D collision frequency grid, same shape as ``ne_grid``, in Hz.
        elevs : np.ndarray, 1D
            Launch elevation angles [°].
        rb : float
            Ground range of the ionospheric profile endpoint [km].
        freqs : np.ndarray, 1D
            Operating frequencies [MHz].
        irreg : np.ndarray
            Irregularity array (passed to PHaRLAP; shape depends on version).
        nhops : int
            Number of ionospheric hops.  Default ``1``.
        tol : float
            ODE solver tolerance.  Default ``1e-7``.
        radius_earth : float
            Earth radius [km].  Default ``6371``.
        irregs_flag : int
            Irregularity flag (0 = no irregularities).  Default ``0``.
        start_height : int
            Starting height for the ray [km].  Default ``50``.
        height_inc, range_inc : int
            Grid increment in height and range [km].  Defaults ``1``.

        Returns
        -------
        ray_data : SimpleNamespace
            Per-ray scalar parameters (group_range, phase_path, …).
        ray_path_data : SimpleNamespace
            Per-ray path arrays (group_range, height, …).
        """
        logger.info("Running Pharlap...")
        self.eng.eval("close all; clear all; clc;", nargout=0)

        self.eng.workspace["ne_grid"] = matlab.double(ne_grid.tolist())
        self.eng.workspace["elevs"] = matlab.double(elevs.tolist())
        self.eng.workspace["rb"] = rb
        self.eng.workspace["collision_freq"] = matlab.double(collision_freq.tolist())
        self.eng.workspace["freqs"] = matlab.double(freqs.tolist())
        self.eng.workspace["nhops"] = nhops
        self.eng.workspace["tol"] = tol
        self.eng.workspace["radius_earth"] = radius_earth
        self.eng.workspace["irregs_flag"] = irregs_flag
        self.eng.workspace["start_height"] = start_height
        self.eng.workspace["height_inc"] = height_inc
        self.eng.workspace["range_inc"] = range_inc
        self.eng.workspace["irreg"] = matlab.double(irreg.tolist())

        self.eng.eval(
            """
            [ray_data, ray_path_data] = ...
                raytrace_2d_sp(elevs, rb, freqs, nhops, tol, ...
                radius_earth, irregs_flag, ne_grid, ne_grid, ...
                collision_freq, start_height, height_inc, range_inc, irreg ...
                );
            """,
            nargout=0,
        )
        try:
            ray_data, ray_path_data = (
                self.eng.workspace["ray_data"],
                self.eng.workspace["ray_path_data"],
            )
        except ValueError:
            # MATLAB Engine cannot directly return non-scalar struct arrays.
            self.eng.eval(
                "ray_data_json = jsonencode(ray_data); "
                "ray_path_data_json = jsonencode(ray_path_data);",
                nargout=0,
            )
            ray_data = json.loads(self.eng.workspace["ray_data_json"])
            ray_path_data = json.loads(self.eng.workspace["ray_path_data_json"])
        logger.info("Pharlap run completed.")
        return to_namespace(ray_data), to_namespace(ray_path_data)

    @staticmethod
    def _as_matlab_double(arr, ensure_row: bool = False):
        a = np.asarray(arr, dtype=float)
        if a.ndim == 0:
            return float(a)
        if ensure_row and a.ndim == 1:
            a = a.reshape(1, -1)
        return matlab.double(a.tolist())

    def _fetch_struct_array(self, keys: list[str]):
        try:
            out = [self.eng.workspace[k] for k in keys]
        except ValueError:
            # MATLAB Engine cannot directly return non-scalar struct arrays.
            enc = "; ".join([f"{k}_json = jsonencode({k})" for k in keys]) + ";"
            self.eng.eval(enc, nargout=0)
            out = [json.loads(self.eng.workspace[f"{k}_json"]) for k in keys]
        return [to_namespace(v) for v in out]

    def run_pharlap_3d(
        self,
        origin_lat: float,
        origin_lon: float,
        origin_ht: float,
        elevs: np.ndarray,
        ray_bearings: np.ndarray,
        freqs: np.ndarray,
        iono_en_grid: np.ndarray,
        iono_en_grid_5: np.ndarray,
        collision_freq: np.ndarray,
        iono_grid_parms: np.ndarray,
        Bx: np.ndarray,
        By: np.ndarray,
        Bz: np.ndarray,
        geomag_grid_parms: np.ndarray,
        OX_mode: int = 1,
        nhops: int = 1,
        tol: float | np.ndarray = 1e-7,
        ray_state_vec_in=None,
    ):
        """
        Wrapper for PHaRLAP `raytrace_3d` (WGS84 formulation).
        """
        logger.info("Running PHaRLAP raytrace_3d...")
        self.eng.eval("close all; clear all; clc;", nargout=0)

        self.eng.workspace["origin_lat"] = float(origin_lat)
        self.eng.workspace["origin_lon"] = float(origin_lon)
        self.eng.workspace["origin_ht"] = float(origin_ht)
        self.eng.workspace["elevs"] = self._as_matlab_double(elevs, ensure_row=True)
        self.eng.workspace["ray_bearings"] = self._as_matlab_double(
            ray_bearings, ensure_row=True
        )
        self.eng.workspace["freqs"] = self._as_matlab_double(freqs, ensure_row=True)
        self.eng.workspace["OX_mode"] = int(OX_mode)
        self.eng.workspace["nhops"] = int(nhops)
        self.eng.workspace["tol"] = self._as_matlab_double(tol, ensure_row=True)

        self.eng.workspace["iono_en_grid"] = self._as_matlab_double(iono_en_grid)
        self.eng.workspace["iono_en_grid_5"] = self._as_matlab_double(iono_en_grid_5)
        self.eng.workspace["collision_freq"] = self._as_matlab_double(collision_freq)
        self.eng.workspace["iono_grid_parms"] = self._as_matlab_double(
            iono_grid_parms, ensure_row=False
        )

        self.eng.workspace["Bx"] = self._as_matlab_double(Bx)
        self.eng.workspace["By"] = self._as_matlab_double(By)
        self.eng.workspace["Bz"] = self._as_matlab_double(Bz)
        self.eng.workspace["geomag_grid_parms"] = self._as_matlab_double(
            geomag_grid_parms, ensure_row=False
        )

        if ray_state_vec_in is None:
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, iono_en_grid, iono_en_grid_5, ...
                    collision_freq, iono_grid_parms, Bx, By, Bz, geomag_grid_parms);
                """,
                nargout=0,
            )
        else:
            self.eng.workspace["ray_state_vec_in"] = ray_state_vec_in
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, iono_en_grid, iono_en_grid_5, ...
                    collision_freq, iono_grid_parms, Bx, By, Bz, ...
                    geomag_grid_parms, ray_state_vec_in);
                """,
                nargout=0,
            )

        ray_data, ray_path_data, ray_state_vec = self._fetch_struct_array(
            ["ray_data", "ray_path_data", "ray_state_vec"]
        )
        logger.info("PHaRLAP raytrace_3d run completed.")
        return ray_data, ray_path_data, ray_state_vec

    def run_pharlap_3d_sp(
        self,
        origin_lat: float,
        origin_lon: float,
        origin_ht: float,
        elevs: np.ndarray,
        ray_bearings: np.ndarray,
        freqs: np.ndarray,
        rad_earth_m: float,
        iono_en_grid: np.ndarray,
        iono_en_grid_5: np.ndarray,
        collision_freq: np.ndarray,
        iono_grid_parms: np.ndarray,
        Bx: np.ndarray,
        By: np.ndarray,
        Bz: np.ndarray,
        geomag_grid_parms: np.ndarray,
        OX_mode: int = 1,
        nhops: int = 1,
        tol: float | np.ndarray = 1e-7,
        ray_state_vec_in=None,
    ):
        """
        Wrapper for PHaRLAP `raytrace_3d_sp` (spherical Earth formulation).
        """
        logger.info("Running PHaRLAP raytrace_3d_sp...")
        self.eng.eval("close all; clear all; clc;", nargout=0)

        self.eng.workspace["origin_lat"] = float(origin_lat)
        self.eng.workspace["origin_lon"] = float(origin_lon)
        self.eng.workspace["origin_ht"] = float(origin_ht)
        self.eng.workspace["elevs"] = self._as_matlab_double(elevs, ensure_row=True)
        self.eng.workspace["ray_bearings"] = self._as_matlab_double(
            ray_bearings, ensure_row=True
        )
        self.eng.workspace["freqs"] = self._as_matlab_double(freqs, ensure_row=True)
        self.eng.workspace["OX_mode"] = int(OX_mode)
        self.eng.workspace["nhops"] = int(nhops)
        self.eng.workspace["tol"] = self._as_matlab_double(tol, ensure_row=True)
        self.eng.workspace["rad_earth_m"] = float(rad_earth_m)

        self.eng.workspace["iono_en_grid"] = self._as_matlab_double(iono_en_grid)
        self.eng.workspace["iono_en_grid_5"] = self._as_matlab_double(iono_en_grid_5)
        self.eng.workspace["collision_freq"] = self._as_matlab_double(collision_freq)
        self.eng.workspace["iono_grid_parms"] = self._as_matlab_double(
            iono_grid_parms, ensure_row=False
        )

        self.eng.workspace["Bx"] = self._as_matlab_double(Bx)
        self.eng.workspace["By"] = self._as_matlab_double(By)
        self.eng.workspace["Bz"] = self._as_matlab_double(Bz)
        self.eng.workspace["geomag_grid_parms"] = self._as_matlab_double(
            geomag_grid_parms, ensure_row=False
        )

        if ray_state_vec_in is None:
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d_sp( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, rad_earth_m, iono_en_grid, ...
                    iono_en_grid_5, collision_freq, iono_grid_parms, ...
                    Bx, By, Bz, geomag_grid_parms);
                """,
                nargout=0,
            )
        else:
            self.eng.workspace["ray_state_vec_in"] = ray_state_vec_in
            self.eng.eval(
                """
                [ray_data, ray_path_data, ray_state_vec] = raytrace_3d_sp( ...
                    origin_lat, origin_lon, origin_ht, elevs, ray_bearings, ...
                    freqs, OX_mode, nhops, tol, rad_earth_m, iono_en_grid, ...
                    iono_en_grid_5, collision_freq, iono_grid_parms, ...
                    Bx, By, Bz, geomag_grid_parms, ray_state_vec_in);
                """,
                nargout=0,
            )

        ray_data, ray_path_data, ray_state_vec = self._fetch_struct_array(
            ["ray_data", "ray_path_data", "ray_state_vec"]
        )
        logger.info("PHaRLAP raytrace_3d_sp run completed.")
        return ray_data, ray_path_data, ray_state_vec