Model Building

Introduction

For choosing a source and interpreting its port output, start with Sources and ports. For repeated acquisitions, terminal/multimode matrices or plane-wave sweeps, see Reusable studies. This page is the command reference.

gprMax has a choice of two methods for building a model to simulate:

1. A text-based (ASCII) input file, which can be created with any text editor, and uses a series of gprMax commands which begin with the hash character (#). This method is recommended for beginners and those not familiar with Python, and is described in this section of the documentation.

2. A Python API, which includes all the functionality of method 1 as well as several more advanced features. This method is recommended for those who prefer to use Python or need access to specific API-only advanced features, and is documented in the Python API section.

The general syntax of the hash commands is:

#command_name: parameter1 parameter2 parameter3 ...

A command and associated parameters should occupy a single line of the input file, and only one command per line is allowed. Hence, the first character of a line containing a command must be the hash character (#). If the line starts with any other character it is ignored by the program. Therefore, user comments or descriptions can be included in the input file. If a line starts with a hash character (#) the program will expect a valid command. If the name of the command is not correct the program will abandon execution and issue an error message. When a command requires more than one parameter then these should be separated using a white space character.

The order of commands in the input file is not important with the exception of object construction commands.

The pre-solve spatial-resolution check can be explicitly overridden for under-resolved research runs using python -m gprMax model.in --allow-underresolved. This is a command-line option, not a hash command. It preserves warnings and does not disable stability or output-validity checks; see spatial-resolution diagnostic. Zero-amplitude waveforms attached to passive sources are excluded from bandwidth estimation, without removing their sources or receiving outputs.

To describe the commands that can be used in the input file and their parameters the following conventions are used:

• f means a real number which can be entered using either a [.] separating the integral from the decimal part, e.g. 1.5, or in scientific notation, e.g. 15e-1 or 0.15e1.

• i means an integer number.

• c means a single character, e.g. y.

• str means a string of characters with no white spaces in between, e.g sand.

• file means a filename.

• [ ] square brackets are used to indicate optional parameters.

Unless otherwise specified, the SI system of units is used throughout gprMax:

• All parameters associated with simulated space (i.e. size of the model, spatial increments, etc…) should be specified in metres.

• All parameters associated with time (i.e. total simulation time, time instants, etc…) should be specified in seconds.

• All parameters denoting frequency should be specified in Hertz.

• All parameters associated with spatial coordinates in the model should be specified in metres. The origin of the coordinate system (0,0) is at the lower left corner of the model.

It is important to note that gprMax converts spatial and temporal parameters given in metres and seconds to integer values corresponding to FDTD cell coordinates and iteration number respectively. Therefore, rounding to the nearest integer number of the user defined values is performed.

The fundamental spatial and temporal discretization steps are denoted as \(\Delta x\) , \(\Delta y\), \(\Delta z\) and \(\Delta t\) respectively.

The commands have been grouped into six categories:

• Essential - required to run any model, such as the domain size and spatial discretization

• General - provide further control over the model

• Material - used to introduce different materials into the model

• Object construction - used to build geometric shapes with different constitutive parameters

• Source and output - used to place source and output points in the model

• PML - provide advanced customisation and optimisation of the absorbing boundary conditions

Essential commands

Most of the commands are optional but there are some essential commands which are necessary in order to construct any model. For example, none of the media and object commands are necessary to run a model. However, without specifying any objects in the model gprMax will simulate free space (air), which on its own, is not particularly useful for GPR modelling. If you have not specified a command which is essential in order to run a model, for example the size of the model, gprMax will terminate execution and issue an appropriate error message.

The essential commands are:

#domain:

Allows you to specify the size of the model. The syntax of the command is:

#domain: f1 f2 f3

where f1 f2 f3 are the size of the model in the x, y, and z directions respectively. For example, to specify a 500 x 500 x 1000 mm model use #domain: 0.5 0.5 1.0.

For an explicitly declared 2D model, use inf for its invariant axis and select the polarisation with #domain_mode. gprMax resolves inf to the internal Yee-grid thickness required by the chosen mode; it does not create an infinite allocation. For example, the following specifies a model invariant in z:

#domain_mode: TM
#domain: 0.5 0.5 inf

The legacy convention of giving one spatial cell on one axis remains supported and is interpreted as TM mode, but #domain_mode with inf is recommended for new models and is required for TE mode.

#domain_mode:

Selects whether the domain is three-dimensional or uses a two-dimensional TM or TE field reduction. The syntax is:

#domain_mode: str1

str1 is TM, TE, or 3D (case-insensitive). 3D is the default when this command is omitted. For TM or TE, exactly one coordinate of #domain must be inf; that coordinate selects the invariant axis. If inf is used without #domain_mode, gprMax defaults to TM for backwards compatibility.

For example, these declarations select TMz and TEz respectively:

#domain_mode: TM
#domain: 0.5 0.5 inf

#domain_mode: TE
#domain: 0.5 0.5 inf

TM uses one internal cell and TE uses two internal cells on the invariant axis. The physical coordinates on that axis may be written as inf in commands that accept points or bounds; gprMax resolves lower and upper bounds to the appropriate faces of the reduced domain. For a single source or receiver, inf on the invariant axis selects the active interior reference layer (rather than a TE boundary layer whose fields are constrained). On an in-plane axis, -inf selects the lower domain face and inf the upper face.

Note

• A 2D model must have exactly one invariant axis.

• Subgrids, symmetry boundaries, transmission-line sources, and magnetic edges are not currently supported in 2D mode.

• A source polarisation must be one of the active components for the selected plane and mode. gprMax rejects incompatible electric and magnetic sources rather than silently creating a zero source.

#dx_dy_dz:

Allows you to specify the discretization of space in the x , y and z directions respectively (i.e. \(\Delta x\) , \(\Delta y\), \(\Delta z\)). The syntax of the command is:

#dx_dy_dz: f1 f2 f3

where f1 is the spatial step in the x direction (\(\Delta x\)), f2 is the spatial step in the y direction (\(\Delta y\)) and f3 is the spatial step in the z direction (\(\Delta z\)). The spatial discretization controls the maximum permissible time step \(\Delta t\) with which the solution advances in time in order to reach the required simulated time window. The relation between \(\Delta t\) and \(\Delta x\) , \(\Delta y\), \(\Delta z\) is:

\[\Delta t \leq \frac{1}{c\sqrt{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2}}},\]

where \(c\) is the speed of light. In gprMax the equality is used to determine \(\Delta t\) from \(\Delta x\) , \(\Delta y\), and \(\Delta z\). Small values of \(\Delta x\) , \(\Delta y\), and \(\Delta z\) result in small values for \(\Delta t\) which means more iterations in order to reach a given simulated time. However, it is important to note that the smaller the values of \(\Delta x\) , \(\Delta y\), \(\Delta z\) and \(\Delta t\) are the more accurate your model will be. See the Guidance on GPR and general electromagnetic modelling section for tips on choosing a spatial discretisation.

#time_window:

Allows you to specify the total required simulated time. The syntax of the command is:

#time_window: f1

or

#time_window: i1

In the first case the f1 parameter determines the required simulated time in seconds. For example, if you want to simulate a GPR trace of 20 nanoseconds then #time_window: 20e-9 can be used. gprMax will perform the necessary number of iterations in order to reach the required simulated time. Alternatively, if the command is specified with an i1 gprMax will interpret this value as a total number of iterations. Hence the command #time_window: 100 means that 100 iterations will be performed. The number of iterations and the total simulated time window are related by:

\[t_w = \Delta t × N_{it},\]

where \(t_w\) is the time window in seconds, \(\Delta t\) the time step, and \(N_{it}\) the number of iterations. gprMax converts the specified time window in seconds to a number of iterations internally using the aforementioned equation. The result of the division is rounded to the nearest integer.

General commands

#include_file:

Allows you to include commands from a file. It will insert the commands from the specified file at the location where the #include_file command is placed. The syntax of the command is:

#include_file: file1

file1 can be the name of the file containing the commands in the same directory as the input file, or file can be the full path to the file containing the commands (allowing you to specify any location).

Includes may be nested. Each relative include path is resolved against the directory of the file containing that command, not the process working directory. Commands retain their inclusion order. Repeated includes are allowed, but cycles and nesting beyond 100 files are rejected. Python blocks inside included files are not supported; keep legacy Python blocks in the top-level input file or use the Python API.

Study and array-codebook preflight use the same include traversal and path rules as model parsing. #study and #array_codebook must be literal commands in the top-level input or an included file, not printed by a #python block: they determine the run before Python preprocessing occurs. Their CSV/JSON paths remain relative to the top-level input file.

Older development versions searched the working directory first and could silently ignore nested includes. Models relying on that search order should use an absolute path or a path relative to the containing file.

#time_step_stability_factor:

Allows you to alter the value of the time step \(\Delta t\) used by gprMax. gprMax uses the equality in the CFL condition, hence the maximum permissible time step. If a smaller time step is required then the syntax of the command is:

#time_step_stability_factor: f1

where f1 can take values \(0 < \textrm{f1} \leq 1\). Then the actual time step used will be \(\textrm{f1} \times \Delta t\), where \(\Delta t\) is calculated using the equality from the CFL condition.

Models declaring a surface-impedance material automatically cap this factor at 0.99 and report the adjustment in the build log. Smaller user factors are preserved. See Automatic time-step margin for the timing and scope of this automatic margin.

Dispersive-material timestep check

The ordinary spatial CFL limit is not sufficient for every dispersive material. After geometry construction has resolved the component materials, gprMax checks the chosen timestep against each defined dispersive material, including imported and averaged mixtures. All poles of a material are checked together. The check uses each grid’s timestep and spatial steps; an invariant direction in a 2D TM or TE model does not contribute a spatial derivative.

If the check fails, gprMax stops before time stepping. It reports the grid, material, timestep, and failed condition. It does not change the timestep or the material parameters automatically. Where possible, the error gives a smaller timestep that passes the check for that material, and a multiplier to apply to your existing #time_step_stability_factor. For example, if your current factor is 0.8 and the reported multiplier is 0.5, use 0.4. Rerun to check all the materials: another material may require a smaller step. The suggestion is a checked candidate, not a computed maximum stable timestep.

The check follows the current-density recurrence of [GIA2014]. It tests the discrete material response for passivity and applies a conservative dispersive extension of the spatial CFL bound. In particular, for positive lossless Lorentz poles in the principal interval \(0 < 2\pi f_p\Delta t < \pi\), the bound is

\[\frac{c_0^2\Delta t^2}{\mu_{r,\min}} \sum_{j\in\mathcal D}\frac{1}{\Delta_j^2} +\sum_p\Delta\epsilon_{rp}\left[\sec(\pi f_p\Delta t)-1\right] <\epsilon_{r\infty}.\]

Here \(\mathcal D\) contains the varying spatial directions and \(\mu_{r,\min}\) is a common lower permeability bound across the grid’s material definitions. Damped Lorentz, Debye, Drude, and inclusive mixed poles are checked through their discrete exponential representation; the required Drude equivalent conductivity is included. No pole-count-dependent work is added to the field-update loop. MPI ranks exchange setup diagnostics so a failure is reported consistently.

Important

A failed conservative certificate is not, by itself, proof that a particular finite geometry is unstable. The diagnostic distinguishes the material passivity test from the spatial bound. Passing these checks is also not a guarantee of accuracy or of stability of every PML, subgrid, or source-device coupling. Continue to check physical resolution, precision, and convergence.

First check units and material parameters for mistakes. Reducing the timestep preserves the intended physical material; changing pole strengths, frequencies, or damping changes that material and should not be done merely to bypass the check. A positive-growth or non-passive fit may need a corrected material representation rather than an arbitrarily small timestep.

#title:

Allows you to include a title for your model. This title is saved in the output file(s). The syntax of the command is:

#title: str1

where str1 can contain white space characters to separate individual words. The title has to be contained in a single line.

#output_dir:

Allows you to control the directory where output file(s) will be stored. The syntax of the command is:

#output_dir: str1

where str1 can be either the absolute path to the directory for the output file(s) or a path relative to the directory of the input files. The default value is the same as the directory of the input files.

Relative paths are based on the top-level input file, including when this command is declared in an included file. The directory is preserved across --geometry-fixed runs; each model still receives its own numbered filename and snapshot directory. Python API OutputDir paths remain relative to the process working directory.

#omp_threads:

Allows you to control how many OpenMP threads (usually the number of physical CPU cores available) are used when running the model. The most computationally intensive parts of gprMax, which are the FDTD solver loops, have been parallelised using OpenMP which supports multi-platform shared memory multiprocessing. The syntax of the command is:

#omp_threads: i1

where i1 is the number of OpenMP threads to use. If #omp_threads is not specified gprMax will first look to see if the environment variable OMP_NUM_THREADS exists, and if not will detect and use all available physical CPU cores on the machine.

For compatibility with input files created for earlier versions of gprMax, #num_threads: i1 is accepted as an exact alias for #omp_threads: i1. #omp_threads is the preferred name because it distinguishes OpenMP CPU threads from other forms of parallel execution. Do not specify both commands in the same input file.

Material commands

Built-in materials

gprMax has three built-in materials which can be used by specifying their identifiers:

• pec is a perfect electric conductor (PEC), represented by infinite electric conductivity.

• pmc is a perfect magnetic conductor (PMC), represented by infinite magnetic conductivity.

• free_space is free space, with \(\epsilon_r = \mu_r = 1\) and \(\sigma = \sigma_* = 0\).

The identifiers grass and water are reserved for internal use and should not be used unless you intentionally want to change their properties.

#material:

Allows you to introduce a material into the model described by a set of constitutive parameters. The syntax of the command is:

#material: f1 f2 f3 f4 str1

• f1 is the relative permittivity, \(\epsilon_r\)

• f2 is the conductivity (Siemens/metre), \(\sigma\)

• f3 is the relative permeability, \(\mu_r\)

• f4 is the magnetic loss (Ohms/metre), \(\sigma_*\)

• str1 is an identifier for the material. It must not contain +, which is reserved for internally generated material IDs. Use _ instead.

For example #material: 3 0.01 1 0 my_sand creates a material called my_sand which has a relative permittivity (frequency independent) of \(\epsilon_r = 3\), a conductivity of \(\sigma = 0.01\) S/m, and is non-magnetic, i.e. \(\mu_r = 1\) and \(\sigma_* = 0\)

#surface_impedance:

Defines a reusable scalar surface-impedance model for closed volumes that occupy mesh cells. This command defines the boundary response, not the volume geometry. Use its identifier directly as the geometry’s sole material identifier. The command accepts either a constant real resistance or a fitted microwave-metal target from a named preset or bulk conductivity:

#surface_impedance: str1 resistance f1
#surface_impedance: str1 preset str2 f2 f3 [auto|i1] [f4] [c1]
#surface_impedance: str1 conductivity f1 f2 f3 [auto|i1] [f4] [c1]

• str1 is the surface-impedance identifier.

• f1 is the positive surface resistance in ohms, or inf for an exact voxel-face PMC boundary. Use PEC geometry for zero surface impedance. In the conductivity form it is instead the positive bulk conductivity in S/m.

• str2 is aluminium, copper, gold, molybdenum, palladium, silver, tungsten, or zinc. Element symbols and aluminum are also accepted.

• f2 and f3 are the mandatory positive, finite fit-band limits in hertz, with f2 < f3.

• auto tests actual runtime pole counts from 1 through 64 and selects the first whose deterministic local fit is independently certified to meet the tolerance. i1 instead requests exactly that many Foster poles in the same range. The default is auto.

• f4 is the positive maximum relative fit-error tolerance; the default is 2e-3. It controls order selection for auto. With an explicit integer order it is a diagnostic threshold; exceeding it produces a warning without changing the requested order.

• c1 controls full-run fit plotting. n (the default) writes the plot only for --geometry-only; y writes it for geometry-only and full runs. To specify a later optional field, include those before it.

Surface-impedance and bulk-material identifiers must be unique. The __impedance_ prefix is reserved for private runtime materials.

Warning

The finite resistance form is a frequency-independent, purely real idealised boundary condition; it does not represent the causal, dispersive response of a physical conductor. gprMax emits a warning when it is built. Use a fitted metal preset or bulk conductivity over the intended frequency band for physically representative conductor loss.

For example, #surface_impedance: resistive_wall resistance 50 defines a 50 Ohm surface. A copper fit over X band is #surface_impedance: copper_wall preset copper 8e9 12e9 auto 2e-3 n. An explicit conductivity and order can use #surface_impedance: alloy_wall conductivity 3.2e7 8e9 12e9 12.

The metal presets use recommended resistivities for pure bulk metals at 293 K from Matula and Desai et al.; exact per-metal provenance is stored in the output. gprMax fits the good-conductor surface impedance with a non-negative Foster expansion. This makes the resulting realisation passive over the complete frequency axis; its specified fit accuracy applies only inside the selected band. Output HDF5 files store the internal continuous model, fit diagnostics, provenance, model hashes, and the exact local runtime f, q, Z0 recurrence coefficients under surface_impedance_models.

This microwave-only implementation requires f3 <= 300e9 and checks \(\sigma/(\omega\epsilon_0)\geq100\) at the upper band edge. It rejects optical/terahertz bands and targets for which the good-conductor law is not a valid approximation.

See Surface impedance and impedance volumes for model selection, geometry, and troubleshooting; Surface impedance: theory and validation gives the boundary equations and validation.

Note

SIBC supports the 3-D and 2-D TE/TM CPU solvers with cell-occupying geometry. In 2-D, extrude geometry through the full invariant dimension. Pass the #surface_impedance ID directly as the sole material identifier of #box, #sphere, #ellipsoid, #cylinder, #cone, or a finite-thickness #triangle or #cylindrical_sector. Sheets, lines, zero-thickness patches, and directional material assignments are rejected. The implementation does not yet support MPI, subgrids, accelerator backends, or thin wires. PEC/PMC symmetry contacts and supported dispersive retained dielectrics are described in Surface impedance and impedance volumes. Axial discrete plane waves are unsupported. A homogeneous vector/angle plane wave is allowed only when the complete impedance boundary lies strictly inside its TFSF box. The common-metal presets describe thick, homogeneous bulk metal in the local good-conductor regime. They are not optical Drude–Lorentz data and are not valid for thin films, rough surfaces, nanoscale conductors, or thicknesses comparable with the skin depth.

Direct 3-D and 2-D TE/TM #eigenmode_port planes are supported when the impedance boundary is propagation-invariant through the modal plane and the modal window contains its required boundary field components. Every modal window has a PEC rim, with or without a virtual waveguide: tangential electric fields there are zero, replacing any SIBC rows on the rim. Extend the window at least one opaque voxel beyond a wall to retain its surface impedance and loss. The FDFD operator uses the exact time-discrete ADE response and clipped Yee Ampere rows for the retained interior walls. Keep every eigenmode anchor and its trapezoidal bilinear-warped evaluation frequency inside the surface model’s declared fit band; extrapolation is rejected. The bulk operator uses the Yee temporal and longitudinal difference symbols; dispersive bulk poles still use their analytic physical-frequency response. The surface ADE reduction uses its exact discrete response; see Eigenmode ports: theory and numerical methods. #virtual_waveguide supports passive SIBC walls, including fitted metal models, uniformly extruded along the propagation direction. Its transverse PEC boundary matches the modal solve. Opaque-voxel padding keeps SIBC walls active inside the window along each physical transverse axis (only one such axis in 2-D). SIBC can continue through a longitudinal PML when its wall and every retained material are uniformly extruded along the absorption direction, including neighbouring stencil cells. At intersecting edges, retained materials may be heterogeneous, lossy, and electrically dispersive. Each must be isotropic, with positive finite high-frequency permittivity and permeability and finite nonnegative conductivities. The surface model may also be dispersive. Domain PML and CPU virtual guides support these hosts; SIBC coupling remains unavailable on MPI, subgrids, and accelerator backends. See Waveguides and absorbers for the slab-coverage and aperture requirements.

#surface_impedance: wall resistance inf selects an exact PMC on the voxel face. It retains the clipped H circulation and electric dual area with zero surface admittance. This also supports longitudinal PML and virtual waveguides under the same restrictions.

#material_from_database:

Creates a material from an official or local versioned JSON material database. See Material databases for the schema, lookup policy, and provenance information. The syntax is:

#material_from_database: str1 str2 [str3]

• str1 is the database name, without .json.

• str2 is the entry key.

• str3 is an optional local material ID. It defaults to str2 and is useful for giving the material a model-specific name, avoiding a name collision, or matching the name expected by imported geometry. It must not contain +, which is reserved for internally generated material IDs.

For example #material_from_database: fundamental vacuum creates the official vacuum entry and uses vacuum as its local ID. To refer to the same material as model_vacuum within the model, use #material_from_database: fundamental vacuum model_vacuum.

#material_density:

Assigns an optional physical mass density to one or more existing materials. The syntax is:

#material_density: f1 str1 [str2 ...]

• f1 is a finite, positive mass density in \(\mathrm{kg\,m^{-3}}\).

• str1, str2, … are existing material identifiers.

For example, #material_density: 1040 brain white_matter assigns the same mass density to the two named materials. Density remains a cell-centred property of the final volumetric material assignment. It is not included in dielectric smoothing and has no effect on ordinary electromagnetic updates. Materials without this command retain an unspecified density; this is valid for ordinary simulations, but a dosimetry request will reject selected cells whose material has no density.

Note

The Debye, Lorentz, and Drude commands below describe electric dispersion. They cannot be applied to PEC or PMC materials, including user-defined materials with infinite electric or magnetic conductivity.

#add_dispersion_debye:

Allows you to add dispersive properties to an already defined #material based on a multiple pole Debye formulation (see Software Features section). For example, the susceptibility function for a single-pole Debye material is given by:

\[\chi_p (t) = \frac{\Delta \epsilon_{rp}}{\tau_p} e^{-t/\tau_p},\]

where \(\Delta \epsilon_{rp} = \epsilon_{rsp} - \epsilon_{r \infty}\), \(\epsilon_{rsp}\) is the zero-frequency relative permittivity for the pole, \(\epsilon_{r \infty}\) is the relative permittivity at infinite frequency, and \(\tau_p\) is the pole relaxation time.

The syntax of the command is:

#add_dispersion_debye: i1 f1 f2 f3 f4 ... str1

• i1 is the number of Debye poles.

• f1 is the difference between the zero-frequency relative permittivity and the relative permittivity at infinite frequency, i.e. \(\Delta \epsilon_{rp1} = \epsilon_{rsp1} - \epsilon_{r \infty}\) , for the first Debye pole.

• f2 is the relaxation time (seconds), \(\tau_{p1}\), for the first Debye pole.

• f3 is the difference between the zero-frequency relative permittivity and the relative permittivity at infinite frequency, i.e. \(\Delta \epsilon_{rp2} = \epsilon_{rsp2} - \epsilon_{r \infty}\) , for the second Debye pole.

• f4 is the relaxation time (seconds), \(\tau_{p2}\), for the second Debye pole.

• …

• str1 identifies the material to add the dispersive properties to.

For example to create a model of water with a single Debye pole, \(\epsilon_{rsp1} = 80.1\), \(\epsilon_{r \infty} = 4.9\) and \(\tau_{p1} = 9.231\times 10^{-12}\) seconds use: #material: 4.9 0 1 0 my_water and #add_dispersion_debye: 1 75.2 9.231e-12 my_water.

Note

• You can continue to add pairs of values for \(\Delta \epsilon_{rp}\) and \(\tau_p\) for as many Debye poles as you have specified with i1.

• The relative permittivity in the #material command should be given as the relative permittivity at infinite frequency, i.e. \(\epsilon_{r \infty}\).

• Temporal values associated with pole frequencies and relaxation times should always be greater than the time step \(\Delta t\) used in the model.

• Very long relaxation times can also require double precision: if \(\exp(-\Delta t/\tau_p)\) rounds to one in the selected precision, gprMax warns that long-time pole decay may be inaccurate. The small exponential differences used for pole coupling are evaluated with cancellation-resistant expm1 arithmetic, but the stored recurrence still has the solver’s precision. Check time-window and precision convergence when this warning occurs. The same consideration applies to the decaying exponential poles of Lorentz, Drude, and mixed materials.

#add_dispersion_lorentz:

Allows you to add dispersive properties to an already defined #material based on a multiple pole Lorentz formulation (see Software Features section). For example, the susceptibility function for a single-pole Lorentz material is given by:

\[\chi_p (t) = \Re \left\{ -j\gamma_p e^{(-\delta_p + j\beta_p)t} \right\},\]

where

\[\beta_p = \sqrt{(2\pi f_p)^2 - \delta_p^2} \quad \textrm{and} \quad \gamma_p = \frac{(2\pi f_p)^2 \Delta \epsilon_{rp}}{\beta_p},\]

where \(\Delta \epsilon_{rp} = \epsilon_{rsp} - \epsilon_{r \infty}\), \(\epsilon_{rsp}\) is the zero-frequency relative permittivity for the pole, \(\epsilon_{r \infty}\) is the relative permittivity at infinite frequency, \(f_p\) is the resonance frequency (Hertz) of the pole pair, \(\delta_p\) is the damping coefficient (per second), and \(j=\sqrt{-1}\). The angular resonance frequency used in the formulation is \(\omega_p=2\pi f_p\).

The syntax of the command is:

#add_dispersion_lorentz: i1 f1 f2 f3 f4 f5 f6 ... str1

• i1 is the number of Lorentz poles.

• f1 is the difference between the zero-frequency relative permittivity and the relative permittivity at infinite frequency, i.e. \(\Delta \epsilon_{rp1} = \epsilon_{rsp1} - \epsilon_{r \infty}\) , for the first Lorentz pole.

• f2 is the resonance frequency (Hertz), \(f_{p1}\), for the first Lorentz pole.

• f3 is the damping coefficient (per second), \(\delta_{p1}\), for the first Lorentz pole.

• f4 is the difference between the zero-frequency relative permittivity and the relative permittivity at infinite frequency, i.e. \(\Delta \epsilon_{rp2} = \epsilon_{rsp2} - \epsilon_{r \infty}\) , for the second Lorentz pole.

• f5 is the resonance frequency (Hertz), \(f_{p2}\), for the second Lorentz pole.

• f6 is the damping coefficient (per second), \(\delta_{p2}\), for the second Lorentz pole.

• …

• str1 identifies the material to add the dispersive properties to.

Note

• You can continue to add triplets of values for \(\Delta \epsilon_{rp}\), \(f_p\) and \(\delta_p\) for as many Lorentz poles as you have specified with i1.

• The relative permittivity in the #material command should be given as the relative permittivity at infinite frequency, i.e. \(\epsilon_{r \infty}\).

• The recursive formulation requires \(f_p < 1/\Delta t\), \(\delta_p < 1/\Delta t\), and the underdamped condition \(\delta_p < 2\pi f_p\). These coefficient limits do not replace the stricter Nyquist and spatial-resolution limits on useful simulation output.

• These per-pole input limits do not establish coupled-update stability. The complete material must also pass the Dispersive-material timestep check.

#add_dispersion_drude:

Allows you to add dispersive properties to an already defined #material based on a multiple pole Drude formulation (see Software Features section). For example, the susceptibility function for a single-pole Drude material is given by:

\[\chi_p (t) = \frac{(2\pi f_p)^2}{\gamma_p} (1-e^{-\gamma_p t}),\]

where \(f_p\) is the plasma frequency (Hertz) of the pole, \(\omega_p=2\pi f_p\) is its angular frequency, and \(\gamma_p\) is the inverse of the pole relaxation time (per second).

The syntax of the command is:

#add_dispersion_drude: i1 f1 f2 f3 f4 ... str1

• i1 is the number of Drude poles.

• f1 is the plasma frequency (Hertz), \(f_{p1}\), for the first Drude pole.

• f2 is the inverse of the relaxation time (per second), \(\gamma_{p1}\), for the first Drude pole.

• f3 is the plasma frequency (Hertz), \(f_{p2}\), for the second Drude pole.

• f4 is the inverse of the relaxation time (per second), \(\gamma_{p2}\) for the second Drude pole.

• …

• str1 identifies the material to add the dispersive properties to.

Note

• You can continue to add pairs of values for \(f_p\) and \(\gamma_p\) for as many Drude poles as you have specified with i1.

• The recursive formulation requires \(f_p < 1/\Delta t\) and \(\gamma_p < 1/\Delta t\). These coefficient limits do not replace the stricter Nyquist and spatial-resolution limits on useful simulation output.

• The complete material must also pass the Dispersive-material timestep check.

#material_range:

Allows you to create a series of materials with properties specified by ranges of relative permittivity, conductivity, relative permeability, and magnetic loss. The command is designed to be used in conjunction with the #fractal_box command for spatial distributions of dielectric properties. The syntax of the command is:

#material_range: f1 f2 f3 f4 f5 f6 f7 f8 str1

• f1 is the lower end of the range of relative permittivity values.

• f2 is the upper end of the range of relative permittivity values.

• f3 is the lower end of the range of conductivity values.

• f4 is the upper end of the range of conductivity values.

• f5 is the lower end of the range of relative permeability values.

• f6 is the upper end of the range of relative permeability values.

• f7 is the lower end of the range of magnetic loss values.

• f8 is the upper end of the range of magnetic loss values.

• str1 is an identifier for the material range.

The number of bins is set by #fractal_box. Each bin uses the midpoint of the corresponding interval in each property range. Material reuse compares the full computed values, not rounded display values: small differences in conductivity or other properties remain distinct. Bins with exactly identical properties can share a generated material, including across material ranges. Generated material names are internal identifiers; use str1 to select the range rather than constructing those names yourself.

For example to create a series of 10 materials with relative permittivity ranging between 2 and 6, \(\sigma=0\), \(\mu_r=1\), and \(\sigma_*=0\), distributed using a fractal approach, use: #material_range: 2 6 0 0 1 1 0 0 er2_6 and #fractal_box: 0 0 0 0.15 0.15 0.15 1.5 1 1 1 10 er2_6 my_frac_box.

#material_list:

Allows you to create a list of pre-defined materials that can be used in conjunction with the #fractal_box command for spatial distributions of dielectric properties. The syntax of the command is:

#material_list: str1 str2 ... str3

• str1 and str2 are identifiers for materials. You can have identifiers for as many pre-defined materials as required.

• str3 is an identifier for the material list.

For example to create a fractal distribution of two different sand materials and water use: #material: 3 0 1 0 sand1, #material: 4 0.1 1 0 sand2, #material: 4.9 0.001 1 0 my_water, #add_dispersion_debye: 1 75.2 9.231e-12 my_water, #material_list: sand1 sand2 my_water my_list, #fractal_box: 0 0 0 0.15 0.15 0.15 1.5 1 1 1 3 my_list my_frac_box.

#material_crim:

Allows you to use the Complex Refractive Index Model (CRIM) to mix a fixed-fraction non-dispersive matrix material with a single-pole Debye dispersive material (e.g. water or brine), with the remaining volume fraction assumed to be air. The command is designed to be used in conjunction with the #fractal_box command for creating media, such as wet soils, wet concrete, or brine-bearing ice, with realistic dielectric and geometric properties. Both the matrix and dispersive materials must already exist, i.e. be defined using #material (and, for the dispersive material, #add_dispersion_debye with exactly one pole) before this command is used.

CRIM is a general power-law mixing formula for the bulk complex permittivity of a multi-phase medium [BIR1974], widely used for estimating the permittivity of soils, snow, ice, and other granular or porous media from the properties and volume fractions of their constituents. It takes the general form

\[\varepsilon_{\mathrm{mix}}(\omega)^{a} = \sum_{i} f_i \, \varepsilon_i(\omega)^{a},\]

where \(f_i\) and \(\varepsilon_i(\omega)\) are the volumetric fraction and complex relative permittivity of constituent \(i\) (with \(\sum_i f_i = 1\)), and \(a\) is an empirical shape factor that encodes the geometry of the mixed phases relative to the applied field: \(a=1\) reduces to a simple volumetric (linear) average of permittivities, \(a=0.5\) is the square-root, or “refractive index”, form most commonly used for soils and other granular GPR media, and as \(a \to 0\) the model approaches logarithmic (Lichtenecker) mixing.

The syntax of the command is:

#material_crim: str1 f1 str2 f2 f3 f4 f5 f6 str3

• str1 is an identifier for an existing non-dispersive material used for the fixed-fraction matrix (solid) phase.

• f1 is the fixed volumetric fraction of the matrix phase.

• str2 is an identifier for an existing single-pole Debye material used for the dispersive phase.

• f2 and f3 define a range for the volumetric fraction of the dispersive phase.

• f4 and f5 are the lower and upper bounds of the frequency range (Hz) used to fit the CRIM mixing curve.

• f6 is the CRIM shape factor, \(a\), described above.

• str3 is an identifier for the CRIM mixing model.

For this command, the mixture has exactly three phases - the fixed-fraction matrix, the dispersive phase, and air - so \(\varepsilon_i\) is specialised to

\[\varepsilon_{\mathrm{mix}}(\omega)^{a} = f_1 \, \varepsilon_{\mathrm{matrix}}^{a} + f_{\mathrm{disp}} \left[\varepsilon_{\infty} + \frac{\Delta\varepsilon}{1+j\omega\tau}\right]^{a} + f_{\mathrm{air}} \, (1)^{a},\]

where \(\varepsilon_{\infty}\), \(\Delta\varepsilon\), and \(\tau\) are the dispersive material’s own Debye parameters (taken directly from its #add_dispersion_debye definition), \(f_{\mathrm{disp}}\) is a value between f2 and f3 for each bin, and the implied air fraction is \(f_{\mathrm{air}} = 1 - f_1 - f_{\mathrm{disp}}\) (so f1 plus f3 must not exceed 1). For each of the n_materials bins created by the associated #fractal_box command, a single-pole Debye material is fitted by linear least-squares over the frequency range f4 to f5 to this dielectric CRIM curve. The dispersive material’s relaxation time \(\tau\) is fixed and reused for every bin; the fit determines the pole weight and infinite-frequency permittivity.

The constituent DC conductivities are not included inside the fractional power in the equation above. They are mixed separately using the volumetric rule

\[\sigma_{\mathrm{mix}} = f_1 \sigma_{\mathrm{matrix}} + f_{\mathrm{disp}} \sigma_{\mathrm{disp}},\]

with zero conductivity assumed for the air phase. This is a defined approximation rather than a fit of the complete conductive complex permittivity. The current implementation is limited to non-magnetic constituents (\(\mu_r=1\) and \(\sigma_m=0\)).

For example, to create a fractal distribution of wet sand, with a fixed sand fraction of 0.6 and a water fraction ranging from 0.02 to 0.35, fitted over 1 MHz to 3 GHz with the standard CRIM shape factor, use: #material: 5 0 1 0 sand, #material: 4.9 0 1 0 water, #add_dispersion_debye: 1 73.3389 8.0994e-12 water, #material_crim: sand 0.6 water 0.02 0.35 1e6 3e9 0.5 wetsand, and then specify the fractal box using #fractal_box: 0 0 0 0.1 0.1 0.08 1.5 1 1 1 10 wetsand my_fractal_box.

Note

Only passive, single-pole Debye materials are accepted for the dispersive phase - CRIM’s dispersive constituent is conventionally assumed to be water or brine, which is well represented by a single relaxation. Using a material with more than one Debye pole, a non-Debye dispersive material, a perfect conductor, or a magnetic material will raise an error.

#soil_peplinski:

Allows you to use a mixing model for soils proposed by Peplinski (http://dx.doi.org/10.1109/36.387598), valid for frequencies in the range 0.3GHz to 1.3GHz. The command is designed to be used in conjunction with the #fractal_box command for creating soils with realistic dielectric and geometric properties. The syntax of the command is:

#soil_peplinski: f1 f2 f3 f4 f5 f6 str1

• f1 is the sand fraction of the soil.

• f2 is the clay fraction of the soil.

• f3 is the bulk density of the soil in grams per centimetre cubed.

• f4 is the density of the sand particles in the soil in grams per centimetre cubed.

• f5 and f6 define a range for the volumetric water fraction of the soil.

• str1 is an identifier for the soil.

For example for a soil with sand fraction 0.5, clay fraction 0.5, bulk density \(2~g/cm^3\), sand particle density of \(2.66~g/cm^3\), and a volumetric water fraction range of 0.001 - 0.25 use: #soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil.

Note

Further information on the Peplinski soil model and our implementation can be found in ‘Giannakis, I. (2016). Realistic numerical modelling of Ground Penetrating Radar for landmine detection. The University of Edinburgh, United Kingdom. (http://hdl.handle.net/1842/20449)’

Object construction commands

Object construction commands are processed in the order they appear in the input file. Therefore space in the model allocated to a specific material using for example the #box command can be reallocated to another material using the same or any other object construction command. Space in the model can be regarded as a canvas in which objects are introduced and one can be overlaid on top of the other overwriting its properties in order to produce the desired geometry. The object construction commands can therefore be used to create complex shapes and configurations.

A #surface_impedance ID can be used directly wherever these commands take one isotropic material identifier for a closed, cell-occupying volume. For example:

#surface_impedance: resistive_wall resistance 50
#box: 0.1 0.1 0.1 0.3 0.3 0.2 resistive_wall n

Normal ordered overwrite semantics apply. Several objects with the same surface-impedance ID can form a union; a later ordinary material object can carve a cavity; and any later geometry can replace earlier cells. Plates, edges, zero-thickness triangles and cylindrical sectors, and boxes that occupy zero cells on any axis reject surface-impedance IDs. A scalar surface impedance also cannot be used as one member of a directional anisotropic material assignment. #fractal_box accepts a surface-impedance ID only as its sole material with one material bin and a roughness, grass, or water modifier; an unmodified one-material fractal box must use #box instead.

The final rasterized geometry is checked after all ordered overwrites. Face-connected staircases are valid. Impedance cells that touch only diagonally across a Yee edge or only at a grid vertex are non-manifold and are rejected, including contacts between different surface-impedance IDs. Within each local 2 x 2 x 2 vertex neighbourhood, both the impedance and retained cell sets must be face-connected when non-empty. Connect the regions through a full voxel face, or separate them so they share neither an edge nor a vertex; mesh refinement, increased feature thickness, or a small position change can also repair the final voxel topology. See Surface impedance and impedance volumes for the complete rule and troubleshooting guidance.

Anisotropy

It is possible to specify objects that have diagonal anisotropy which allows materials such as wood and fibre-reinforced composites, often imaged with GPR, to be more accurately modelled.

\[\begin{split}\bar{\bar{\epsilon}} = \left[ \begin{array}{ccc} \epsilon_{xx} & 0 & 0 \\ 0 & \epsilon_{yy} & 0 \\ 0 & 0 & \epsilon_{zz} \end{array} \right],\quad \bar{\bar{\sigma}}= \left[ \begin{array}{ccc} \sigma_{xx} & 0 & 0 \\ 0 & \sigma_{yy} & 0 \\ 0 & 0 & \sigma_{zz} \end{array} \right]\end{split}\]

Standard isotropic objects specify one material identifier that defines the same properties in x, y, and z directions. However, every volumetric object building command can also be specified with three material identifiers, which allows properties for the x, y, and z directions to be separately defined. The #plate command, which defines a surface, can specify up to two material identifiers, and the #edge command, which defines a line, continues to take one material identifier. For example to create a box with different material properties in each of the x, y, and z directions use:

#material: 41 10 1 0 matX
#material: 35 10 1 0 matY
#material: 33 1 1 0 matZ
#box: 0 0 0 0.1 0.1 0.1 matX matY matZ

As another example, to create a cylinder of radius 10 mm that has the same properties in the x and y directions but different properties in the z direction use:

#material: 41 10 1 0 matXY
#material: 33 1 1 0 matZ
#cylinder: 0.1 0.1 0.1 0.5 0.1 0.1 0.01 matXY matXY matZ

Dielectric smoothing

At the boundaries between different materials in the model there is the question of which material properties to use. Should the last object to be defined at that location dictate the properties? Should an average set of properties of the materials of the objects that share that location be used? This latter option is often referred to as dielectric smoothing and has been shown to result in more accurate simulations [LUE1994] [BOU1996]. To address this question gprMax includes an option to turn dielectric smoothing on or off for volumetric object building commands. The default behaviour (if no option is specified) is for dielectric smoothing to be on. The option can be specified with a single character y (on) or n (off) given after the material identifier in each object command. For example to specify a sphere of material sand with dielectric smoothing turned off use: #sphere: 0.5 0.5 0.5 0.1 sand n.

Note

• Dispersive materials are not smoothed by default. Use #dispersive_averaging: y to include Debye, Lorentz, and Drude media in electric-edge smoothing.

• If an object is anistropic then dielectric smoothing is automatically turned off for that object.

• Non-volumetric object building commands, #edge, #plate, and #triangle (applies to triangular patch not triangular prism) cannot have dielectric smoothing.

#magnetic_averaging:

Selects the mixing rule used for magnetic-field components at smoothed material interfaces. Each H component is constructed from the two cells stacked along its own axis. Because the normal component of magnetic flux density is continuous across an interface, the harmonic mean of relative permeability (\(\mu_r\)) and magnetic loss (\(\sigma_*\)) is used by default. Electric-field smoothing is unchanged and continues to use its arithmetic four-cell average. The syntax is:

#magnetic_averaging: str1

• str1 is harmonic or arithmetic (case-insensitive).

Note

• This command is optional; the default is harmonic.

• Earlier versions of gprMax used an arithmetic magnetic average. #magnetic_averaging: arithmetic selects that earlier mixing rule, but does not undo source-timing, magnetic-component placement or PEC-interface corrections. It does not guarantee exact reproduction of v3 results; see 5. Revalidate numerical reference results.

• The command chooses the magnetic mixing rule only; it does not enable or disable dielectric smoothing.

#dispersive_averaging:

Enables or disables interface averaging for dispersive materials. When enabled, Debye, Lorentz, and Drude media use the same arithmetic four-cell electric-edge average as nondispersive dielectrics. The Debye case follows the contour-path formulation developed by [HAR2020]; its extension to all three dispersion families uses the inclusive susceptibility representation of [GIA2014]. The syntax is:

#dispersive_averaging: c1

• c1 is y to enable dispersive averaging or n to disable it.

For surrounding-cell weights \(w_m\), the effective high-frequency permittivity and conductivity are

\[\begin{split}\begin{aligned} \epsilon_{\infty,\mathrm{eff}} &= \sum_m w_m\epsilon_{\infty,m}, \\ \sigma_{\mathrm{eff}} &= \sum_m w_m\sigma_m. \end{aligned}\end{split}\]

Writing a constituent susceptibility in inclusive pole-residue form as

\[\chi_m(\omega) = \sum_p \frac{W_{m,p}}{j\omega-Q_{m,p}},\]

where a complex Lorentz term implicitly includes its conjugate contribution, the effective response is

\[\epsilon_{r,\mathrm{eff}}(\omega) = \sum_m w_m\left[ \epsilon_{\infty,m}+ \frac{\sigma_m}{j\omega\epsilon_0}+ \chi_m(\omega) \right].\]

Consequently, every pole location \(Q_{m,p}\) is retained and its residue \(W_{m,p}\) is multiplied by the corresponding cell weight. For a Debye term this is equivalent to retaining its relaxation time and scaling \(\Delta\epsilon\):

\[\chi_{m,p}(\omega) = \frac{w_m\Delta\epsilon_{m,p}} {1+j\omega\tau_{m,p}}.\]

The four surrounding cells each contribute \(1/4\); repeated materials therefore naturally produce weights of \(1/2\) or \(3/4\). Terms with identical pole locations are combined exactly. Two different single-pole materials generally produce a two-pole effective material, even if each bulk material uses only one pole.

Note

• This command is optional; the default is n. Enabling it is recommended when interface accuracy is more important than the possible increase in memory use.

• Nondispersive dielectric smoothing remains controlled by the optional y/n argument of each volumetric geometry command.

• Dispersive state arrays use the model-wide maximum pole count for computational efficiency. An averaged interface containing additional distinct pole locations can therefore increase memory use throughout the grid. The resolved maximum is included in gprMax’s memory estimate.

• The production solver preserves the exact pole union. It does not apply automatic pole reduction because any reduced model is valid only over a chosen frequency band and error tolerance.

#geometry_view:

Allows you to output file(s) containing information about the geometry of a model. The file(s) use the open source Visualization ToolKit (VTK) VTKHDF format and require ParaView 5.12 or newer for viewing. The command can be used to create several 3D views of the model which are useful for checking that it has been constructed as desired. The syntax of the command is:

#geometry_view: f1 f2 f3 f4 f5 f6 f7 f8 f9 file1 c1

• f1 f2 f3 are the lower left (x,y,z) coordinates of the volume of the geometry view in metres.

• f4 f5 f6 are the upper right (x,y,z) coordinates of the volume of the geometry view in metres.

• f7 f8 f9 are the spatial discretisation of the geometry view in metres. Typically these will be the same as the spatial discretisation of the model but they can be courser if desired.

• file1 is the filename of the file where the geometry view will be stored in the same directory as the input file.

• c1 can be either n (normal) or f (fine) which specifies whether to output the geometry information on a per-cell basis (n) or a per-cell-edge basis (f). The fine mode should be reserved for viewing detailed parts of the geometry that occupy small volumes, as using this mode can generate geometry files with large file sizes.

Tip

When you want to just check the geometry of your model, run gprMax using the optional command line argument --geometry-only. This will build the model and produce any geometry view files, but will not run the simulation.

#edge:

Allows you to introduce a wire with specific properties into the model. A wire is an edge of a Yee cell and it can be useful to model resistors or thin wires. The syntax of the command is:

#edge: f1 f2 f3 f4 f5 f6 str1

• f1 f2 f3 are the starting (x,y,z) coordinates of the edge, and f4 f5 f6 are the ending (x,y,z) coordinates of the edge. The coordinates should define a single line.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

For example to specify a x-directed wire that is a perfect electric conductor, use: #edge: 0.5 0.5 0.5 0.7 0.5 0.5 pec. Note that the y and z coordinates are identical.

#thin_wire:

Allows a conducting wire whose physical radius is smaller than the Yee cell to be represented by an axis-aligned thin-wire model. The logarithmic radius correction is based on Umashankar, Taflove, and Beker [UMA1987], with the improved electric/magnetic contour factors proposed by Mäkinen, Juntunen, and Kivikoski [MAK2002]. The wire occupies electric edges like a PEC #edge, while the magnetic updates on the four surrounding Yee edges are corrected to represent the specified sub-cell radius. The syntax is:

#thin_wire: f1 f2 f3 f4 f5 f6 f7

• f1 f2 f3 and f4 f5 f6 are the start and end coordinates of one non-zero, axis-aligned wire in metres.

• f7 is the physical wire radius \(a\) in metres. It must be positive and smaller than half the minimum transverse cell size.

For a wire along \(w\), consider a surrounding \(H_v\) edge whose radial direction is \(u\) (so \(u\), \(v\), and \(w\) are the three Cartesian axes). The Umashankar radius factor is

\[F_u = \frac{2}{\ln(\Delta u/a)},\]

and the Mäkinen contour factors are

\[k_{H_v} = \frac{\Delta u}{\Delta v} \tan^{-1}\!\left(\frac{\Delta v}{\Delta u}\right), \qquad k_{E_u} = \frac{1}{k_{H_v}}.\]

gprMax stores the projected Yee-edge value \(\widetilde{H}_v=k_{H_v}H_v\). In that representation, Mäkinen’s magnetic update is implemented by multiplying the radial curl coefficient by \(F_u k_{H_v}\). The coefficient for the derivative along the wire remains the background value because \(k_{H_v}k_{E_u}=1\). The ordinary electric update then consumes \(\widetilde{H}_v\) directly, thereby applying the required \(k_H\) correction to the radial electric field without a special runtime kernel. On a square transverse mesh, \(k_H=\pi/4\) and \(k_E=4/\pi\).

The magnetic permeability and magnetic conductivity at each affected H component are inherited from its already-resolved background material, including magnetic material averaging. The magnetic-source coefficient is also multiplied by \(k_H\); therefore a co-located #magnetic_frill_source must not apply \(k_H\) a second time. It does, however, apply the feed-cell factor \(F_u\), as required by Hyun’s discrete magnetic-current equation.

Important

Receiver samples on the four H edges immediately surrounding the wire are the stored projected values \(\widetilde{H}\). Divide such a sample by its orientation-specific \(k_H\) if the unprojected point value is required. Current loops and ordinary electric updates should use the stored values directly.

For example, a z-directed wire of radius 0.1 mm is specified by: #thin_wire: 0.05 0.05 0.02 0.05 0.05 0.12 0.0001.

The wire may lie on a transverse domain face only when that face is a PMC symmetry boundary. A wire and its surrounding magnetic stencil must not touch a PML region. MPI domain decomposition is supported: each electric edge and each member of the surrounding magnetic stencil is constructed on its owning rank, including a stencil that crosses an internal rank boundary. Thin wires in 2D models and overlapping sub-cell wire junctions are not currently supported. The charge-based end-cap treatment from [MAK2002] is not implemented, so an isolated open end retains the usual staircasing/electrically-long end error; the improved straight-section update is used up to the final wire edge. No special runtime solver is used: the corrected material coefficients are consumed by the normal CPU, CUDA, OpenCL, and Metal field-update kernels.

#magnetic_edge:

Allows you to introduce a single magnetic-field edge with specific properties into the model. It is the magnetic dual of #edge. The syntax is:

#magnetic_edge: f1 f2 f3 f4 f5 f6 str1

• f1 f2 f3 are the starting (x,y,z) coordinates of the edge, and f4 f5 f6 are its ending coordinates. The coordinates must define a single axis-aligned line.

• str1 is a material identifier that has already been defined, or one of the built-in materials.

For example, an x-directed perfect magnetic conductor is specified with #magnetic_edge: 0.5 0.5 0.5 0.7 0.5 0.5 pmc.

Note

#magnetic_edge is not currently supported in 2D mode.

#plate:

Allows you to introduce a plate with specific properties into the model. A plate is a surface of a Yee cell and it can be useful to model objects thinner than a Yee cell. The syntax of the command is:

#plate: f1 f2 f3 f4 f5 f6 str1

• f1 f2 f3 are the lower left (x,y,z) coordinates of the plate, and f4 f5 f6 are the upper right (x,y,z) coordinates of the plate. The coordinates should define a surface and not a 3D object like the #box command.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

For example to specify a xy oriented plate that is a perfect electric conductor, use: #plate: 0.5 0.5 0.5 0.7 0.8 0.5 pec. Note that the z coordinates are identical.

#triangle:

The volumetric commands #triangle (non-zero thickness), #box, #sphere, #cylinder, #cylindrical_sector, #cone, #ellipsoid, and #fractal_box accept an optional final positional tag_name string. Geometry tags are cell-centred semantic metadata and are independent of materials. They follow normal ordered geometry semantics: a later tagged volume replaces the older tag in its cells, while a later untagged volume clears those cells to tag ID zero. Consequently #cylinder: ... plastic n container followed by a smaller untagged free-space cylinder creates a tagged hollow shell. The tag must always be the last argument. To keep the positional grammar unambiguous and backward compatible, hash-command users must explicitly provide the preceding y or n smoothing argument whenever a tag is present. A tagged #fractal_box must therefore provide its seed and smoothing arguments before the tag. Tags may contain letters, digits, _, -, ., and :; spaces are not permitted. The Python API does not require an explicit smoothing argument when tag is used.

Directional volumes use the three ordered x/y/z material IDs followed by the smoothing marker and tag, for example:

#box: 0.01 0.01 0.01 0.02 0.02 0.02 tissue_x tissue_y tissue_z n body

Directional assignments remain non-averaged even if the marker is y. Their tags can be selected by SAR and radiometry outputs, which use the loss of each directional material with the corresponding electric component.

Allows you to introduce a triangular patch or a triangular prism with specific properties into the model. The patch is just a triangular surface made as a collection of staircased Yee cells, and the triangular prism extends the triangular patch in the direction perpendicular to the plane. The syntax of the command is:

#triangle: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 str1 [c1 [tag_name]]

• f1 f2 f3 are the coordinates (x,y,z) of the first apex of the triangle, f4 f5 f6 the coordinates (x,y,z) of the second apex, and f7 f8 f9 the coordinates (x,y,z) of the third apex.

• f10 is the thickness of the triangular prism. If the thickness is zero then a triangular patch is created.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing. For use only when creating a triangular prism, not a triangular patch.

For example, to specify a xy orientated triangular patch that is a perfect electric conductor, use: #triangle: 0.5 0.5 0.5 0.6 0.4 0.5 0.7 0.9 0.5 0.0 pec. Note that the z coordinates are identical and the thickness is zero.

#box:

Allows you to introduce an orthogonal parallelepiped with specific properties into the model. The syntax of the command is:

#box: f1 f2 f3 f4 f5 f6 str1 [c1 [tag_name]]

• f1 f2 f3 are the lower left (x,y,z) coordinates of the parallelepiped, and f4 f5 f6 are the upper right (x,y,z) coordinates of the parallelepiped.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing.

#sphere:

Allows you to introduce a spherical object with specific parameters into the model. The syntax of the command is:

#sphere: f1 f2 f3 f4 str1 [c1 [tag_name]]

• f1 f2 f3 are the coordinates (x,y,z) of the centre of the sphere.

• f4 is its radius.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing.

For example, to specify a sphere with centre at (0.5, 0.5, 0.5), radius 100 mm, and with constitutive parameters of my_sand, use: #sphere: 0.5 0.5 0.5 0.1 my_sand.

Note

• Sphere objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created.

#cylinder:

Allows you to introduce a circular cylinder into the model. The orientation of the cylinder axis can be arbitrary, i.e. it does not have align with one of the Cartesian axes of the model. The syntax of the command is:

#cylinder: f1 f2 f3 f4 f5 f6 f7 str1 [c1 [tag_name]]

• f1 f2 f3 are the coordinates (x,y,z) of the centre of one face of the cylinder, and f4 f5 f6 are the coordinates (x,y,z) of the centre of the other face.

• f7 is the radius of the cylinder.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing.

For example, to specify a cylinder with its axis in the y direction, a length of 0.7 m, a radius of 100 mm, and that is a perfect electric conductor, use: #cylinder: 0.5 0.1 0.5 0.5 0.8 0.5 0.1 pec.

Note

• Cylinder objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created.

#cylindrical_sector:

Allows you to introduce a cylindrical sector (shaped like a slice of pie) into the model. The syntax of the command is:

#cylindrical_sector: c1 f1 f2 f3 f4 f5 f6 f7 str1 [c1 [tag_name]]

• c1 is the direction of the axis of the cylinder from which the sector is defined and can be x, y, or z.

• f1 f2 are the coordinates of the centre of the cylindrical sector.

• f3 f4 are the lower and higher coordinates of the axis of the cylinder from which the sector is defined (in effect they specify the thickness of the sector).

• f5 is the radius of the cylindrical sector.

• f6 is the starting angle (in degrees) for the cylindrical sector (with zero degrees defined on the positive first axis of the plane of the cylindrical sector).

• f7 is the angle (in degrees) swept by the cylindrical sector (the finishing angle of the sector is always anti-clockwise from the starting angle).

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing.

For example, to specify a cylindrical sector with its axis in the z direction, radius of 0.25 m, thickness of 2 mm, a starting angle of 330 \(^\circ\), a sector angle of 60 \(^\circ\), and that is a perfect electric conductor, use: #cylindrical_sector: z 0.34 0.24 0.500 0.502 0.25 330 60 pec.

Note

• Cylindrical sector objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created.

#cone:

Allows you to introduce a cone into the model. The orientation of the cylinder axis can be arbitrary, i.e. it does not have align with one of the Cartesian axes of the model. The syntax of the command is:

#cone: f1 f2 f3 f4 f5 f6 f7 f8 str1 [c1 [tag_name]]

• f1 f2 f3 are the coordinates (x,y,z) of the centre of the first face of the cone, and f4 f5 f6 are the coordinates (x,y,z) of the centre of the other face.

• f7 is the radius of the first face of the cone, and f8 is the radius of the other face of the cone.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing.

For example, to specify a cone with centres at (0.02, 0.075, 0.075) and (0.08, 0.075, 0.075), radii of 30 mm and 0 mm, and with constitutive parameters of my_sand, use: cone: 0.02 0.075 0.075 0.08 0.075 0.075 0.03 0 my_sand.

Note

• Cone objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created.

#ellipsoid:

Allows you to introduce an ellipsoid into the model. The syntax of the command is:

#ellipsoid: f1 f2 f3 f4 f5 f6 str1 [c1 [tag_name]]

• f1 f2 f3 are the coordinates (x,y,z) of the centre of the ellipsoid.

• f4 f5 f6 are the coordinates (x,y,z) of the semi-axes of the ellipsoid.

• str1 is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing.

For example, to specify an ellipsoid with centre at (0.045, 0.045, 0.045), and semi-axes (0.03, 0.02, 0.03), and with constitutive parameters of my_sand, use: #ellipsoid: 0.045 0.045 0.045 0.03 0.02 0.03 my_sand.

Note

• Ellipsoidal objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created.

#fractal_box:

Allows you to introduce an orthogonal parallelepiped with fractal distributed properties which are related to a mixing model or normal material into the model. The syntax of the command is:

#fractal_box: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 i1 str1 str2 [i2 [c1 [tag_name]]]

• f1 f2 f3 are the lower left (x,y,z) coordinates of the parallelepiped, and f4 f5 f6 are the upper right (x,y,z) coordinates of the parallelepiped.

• f7 is the fractal dimension which, for an orthogonal parallelepiped, should take values between zero and three.

• f8 is used to weight the fractal in the x direction.

• f9 is used to weight the fractal in the y direction.

• f10 is used to weight the fractal in the z direction.

• i1 is the number of materials to use for the fractal distribution (defined according to the associated mixing model). This should be set to one if using a normal material instead of a mixing model.

• str1 is an identifier for the associated mixing model or material.

• str2 is an identifier for the fractal box itself.

• i2 is an optional parameter which controls the seeding of the random number generator used to create the fractals. By default (if you don’t specify this parameter) the random number generator will be seeded by trying to read data from /dev/urandom (or the Windows analogue) if available or from the clock otherwise.

• c1 is an optional parameter which can be y or n, used to switch on and off dielectric smoothing. If c1 is specified then a value for i2 must also be present.

For example, to create an orthogonal parallelepiped with fractal distributed properties using a Peplinski mixing model for soil, with 50 different materials over a range of water volumetric fractions from 0.001 - 0.25, you should first define the mixing model using: #soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil and then specify the fractal box using #fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_fractal_box.

Several fractal boxes may reference the same mixing model with different bin counts. Each box retains the bin-to-material mapping calculated for its own i1 value, including when surface modifiers are applied. Reusing a mixing model does not change the materials assigned to previously prepared boxes. For a material list, each box uses the first i1 entries of the list; i1 must not exceed the list length.

Note

• We are not aware of a formulation of Perfectly Matched Layer (PML) absorbing boundary that can specifically handle distributions of material properties (such as those created by fractals) throughout the thickness of the PML, i.e. this is a required area of research. Our PML formulations can work to an extent depending on your modelling scenario and requirements. You may need to increase the thickness of the PML and/or consider tuning the parameters of the PML (PML functions) to improve performance for your specific model.

#add_surface_roughness:

Allows you to add rough surfaces to a #fractal_box in the model. A fractal distribution is used for the profile of the rough surface. The syntax of the command is:

#add_surface_roughness: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 str1 [i1]

• f1 f2 f3 are the lower left (x,y,z) coordinates of a surface on a #fractal_box, and f4 f5 f6 are the upper right (x,y,z) coordinates of a surface on a #fractal_box. The coordinates must locate one of the six surfaces of a #fractal_box but do not have to extend over the entire surface.

• f7 is the fractal dimension which, for an orthogonal parallelepiped, should take values between zero and three.

• f8 is used to weight the fractal in the first direction of the surface.

• f9 is used to weight the fractal in the second direction of the surface.

• f10 f11 define lower and upper limits for a range over which the roughness can vary. These limits should be specified relative to the dimensions of the #fractal_box that the rough surface is being applied.

• str1 is an identifier for the #fractal_box that the rough surface should be applied to.

• i1 is an optional parameter which controls the seeding of the random number generator used to create the fractals. By default (if you don’t specify this parameter) the random number generator will be seeded by trying to read data from /dev/urandom (or the Windows analogue) if available or from the clock otherwise.

Up to six #add_rough_surface commands can be given for any #fractal_box corresponding to the six surfaces.

For example, if a #fractal_box has been specified using: #fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_fractal_box then to apply a rough surface that varies between 85 mm and 110 mm (i.e. valleys that are up to 15 mm deep and peaks that are up to 10 mm tall) to the surface that is in the positive z direction, use #add_surface_roughness: 0 0 0.1 0.1 0.1 0.1 1.5 1 1 0.085 0.110 my_fractal_box.

#add_surface_water:

Allows you to add surface water to a #fractal_box in the model that has had a rough surface applied. The syntax of the command is:

#add_surface_water: f1 f2 f3 f4 f5 f6 f7 str1

• f1 f2 f3 are the lower left (x,y,z) coordinates of a surface on a #fractal_box, and f4 f5 f6 are the upper right (x,y,z) coordinates of a surface on a #fractal_box. The coordinates must locate one of the six surfaces of a #fractal_box but do not have to extend over the entire surface.

• f7 defines the depth of the water, which should be specified relative to the dimensions of the #fractal_box that the surface water is being applied.

• str1 is an identifier for the #fractal_box that the surface water should be applied to.

For example, to add surface water that is 5 mm deep to an existing #fractal_box that has been specified using #fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_fractal_box and has had a rough surface applied using #add_surface_roughness: 0 0 0.1 0.1 0.1 0.1 1.5 1 1 0.085 0.110 my_fractal_box, use #add_surface_water: 0 0 0.1 0.1 0.1 0.1 0.105 my_fractal_box.

Note

• The water is modelled using a single-pole Debye formulation with properties \(\epsilon_{rs} = 80.1\), \(\epsilon_{\infty} = 4.9\), and a relaxation time of \(\tau = 9.231 \times 10^{-12}\) seconds (http://dx.doi.org/10.1109/TGRS.2006.873208). If you prefer, gprMax will use your own definition for water as long as it is named water.

#add_grass:

Allows you to add grass with roots to a #fractal_box in the model. The blades of grass are randomly distributed over the specified surface area and a fractal distribution is used to vary the height of the blades of grass and depth of the grass roots. The syntax of the command is:

#add_grass: f1 f2 f3 f4 f5 f6 f7 f8 f9 i1 str1 [i2]

• f1 f2 f3 are the lower left (x,y,z) coordinates of a surface on a #fractal_box, and f4 f5 f6 are the upper right (x,y,z) coordinates of a surface on a #fractal_box. The coordinates must locate one of three surfaces (in the positive axis direction) of a #fractal_box but do not have to extend over the entire surface.

• f7 is the fractal dimension which, for an orthogonal parallelepiped, should take values between zero and three.

• f8 f9 define lower and upper limits for a range over which the height of the blades of grass can vary. These limits should be specified relative to the dimensions of the #fractal_box that the grass is being applied.

• i1 is the number of blades of grass that should be applied to the surface area.

• str1 is an identifier for the #fractal_box that the grass should be applied to.

• i2 is an optional parameter which controls the seeding of the random number generator used to create the fractals. By default (if you don’t specify this parameter) the random number generator will be seeded by trying to read data from /dev/urandom (or the Windows analogue) if available or from the clock otherwise.

For example, to apply 100 blades of grass that vary in height between 100 and 150 mm to the entire surface in the positive z direction of a #fractal_box that had been specified using #fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 50 my_soil my_fractal_box, use #add_grass: 0 0 0.1 0.1 0.1 0.1 1.5 0.2 0.25 100 my_fractal_box.

Note

• The grass is modelled using a single-pole Debye formulation with properties \(\epsilon_{rs} = 18.5087\), \(\epsilon_{\infty} = 12.7174\), and a relaxation time of \(\tau = 1.0793 \times 10^{-11}\) seconds (http://dx.doi.org/10.1007/BF00902994). If you prefer, gprMax will use your own definition for grass if you use a material named grass. The geometry of the blades of grass are defined by the parametric equations: \(x = x_c +s_x {\left( \frac{t}{b_x} \right)}^2\), \(y = y_c +s_y {\left( \frac{t}{b_y} \right)}^2\), and \(z=t\), where \(s_x\) and \(s_y\) can be -1 or 1 which are randomly chosen, and where the constants \(b_x\) and \(b_y\) are random numbers based on a Gaussian distribution.

#geometry_objects_read:

Allows you to insert pre-defined geometry into a model. New geometry files store compact integer arrays in HDF5, stable entry names in /material_keys, and constitutive properties in a companion JSON material database. The syntax of the command is:

#geometry_objects_read: f1 f2 f3 file1 file2 [y/n]

• f1 f2 f3 are the lower left (x,y,z) coordinates in the domain where the lower left corner of the geometry array should be placed.

• file1 is the path to and filename of the HDF5 file that contains an integer array which defines the geometry.

• file2 is the material database name without .json. It is resolved beside file1.

• The optional y/n parameter controls dielectric interface averaging when reconstructing a voxel-only file. The default is n. Use y for voxelised tissue or other dielectric models when smoothed interfaces are desired.

Warning

Geometry files and simulation results both use .h5. For a single run, model.in produces model.h5 by default. If the geometry input is that same file, saving the simulation results overwrites it without confirmation. Importing the geometry does not protect the original file from later output writes.

Use distinct paths for geometry and results: for example, name the geometry file model_geometry.h5, choose a separate output directory, or set a different result basename with -o on the command line or outputfile in the Python API. The same precaution applies when choosing a geometry-export filename.

Note

• Store the cell-material indices in an integer array at the root named data. The gprMax writer uses np.int16; the reader also accepts wider integer arrays and legacy unsigned arrays. Non-negative values index the file’s material table, not the destination grid’s material IDs. Imported indices are mapped without narrowing the destination IDs to 16 bits.

• /material_keys maps each non-negative integer to an entry key in the JSON database.

• In a signed array, -1 means leave the existing geometry unchanged; it does not mean free space. Other negative indices are invalid. For complete component meshes, transparency in /ID also preserves the corresponding existing rigidity flags. Cell materials and tags follow /data independently, allowing explicit component-only objects without overwriting the surrounding cells.

• The spatial resolution of the geometry objects must match the spatial resolution defined in the model.

• The spatial resolution must be specified as a root attribute of the HDF5 file with the name dx_dy_dz equal to a tuple of floats, e.g. (0.002, 0.002, 0.002)

• Legacy .txt material files are not accepted by the simulator. First convert the original HDF5/text pair using the Migrating legacy geometry and text materials utility, then supply the converted HDF5 file and its JSON database name. Do not simply rename a text file.

• Averaging is applied only when the imported file contains cell material indices but no complete Yee-component arrays. A file containing /ID, /rigidE, and /rigidH records an authoritative component mesh, so the optional averaging flag is ignored.

• For imported dispersive dielectrics, also enable #dispersive_averaging: y. Without that global option, dispersive materials remain non-smoothable even if this command requests averaging. Cell-centred material density and geometry tags are never averaged.

For example, after converting a downloaded 2 mm AustinMan model to the current HDF5/JSON format as described in the AustinMan/AustinWoman toolbox, insert it with its lower-left corner 40 mm from the domain origin and enable tissue-interface averaging using #geometry_objects_read: 0.04 0.04 0.04 AustinMan_v2.6_2x2x2_gprmax.h5 AustinMan_v2_6_2mm_materials y.

#geometry_objects_write:

Allows you to write geometry generated in a model to file. The file can be read back into gprMax using the #geometry_objects_read command. This allows complex geometry that can take some time to generate to be saved to file and more quickly imported into subsequent models. Geometry arrays, cell-centred semantic tags (when present), their ID-to-name table, and stable material keys are saved in HDF5. Corresponding material definitions are saved in a versioned _materials.json database. The syntax of the command is:

#geometry_objects_write: f1 f2 f3 f4 f5 f6 file1

• f1 f2 f3 are the lower left (x,y,z) coordinates of the parallelepiped, and f4 f5 f6 are the upper right (x,y,z) coordinates of the parallelepiped.

• file1 is the basename for the files where geometry and material information will be stored.

Note

• The structure of the HDF5 file is the same as that described for the #geometry_objects_read command.

• Objects are stored using spatial resolution defined in the model.

• The exported material table includes both cell-centred materials and Yee-edge materials, including generated interface averages. Indices are compacted before writing; the int16 geometry format permits at most 32768 distinct exported materials. Larger catalogues raise an error rather than wrapping material indices. Large global material IDs are allowed when the exported subset fits this limit.

Source and output commands

#waveform:

Allows you to specify common waveform shapes to use with sources in the model. The syntax of the command is:

#waveform: str1 f1 f2 str2

• str1 is the type of waveform which can be:

  • gaussian which is a Gaussian waveform.

  • gauspulse which is a Gaussian-modulated cosine with fractional bandwidth 0.5 at -6 dB, matching MATLAB’s default gauspuls bandwidth. It is shifted to start at the -60 dB envelope level; see gauspulse (modulated Gaussian) for its delay, formula and plots.

  • gaussiandot which is the first derivative of a Gaussian waveform.

  • gaussiandotnorm which is the normalised first derivative of a Gaussian waveform.

  • gaussiandotdot which is the second derivative of a Gaussian waveform.

  • gaussiandotdotnorm which is the normalised second derivative of a Gaussian waveform.

  • ricker which is a Ricker (or Mexican hat) waveform, i.e. the negative, normalised second derivative of a Gaussian waveform.

  • gaussianprime which is the first derivative of a Gaussian waveform, directly derived from the aforementioned gaussian (see notes below).

  • gaussiandoubleprime which is the second derivative of a Gaussian waveform, directly derived from the aforementioned gaussian (see notes below).

  • sine which is a single cycle of a sine waveform.

  • contsine which is a continuous sine waveform with a linear amplitude ramp over the first four cycles.

  • impulse which is one nonzero sample on the source’s update lattice, not a continuous-time unit-area Dirac delta. The required positive f2 value is unused by this waveform; use 1. A hard source continues to clamp its edge to zero after the impulse while its start/stop window remains active; see #voltage_source:.

• f1 is the scaling of the maximum amplitude of the waveform (for a #hertzian_dipole the units will be Amps, for a #voltage_source or #transmission_line the units will be Volts).

• f2 is the centre frequency of the waveform (Hertz). In the case of the Gaussian waveform it is related to the pulse width.

• str2 is an identifier for the waveform used to assign it to a source.

For example, to specify the normalised first derivative of a Gaussian waveform with an amplitude of one and a centre frequency of 1.2GHz, use: #waveform: gaussiandotnorm 1 1.2e9 my_gauss_pulse.

For a Gaussian-modulated cosine at 1 GHz, use #waveform: gauspulse 1 1e9 pulse. Its intrinsic peak delay is about 2.781 ns; allow about 5.6 ns for the pulse itself, plus the required response-recording time. The source’s optional start time adds to that delay. For another bandwidth, use a user-defined waveform through the Python API or #excitation_file; the built-in command has no bandwidth parameter.

Note

• Waveforms are implemented in gprMax/waveforms.py. The gprMax/toolboxes/Plotting package provides waveform visualisation; see plot_source_wave.py.

• gaussiandot, gaussiandotnorm, gaussiandotdot, gaussiandotdotnorm, ricker waveforms have their centre frequencies specified by the user, i.e. they are not derived to the ‘base’ gaussian

• gaussianprime and gaussiandoubleprime waveforms are the first derivative and second derivative of the ‘base’ gaussian waveform, i.e. the centre frequencies of the waveforms will rise for the first and second derivatives.

#excitation_file:

Allows you to specify an ASCII file that contains amplitude values that specify custom waveform(s) that can be used with sources in the model.

The first row of each column must begin with a identifier string that will be used as the name of each waveform. Subsequent rows should contain amplitude values for the custom waveform you want to use. You can import multiple different waveforms (as columns of amplitude data) in a single file.

Ideally, there should be the same number of amplitude values as number of iterations in your model. If there are less amplitude values than the number of iterations in the model, the end of the sequence of amplitude values will be padded with zero values up to the number of iterations. If extra amplitude values are specified than needed then they are ignored.

Optionally, in the first column of the file you may specify your own time vector of values (which must use the identifier time) to use with the amplitude values of the waveform.

The amplitude values will be interpolated using either the aforementioned user specified time vector, or if none was supplied, a vector of time values corresponding to the simulation time step and number of iterations will be used. Key parameters used for the interpolation can be specified in the command.

The syntax of the command is:

#excitation_file: file1 [str1 str2]

• file1 can be the name of the file containing the specified waveform in the same directory as the input file, or file can be the full path to the file containing the specified waveform (allowing you to specify any location).

• str1 and str2 are an optional parameter pair that allow values for kind and fill_value to be passed to the interpolation function (scipy.interpolate.interp1d). If omitted, gprMax uses linear interpolation and zero outside the supplied time axis. A numeric fill value applies at both ends; extrapolation occurs only when explicitly requested with linear extrapolate (or another interpolation kind).

Without a time column, the time axis is exactly arange(iterations) * dt. Electric-current sources (resistive voltage sources and Hertzian dipoles) retain their half-step sampling and can request a sample past the final supplied whole-step value; the selected fill policy applies there. Magnetic dipoles retain their whole-step sampling. A zero-resistance voltage source evaluates only its required whole-step electric-field samples. These rules do not change source start/stop gating.

For example, to specify the file my_waves.txt, which contains two custom waveform shapes, use: #excitation_file: my_waves.txt. The contents of the file my_waves.txt would take the form:

time my_pulse1 my_pulse2
0 0 0
1.926e-12 1.2e-6 0
3.852e-12 1.3e-6 1.0e-1
5.778e-12 5.0e-6 1.5e-1
...       ...    ...
...       ...    ...
...       ...    ...

Then to use my_pulse1 custom waveform shape with, for example, a z-polarised Hertzian dipole source:

#hertzian_dipole: z 0.5 0.5 0.5 my_pulse1

Note

• The #waveform command is not necessary when using a custom waveform excitation, only the #excitation_file command and whatever source is going to be used with the custom waveform excitation.

The modulated Gaussian example includes the built-in gauspulse, a custom-bandwidth callable, and a plotting script that can export a table for #excitation_file.

#hertzian_dipole:

Outputs (no circuit port)

This ideal impressed-current source does not create a circuit port, S11 or input impedance. Add Rx / #rx for field or Ampere-loop observations, or an NTFF request for transformed fields. Raw field/source ratios are not power-wave S-parameters. See Sources and ports for source selection and GPR acquisitions for point-source acquisition schedules.

Allows you to specify a current density term at an electric field location - the simplest excitation, often referred to as an additive or soft source.

Source positions must satisfy the physical Yee-component bounds, including after #src_steps or study updates. Valid terminal components are retained for symmetry-boundary excitations; padded array entries outside the physical component are rejected.

\[J_s = \frac{I \Delta l}{\Delta x \Delta y \Delta z},\]

where \(J_s\) is the current density, \(I\) is the current, \(\Delta l\) is the length of the infinitesimal electric dipole, and \(\Delta x\), \(\Delta y\), and \(\Delta z\) are the spatial resolution of the grid.

Note

• \(\Delta l\) is set equal to \(\Delta x\), \(\Delta y\), or \(\Delta z\) depending on the specified polarisation.

The syntax of the command is:

#hertzian_dipole: c1 f1 f2 f3 str1 [f4 f5]

• c1 is the polarisation of the source and can be x, y, or z.

• f1 f2 f3 are the coordinates (x,y,z) of the source in the model.

• f4 f5 are optional parameters. f4 is a time delay in starting the source. f5 is a time to remove the source. If the time window is longer than the source removal time then the source will stop after the source removal time. If the source removal time is longer than the time window then the source will be active for the entire time window. If f4 f5 are omitted the source will start at the beginning of time window and stop at the end of the time window.

• str1 is the identifier of the waveform that should be used with the source.

For example, to use a x-polarised Hertzian dipole with unit amplitude and a 600 MHz centre frequency Ricker waveform, use: #waveform: ricker 1 600e6 my_ricker_pulse and #hertzian_dipole: x 0.05 0.05 0.05 my_ricker_pulse.

Note

• When a #hertzian_dipole is used in a 2D simulation it acts as a line source of current in the invariant (geometry) direction of the simulation.

#magnetic_dipole:

Outputs (no circuit port)

This ideal impressed-current source does not create a circuit port, S11 or input impedance. Add Rx / #rx for field or Ampere-loop observations, or an NTFF request for transformed fields. Raw field/source ratios are not power-wave S-parameters. See Sources and ports for source selection and GPR acquisitions for point-source acquisition schedules.

This will simulate an infinitesimal magnetic dipole. This is often referred to as an additive or soft source.

The physical Yee-component bounds for magnetic components differ from those for electric components: a magnetic dipole is half-cell shifted in its two transverse directions, not its own axis.

The syntax of the command is:

#magnetic_dipole: c1 f1 f2 f3 str1 [f4 f5]

• c1 is the polarisation of the source and can be x, y, or z.

• f1 f2 f3 are the coordinates (x,y,z) of the source in the model.

• f4 f5 are optional parameters. f4 is a time delay in starting the source. f5 is a time to remove the source. If the time window is longer than the source removal time then the source will stop after the source removal time. If the source removal time is longer than the time window then the source will be active for the entire time window. If f4 f5 are omitted the source will start at the beginning of time window and stop at the end of the time window.

• str1 is the identifier of the waveform that should be used with the source.

#rational_network:, #network_terminal:, #network_excitation:

Port and outputs

Output is explicit: add #network_port: feed or scene.add(gprMax.NetworkPort(terminal_id='feed')) for terminal feed. The network and its excitation alone do not request port results. Output is /ports/feed, or /subgrids/<subgrid ID>/ports/feed for a CPU or CUDA HSG port.

time, Vgenerator, Vtotal and Inetwork describe the generator, gap voltage (V) and network current (A). Terminal-current spectra remove the background Yee-gap admittance and are distinct from the external-network current. frequency, S11, Zin and Yin have independent validity masks. Reference impedance defaults to 50 Ohms; it does not replace the physical network admittance.

Keep the NetworkPort object to read port.result.total_voltage, port.result.s11, port.result.zin and the other result arrays after solving. Omitting NetworkExcitation leaves a passive measurable load: its received voltage/current remain useful, but source-normalised S11 is undefined. Impedance/admittance are available where their own masks permit.

Use python -m gprMax.toolboxes.Plotting.plot_port model.h5 --port feed --save. See Rational-network port output for the schema, Fixed-topology terminal-source studies for fixed-network drive studies and Ports, antenna power and NTFF for NTFF use.

These commands connect a linear one-port network to one electric Yee edge. The reusable network model is expressed as a rational driving-point admittance

\[Y(s) = G + sC + \sum_{m=1}^{M}\frac{r_m}{s-p_m},\]

where \(G\) is the direct conductance, \(C\) is the direct capacitance, and \(p_m\) and \(r_m\) are pole-residue pairs. The syntax is

#rational_network: str1 f1 f2 i1 [f3 f4 f5 f6 ...] [c1]
#network_terminal: c2 f7 f8 f9 str1 str2
#network_excitation: str2 str3 [f10 f11]

• str1 is the reusable network-model ID.

• f1 and f2 are \(G\) in siemens and \(C\) in farads.

• i1 is the number of poles. Every pole then has four values: pole_real pole_imag residue_real residue_imag. Poles are in rad/s and residues are in S/s. Non-real poles and residues must occur in conjugate pairs so that the time-domain current is real.

• c1 is an optional y/n flag permitting an active model. It is n by default. Passive models are checked over the FDTD band and all poles must be stable.

• c2 is the terminal polarisation x, y, or z; f7 f8 f9 are its coordinates; and str2 is its unique terminal ID.

• str3 is an existing waveform ID. The optional f10 f11 pair gives its start and stop times. Omitting #network_excitation creates a passive load or receiving terminal.

For common elements, a resistor uses \(G=1/R\), a capacitor uses the direct \(C\) term, and an inductor uses \(p=0\), \(r=1/L\). A series \(RL\) branch has \(p=-R/L\) and \(r=1/L\). For example, a driven 50 Ohm terminal is

#waveform: ricker 1 2e9 pulse
#rational_network: source50 0.02 0 0
#network_terminal: z 0.05 0.05 0.02 source50 feed
#network_excitation: feed pulse
#network_port: feed 50 10

The optional excitation is a Thévenin open-circuit voltage applied through the complete rational impedance. The arbitrary linear circuit-to-edge formulation follows [PER1999] and [CHE2007]. gprMax improves the underlying classic PLRC time placement using the exponential recursive-convolution approach of Giannakis and Giannopoulos [GIA2014]. In particular, every pole state is evaluated analytically at the electric half-step under a linearly varying terminal voltage; it is not estimated by averaging its two integer- time values. Each terminal stores only its own poles and applies one locally implicit edge correction, rather than allocating dispersive state throughout the mesh.

A driven network with \(G=1/R\), no direct capacitance, and no poles is discretely identical to a finite-resistance #voltage_source having the same \(R\), waveform, position, and polarisation. A zero-resistance hard source is not equivalent.

This implementation supports 3-D models on the CPU, CUDA, OpenCL, and Metal solvers, including domain-decomposed MPI CPU models, and a nondispersive terminal edge; dispersive materials may exist elsewhere in the model. In an MPI model the sparse terminal state is advanced only by the rank that owns its electric edge, then gathered for final port processing. A terminal may be placed in a CPU or CUDA HSG subgrid through the Python API, where it uses the local spatial and temporal steps. On an accelerator the complete recurrence and local field correction remain device-resident during time stepping. Several independent terminals may be used, but a coupled multiport admittance matrix is not yet supported. See [CHE2007] for the general PLRC lumped-network formulation and Analytical comparisons for a complete loaded-guide comparison.

#voltage_source:

Port and outputs

In a supported 3-D model this source automatically writes /ports/<id> (default ID portN); no extra Rx is needed. A CPU or CUDA HSG port is below /subgrids/<subgrid ID>/ports/<id>. A hard source on a domain-minimum transverse boundary remains an excitation but has no automatic port because its full magnetic-current loop is unavailable. Reduced 2-D sources do not provide this port output.

Vtotal is gap voltage in volts, with sign opposite to the electric-field component along the source edge. Use its stored time axis. Finite-resistance ports use half-step-aligned voltage histories; hard ports retain voltage at time zero and a separately staggered Iloop history on time_current. frequency, S11, Zin and Yin are the corrected terminal spectra, with separate valid_S11, valid_Zin and valid_Yin masks.

The Python source’s result provides time, total_voltage, frequency, s11, zin, yin and their lowercase validity masks after solving. The finite source’s resistance is its reference impedance; a hard source’s independent reference impedance defaults to 50 Ohms. A zero-amplitude finite-resistance source remains a passive load and reports received voltage, but its own source-normalised S11/impedance is not a valid input characterisation. Zero amplitude with zero resistance is a clamp.

python -m gprMax.toolboxes.Plotting.plot_port model.h5 --port feed --validity --save

Replace feed with the actual ID. See Voltage-port equations and timing for correction, timing and support details, Voltage-source S11 and impedance output for the schema, and Finite-resistance voltage-port studies for a complete multiport matrix.

Allows you to introduce a voltage source at an electric field location. It can be a hard source if it’s resistance is zero, i.e. the time variation of the specified electric field component is prescribed, or if it’s resistance is non-zero it behaves as a resistive voltage source. It is useful for exciting antennas when the physical properties of the antenna are included in the model. The syntax of the command is:

#voltage_source: c1 f1 f2 f3 f4 str1 [f5 f6 [f7 | str2 str3 [f7]]]

• c1 is the polarisation of the source and can be x, y, or z.

• f1 f2 f3 are the coordinates (x,y,z) of the source in the model.

• f4 is the internal resistance of the voltage source in Ohms. If f4 is set to zero then the voltage source is a hard source. That means it prescribes the value of the electric field component. If the waveform becomes zero then the source is perfectly reflecting.

• f5 f6 are optional parameters. f5 is a time delay in starting the source. f6 is a time to remove the source. If the time window is longer than the source removal time then the source will stop after the source removal time. If the source removal time is longer than the time window then the source will be active for the entire time window. If f5 f6 are omitted the source will start at the beginning of time window and stop at the end of the time window.

• str2 is an optional port/output identifier. str3 is then required and is either the minimum cells per shortest material wavelength (default 10) or nyquist for every native non-negative FFT bin. Because hash-command arguments are positional, f5 and f6 must be supplied before these port options. If omitted, port1, port2, and so on are assigned.

• f7 is the optional positive wave-reference impedance in Ohms for a hard source only; it is always the final positional value. A hard source defaults to 50 Ohms. A finite-resistance source uses f4 and must not supply f7.

• str1 is the identifier of the waveform that should be used with the source.

For a zero-resistance source, the initial electric field is prescribed at \(t=0\) before receiver/snapshot sampling or the first magnetic update. After each electric update the source prescribes \(E^{n+1}\) using the waveform at \((n+1)\Delta t\). Its start/stop interval is inclusive and tested at this physical electric-field time, including the initial sample. Other source types retain their own staggered update times.

An impulse waveform is nonzero for one sampled time step, but a hard source continues to clamp its edge to zero while it remains active. To prescribe only the initial impulse and then release the edge, use f5=0 and \(0<f6<\Delta t\) (the source’s local time step, if inside a subgrid). f6=\Delta t includes a second, zero-valued prescription. This initialisation occurs once per run, including geometry reuse, on CPU, CUDA, OpenCL, Metal, MPI CPU grids, and supported CPU or CUDA HSG subgrids.

Every supported 3-D voltage port stores its terminal voltage and frequency-domain S11, Zin, and Yin. No separate receiver-port command is required. For example:

#voltage_source: z 0.050 0.050 0.020 50 source_wave 0 10e-9 feed 10
#voltage_source: z 0.060 0.050 0.020 0 source_wave 0 10e-9 ideal_feed nyquist 75

For example, to specify a y directed voltage source with an internal resistance of 50 Ohms, an amplitude of five, and a 1.2 GHz centre frequency Gaussian waveform use: #waveform: gaussian 5 1.2e9 my_gauss_pulse and #voltage_source: y 0.05 0.05 0.05 50 my_gauss_pulse.

Note

• A resistive voltage source adds conductivity to the existing electric edge without replacing its permittivity. The added conductivity is dl / (R * A), where dl is the edge length and A the transverse cell area. For a cubic 1 mm cell and 50 Ohms this is 20 S/m, so a background conductivity of 0.1 S/m becomes 20.1 S/m.

• A finite-resistance source on a dispersive edge includes the complete complex background permittivity in its Yee-gap correction. A hard source on a dispersive edge is not yet supported.

#transmission_line:

Port and outputs

The source automatically writes /tls/tlN in creation order, or /subgrids/<subgrid ID>/tls/tlN for a CPU or CUDA HSG source. Vinc and Vtotal are incident and total line voltages (V); Iinc and Itotal are line currents (A). These are terminal quantities, not a point receiver’s local electric/magnetic fields. The line resistance is the reference impedance.

frequency, S11, Zin, Yin and their validity masks are already calculated. Zin_current is an independent current-wave de-embedding check; do not reconstruct it by simply dividing raw staggered histories. Use python -m gprMax.toolboxes.Plotting.plot_port model.h5 --port tls/tl1 --save to inspect the stored result and its timing. See Transmission-line S11 and impedance output for the line’s spatial/temporal reference plane and masks, and Fixed-topology terminal-source studies for changing generator drives while retaining the passive line termination.

Allows you to introduce a one-dimensional transmission line model [MAL1994] at an electric field location. The transmission line can have a specified resistance greater than zero and less than the impedance of free space (376.73 Ohms). It is useful for exciting antennas when the physical properties of the antenna are included in the model. Transmission lines are supported by the CPU, CUDA, OpenCL, and Metal solvers. The syntax of the command is:

#transmission_line: c1 f1 f2 f3 f4 str1 [f5 f6]

• c1 is the polarisation of the transmission line and can be x, y, or z.

• f1 f2 f3 are the coordinates (x,y,z) of the transmission line in the model.

• f4 is the characteristic resistance of the transmission line source in Ohms. It can be any value greater than zero and less than the impedance of free space (376.73 Ohms).

• f5 f6 are optional parameters. f5 is a time delay in starting the excitation of the transmission line. f6 is a time to remove the excitation of the transmission line. If the time window is longer than the excitation of the transmission line removal time then the excitation of the transmission line will stop after the excitation of the transmission line removal time. If the excitation of the transmission line removal time is longer than the time window then the excitation of the transmission line will be active for the entire time window. If f5 f6 are omitted the excitation of the transmission line will start at the beginning of time window and stop at the end of the time window.

• str1 is the identifier of the waveform that should be used with the source.

Time histories of incident and total voltage and current are saved to the output file. gprMax also calculates S11, input impedance, and input admittance automatically after the simulation. The line resistance is used as the S11 reference impedance; Zin is derived from S11, while an independently de-embedded current result is saved as Zin_current for verification. The frequency axis, validity masks, and lambda/10 mesh limit are stored with the results. No separate receiver-port command is required for a transmission-line source. The complete schema and equations are documented in the Simulation Output section.

For example, to specify a z directed transmission line source with a resistance of 75 Ohms, an amplitude of five, and a 1.2 GHz centre frequency Gaussian waveform use: #waveform: gaussian 5 1.2e9 my_gauss_pulse and #transmission_line: z 0.05 0.05 0.05 75 my_gauss_pulse.

An example antenna model using a transmission line can be found in the examples section.

#magnetic_frill_source:

Port and outputs

The source automatically writes /frills/frillN in creation order, or /subgrids/<subgrid ID>/frills/frillN for a CPU or CUDA HSG source. Vinc is half the generator waveform (V), Vtotal is the solved terminal voltage (V), and Itot is the feed current (A), centred using adjacent magnetic half steps. Zcoax is the wave-reference impedance. This continuous-wire, through-ground-plane feed is not an electric gap source.

frequency, S11, Zin, Yin and their per-quantity validity masks are stored directly. Use python -m gprMax.toolboxes.Plotting.plot_port model.h5 --port frills/frill1 --save; see Magnetic-frill-source S11 and impedance output for the sampling convention and full schema. A zero generator leaves the frill termination present; Fixed-topology terminal-source studies can vary the drive without rebuilding its physical definition.

Allows you to introduce a magnetic-frill (equivalent-feed) source [HYU2009] at an electric field location, for an antenna driven through a PEC ground plane by a coaxial line - the antenna’s own inner conductor passes continuously through the plane, unlike #voltage_source/#transmission_line’s gap-feed model. Complements #transmission_line: it is a different, well-established formulation (not a variant of the two-wire line, and not a general-purpose alternative to #voltage_source), building on the magnetic frill generator of King and Harrison and on Maloney, Smith, & Scott’s FDTD implementation of it. There is no explicit one-dimensional line, no absorbing boundary, and no “magic time step”. The coax’s sub-cell aperture is represented by an equivalent magnetic surface current entering Faraday’s law at the four Yee magnetic-field components immediately surrounding the feed point. The corrected Hyun feed-cell formulation is supported by the CPU, CUDA, OpenCL, and Metal solvers, and by domain-decomposed MPI CPU models. In MPI, the four magnetic feed edges may cross internal rank boundaries: their Ampere-loop terms are combined before the common terminal state is advanced, and each field deposit is applied by its owning rank. The syntax is:

#magnetic_frill_source: c1 f1 f2 f3 f4 str1 [str2 | f5 f6 [str2]]

• c1 is the polarisation of the source and can be x, y, or z - the antenna axis the source drives current along, following the same electrical sign convention already used by gprMax’s Ix/Iy/Iz current output.

• f1 f2 f3 are the coordinates (x,y,z) of the feed point in the model - the position of the antenna’s own inner conductor where it passes through the ground plane.

• f4 is the coax’s characteristic impedance Zcoax (Ohms), and must be finite and greater than zero.

• str1 is the identifier of the waveform that should be used with the source.

• f5 f6 are optional parameters specifying a delay before the incident waveform starts and a time at which that waveform stops. They gate only the incident wave; the coaxial terminal relation remains connected for the rest of the simulation so that late antenna reflections are treated correctly.

• str2 is an optional spectrum limit: either a numeric minimum number of cells per shortest material wavelength or nyquist. It may follow the base command directly, or follow f5 f6 when start and stop times are supplied. The default is 10.

The source must be co-located with a Yee edge of a #thin_wire of the same orientation. gprMax obtains the inner-conductor radius \(a\) from that wire. This is necessary because the radius occurs explicitly in Hyun’s discrete feed-cell equation. An ordinary PEC edge has no unambiguous physical radius and is therefore rejected. In particular, \(a\) cannot be assumed to equal the cell size: that would make the logarithmic correction singular. The normal #thin_wire condition \(a < \min(\Delta u,\Delta v)/2\) applies in the two transverse directions.

The physical outer-conductor/aperture radius \(b\) and the coax filler are represented through Zcoax and are not separate numerical inputs. For a lossless TEM coax with inner-conductor radius \(a\), outer-conductor inner radius \(b\), and filler relative permittivity and permeability \(\varepsilon_{r,c}\) and \(\mu_{r,c}\), respectively,

\[Z_\mathrm{coax} = \frac{1}{2\pi} \sqrt{\frac{\mu_0\mu_{r,c}}{\varepsilon_0\varepsilon_{r,c}}} \ln\!\left(\frac{b}{a}\right) = \frac{\eta_0}{2\pi} \sqrt{\frac{\mu_{r,c}}{\varepsilon_{r,c}}} \ln\!\left(\frac{b}{a}\right).\]

For the usual nonmagnetic filler (\(\mu_{r,c}=1\)), this is commonly written

\[Z_\mathrm{coax} \simeq \frac{60}{\sqrt{\varepsilon_{r,c}}} \ln\!\left(\frac{b}{a}\right)\ \Omega, \qquad b = a\exp\!\left( \frac{Z_\mathrm{coax}\sqrt{\varepsilon_{r,c}}}{60} \right).\]

Here \(\eta_0=\sqrt{\mu_0/\varepsilon_0}\) is the impedance of free space. The permittivity in these equations is that of the coax filler, not necessarily the material above the ground plane surrounding the antenna. Thus a measured or datasheet value of Zcoax may be supplied directly; otherwise it can be calculated from \(a\), \(b\), and the filler. gprMax obtains \(a\) from the attached #thin_wire but does not infer or check \(b\).

Warning

This formulation is only strictly valid while the coax’s aperture radius b is sub-cell: smaller than the discretisation (\(\Delta\)) in the plane perpendicular to the polarisation axis - \(b<\min(\Delta y,\Delta z)\) for x polarisation, \(b<\min(\Delta z,\Delta x)\) for y, or \(b<\min(\Delta x,\Delta y)\) for z. This is a model-validity boundary, not a refinement axis: a physical aperture that is not sub-cell does not converge toward the true coax-fed antenna by refining the mesh, it converges toward a different problem (a continuous PEC plane carrying a fictitious current sheet) - mesh the coax explicitly with #cylinder/#box commands instead in that case. Because gprMax is never given b (only Zcoax), this is not checked automatically - it is the user’s responsibility to confirm it holds for their coax and mesh before trusting the result.

The tangential electric edges of the ground plane at the feed must already be PEC (for example via #plate or #box). The axial edge is supplied by the attached #thin_wire. The sub-cell aperture is invisible to the grid and is represented by the frill term; no gap is cut in the PEC plane.

In Hyun’s notation, the equivalent magnetic current and coax load relation are

\[M_\phi^n = \frac{-2V_\mathrm{inc}^n + Z_0 I_\mathrm{tot}^n} {(\Delta\rho/2)\ln(\Delta\rho/a)}, \qquad V_\mathrm{ab}^n = 2V_\mathrm{inc}^n-Z_0 I_\mathrm{tot}^n.\]

gprMax generalises the cylindrical feed cell to a rectangular Cartesian Yee cell. For each transverse radial direction \(u\), the frill applies

\[F_u = \frac{2}{\ln(\Delta u/a)}\]

to the magnetic-source coefficient. The attached improved thin-wire material already supplies Mäkinen’s orientation-specific \(k_H\) projection, so the complete source coefficient contains \(F_u k_H\). The frill must supply \(F_u\); omitting it would model the wrong inner-conductor radius.

The leapfrog current is evaluated using the time-average approximation recommended in [HYU2009],

\[I_\mathrm{tot}^n = \frac{1}{2} \left(I^{n-1/2}+I^{n+1/2}\right).\]

Because \(I^{n+1/2}\) depends on the frill voltage applied during the same update, gprMax solves this small implicit relation in closed form. If \(G_f\) is the precomputed feed-cell self-admittance and \(I_\mathrm{bulk}^{n+1/2}\) is the current after the ordinary magnetic update but before the new frill deposit, then

\[I^{n+1/2} = \frac{I_\mathrm{bulk}^{n+1/2} + 2G_f V_\mathrm{inc}^n - (G_f Z_0/2) I^{n-1/2}} {1+G_f Z_0/2}.\]

This avoids a forward-time predictor iteration and implements equations (8)–(11) of [HYU2009] directly. The command waveform follows gprMax’s Thevenin-generator convention, so the stored incident wave is one half of the specified waveform amplitude.

Time histories of incident and total voltage (\(V_\mathrm{inc}\), \(V_\mathrm{ab}\)) and total current (\(I_\mathrm{tot}\)) are saved to the output file, along with automatically-calculated S11, input impedance, and input admittance, following the same conventions as #transmission_line. No separate receiver-port command is required. Supply the optional spectrum_limit directly on the magnetic-frill source when needed. The complete schema and equations are documented in the Simulation Output section.

For example, a z-directed, 0.1 mm radius inner conductor driven through a ground plane at z = 0 by a 50 Ohm coax is:

#waveform: ricker 1 10e9 my_pulse
#plate: 0 0 0 0.1 0.1 0 pec
#thin_wire: 0.05 0.05 0 0.05 0.05 0.04 0.0001
#magnetic_frill_source: z 0.05 0.05 0 50 my_pulse

The hash command creates a main-grid source and remains valid in a model that also contains subgrids. To place the frill, its attached thin wire, and its ground plane inside a SubGridHSG, use the Python API and add all four objects (including the waveform) to the same subgrid object.

Note

• This source can be placed at a symmetry-plane corner declared with #symmetry_boundary (for example a monopole fed at the domain origin, to simulate only a quarter of the structure) - but only at the ‘0’-type faces transverse to c1 (y0/z0 for x polarisation, z0/x0 for y, or x0/y0 for z); the corresponding ‘max’-type symmetry corners are not yet supported.

• A feed point placed exactly at a domain-minimum boundary without the matching #symmetry_boundary declared there is rejected outright, since gprMax’s underlying current-loop calculation cannot otherwise distinguish “domain edge” from “symmetry plane”.

• Two frill sources may not share a surrounding H edge. Such adjacent or duplicate feeds form a coupled feed-cell system and cannot be advanced as independent scalar terminal relations.

• MPI symmetry boundaries are supported. A frill and its thin wire may cross internal rank boundaries or use PMC image completion at supported minimum-face symmetry corners.

• This source is a “Path A” (through-ground-plane, continuous-conductor) feed model. It is not intended for a dipole/bow-tie style gap feed (\(E_z \neq 0\) at the feed) - use #voltage_source for that case.

#plane_wave_angles:

Outputs (no circuit port)

A discrete plane wave does not create voltage/current ports or S-parameters. Use receivers inside/outside the TFSF box for total/scattered fields and NTFF outputs for scattering/RCS. The incident-wave normalisation is not antenna-port accepted-power normalisation. See NTFF field-transformation output and Plane-wave and RCS studies for angle/polarisation studies; every observation direction must be requested in the fixed output definition.

Allows you to introduce a discrete plane wave source [TAN2010]. Plane wave sources are a useful tool in multiple different scenarios of electromagnetic simulations, especially when the wave is emitted by a source that is quite far away from the target. The plane wave can originate from any direction and it is assumed that it propagates in a homogeneous background medium. The syntax of the command is:

#plane_wave_angles: f1 f2 f3 f4 f5 f6 f7 f8 f9 str1 [str2 f10 f11]

• f1 f2 f3 are the lower left (x,y,z) coordinates of the total field, scattered field (TFSF) box, and f4 f5 f6 are the upper right (x,y,z) coordinates of the total field, scattered field (TFSF) box.

• f7 is theta which defines the polar propagation angle (degrees) of the incident plane wave.

• f8 is phi which defines the azimuthal propagation angle (degrees) of the incident plane wave.

• f9 is psi which defines the polarisation angle (degrees) of the incident plane wave.

• str1 is the identifier of the waveform that should be used with the source.

• str2 f10 f11 are optional parameters. str2 is a material identifier that is the background material that the plane wave propagates through. The default value is free_space. This material must also be the background material of your full model. f10 is a time delay in starting the excitation of the discrete plane wave. f11 is a time to remove the excitation of the discrete plane wave. If the time window is longer than the excitation of the discrete plane wave removal time then the excitation of the discrete plane wave will stop after the excitation of the discrete plane wave removal time. If the excitation of the discrete plane wave removal time is longer than the time window then the excitation of the discrete plane wave will be active for the entire time window. If f10 f11 are omitted the excitation of the discrete plane wave will start at the beginning of time window and stop at the end of the time window.

For example, to specify a discrete plane wave in a TFSF box (0.010, 0.010, 0.010 to 0.040, 0.040, 0.040) with a polarisation angle \(\psi\) of 90 degrees, azimuthal propagation angle \(\phi\) of 63.4 degrees, polar propagation angle \(\theta\) of 36.7 degrees, and using the waveform defined by the identifier mypulse use: #plane_wave_angles: 0.010 0.010 0.010 0.040 0.040 0.040 36.7 63.4 90.0 mypulse.

Note

• Plane waves support non-dispersive dielectric backgrounds and multi-pole Debye, Lorentz, and Drude media. They do not currently support user-defined waveforms.

• The plane-wave command must be defined on the main grid. Its TFSF box may contain a complete subgrid, but must strictly enclose the subgrid’s HSG outer coupling surface wherever the two regions overlap.

• MPI domain decomposition is supported. Every rank advances an identical, small auxiliary one-dimensional DPW, while TFSF corrections are restricted to the Yee components owned by that rank. This adds no per-timestep plane-wave communication beyond the normal field-halo exchange.

• This plane wave implementation was based on an initial implementation made possible by a Google Summer of Code (GSoC) project and more details can be found in the original pull request.

• Internally, theta and phi are approximated by an integer direction vector (Mx, My, Mz) found to within a maximum acceptable angular difference of 3 arc minutes (0.05 degrees) by default. This tolerance can be relaxed or tightened using the max_angle_diff parameter (in degrees) when using the Python API.

#plane_wave_vector:

Outputs (no circuit port)

A discrete plane wave does not create voltage/current ports or S-parameters. Use receivers inside/outside the TFSF box for total/scattered fields and NTFF outputs for scattering/RCS. The incident-wave normalisation is not antenna-port accepted-power normalisation. See NTFF field-transformation output and Plane-wave and RCS studies for angle/polarisation studies; every observation direction must be requested in the fixed output definition.

Allows you to introduce a discrete plane wave source [TAN2010]. Plane wave sources are a useful tool in multiple different scenarios of electromagnetic simulations, especially when the wave is emitted by a source that is quite far away from the target. The plane wave can originate from any direction and it is assumed that it propagates in a homogeneous background medium. The syntax of the command is:

#plane_wave_vector: f1 f2 f3 f4 f5 f6 i1 i2 i3 f7 str1 [str2 f10 f11]

• f1 f2 f3 are the lower left (x,y,z) coordinates of the total field, scattered field (TFSF) box, and f4 f5 f6 are the upper right (x,y,z) coordinates of the total field, scattered field (TFSF) box.

• i1 i2 i3 are integers that specify the direction of the wave vector (Mx, My, Mz) of the incident plane wave.

• f7 is psi which defines the polarisation angle (degrees) of the incident plane wave.

• str1 is the identifier of the waveform that should be used with the source.

• str2 f10 f11 are optional parameters. str2 is a material identifier that is the background material that the plane wave propagates through. The default value is free_space. This material must also be the background material of your full model. f10 is a time delay in starting the excitation of the discrete plane wave. f11 is a time to remove the excitation of the discrete plane wave. If the time window is longer than the excitation of the discrete plane wave removal time then the excitation of the discrete plane wave will stop after the excitation of the discrete plane wave removal time. If the excitation of the discrete plane wave removal time is longer than the time window then the excitation of the discrete plane wave will be active for the entire time window. If f10 f11 are omitted the excitation of the discrete plane wave will start at the beginning of time window and stop at the end of the time window.

For example, to specify a discrete plane wave in a TFSF box (0.010, 0.010, 0.010 to 0.040, 0.040, 0.040) propagating along the diagonal of your grid using a polarisation angle \(\psi\) of 90 degrees, you can use as a vector (1,1,1) resulting in an azimuthal propagation angle \(\phi\) of 45.0 degrees, polar propagation angle \(\theta\) of approximately 54.736 degrees, and using the waveform defined by the identifier mypulse use: #plane_wave_vector: 0.010 0.010 0.010 0.040 0.040 0.040 1 1 1 90.0 mypulse.

Note

• Plane waves support non-dispersive dielectric backgrounds and multi-pole Debye, Lorentz, and Drude media. They do not currently support user-defined waveforms.

• The plane-wave command must be defined on the main grid. Its TFSF box may contain a complete subgrid, but must strictly enclose the subgrid’s HSG outer coupling surface wherever the two regions overlap.

• MPI domain decomposition is supported. Every rank advances an identical, small auxiliary one-dimensional DPW, while TFSF corrections are restricted to the Yee components owned by that rank. This adds no per-timestep plane-wave communication beyond the normal field-halo exchange.

• This plane wave implementation was based on an initial implementation made possible by a Google Summer of Code (GSoC) project and more details can be found in the original pull request.

#plane_wave_axial:

Outputs (no circuit port)

A discrete plane wave does not create voltage/current ports or S-parameters. Use receivers inside/outside the TFSF box for total/scattered fields and NTFF outputs for scattering/RCS. The incident-wave normalisation is not antenna-port accepted-power normalisation. See NTFF field-transformation output and Plane-wave and RCS studies for angle/polarisation studies; every observation direction must be requested in the fixed output definition.

Allows you to introduce a discrete plane wave source [TAN2010]. Plane wave sources are a useful tool in multiple different scenarios of electromagnetic simulations, especially when the wave is emitted by a source that is quite far away from the target. This command introduces a plane wave that propagates along one of the three grid axes and can be normally incident on multi-layer setups that span the entire model domain perpendicular to the direction of propagation. It takes its media properties from the background materials of the grid at the direction of the axis that it propagates. This allows for half-space simulations but only for normally incident plane waves. The syntax of the command is:

#plane_wave_axial: f1 f2 f3 f4 f5 f6 f7 c1 str1 [f10 f11]

• f1 f2 f3 are the lower left (x,y,z) coordinates of the total field, scattered field (TFSF) box, and f4 f5 f6 are the upper right (x,y,z) coordinates of the total field, scattered field (TFSF) box.

• c1 is a character that specifies the axis along which the incident plane wave propagates and can be x, y, or z. originating at the lower left corner of the TFSF box and propagating in the positive axis direction.

• f7 is psi which defines the polarisation angle (degrees) of the incident plane wave.

• str1 is the identifier of the waveform that should be used with the source.

• f10 f11 are optional parameters. f10 is a time delay in starting the excitation of the discrete plane wave. f11 is a time to remove the excitation of the discrete plane wave. If the time window is longer than the excitation of the discrete plane wave removal time then the excitation of the discrete plane wave will stop after the excitation of the discrete plane wave removal time. If the excitation of the discrete plane wave removal time is longer than the time window then the excitation of the discrete plane wave will be active for the entire time window. If f10 f11 are omitted the excitation of the discrete plane wave will start at the beginning of time window and stop at the end of the time window.

For example, to specify a discrete plane wave in a TFSF box (0.010, 0.010, 0.010 to 0.040, 0.040, 0.040) propagating along the positive x direction using a polarisation angle \(\psi\) of 90 degrees and using the waveform defined by the identifier mypulse use: #plane_wave_axial: 0.010 0.010 0.010 0.040 0.040 0.040 90.0 x mypulse.

Note

• For simulations that do not involve half-space setups it is recommended to use either the #plane_wave_angles or #plane_wave_vector commands instead as the formulations are more efficient and faster if the background medium of propagation for the plane wave is homogeneous.

• Plane waves support non-dispersive dielectric layers and multi-pole Debye, Lorentz, and Drude layers. They do not currently support user-defined waveforms.

• The plane-wave command must be defined on the main grid. Its TFSF box may contain a complete subgrid, but must strictly enclose the subgrid’s HSG outer coupling surface wherever the two regions overlap.

• MPI domain decomposition is supported. The layered one-dimensional material profile is assembled collectively once during model construction using the actual electric, magnetic, and dispersive update-coefficient rows. The compact auxiliary profile is then replicated, and the timestep loop requires no plane-wave-specific MPI communication.

• This plane wave implementation was based on an initial implementation made possible by a Google Summer of Code (GSoC) project and more details can be found in the original pull request.

#eigenmode_band:

Defines the single frequency band shared by every eigenmode port in the model:

#eigenmode_band: str1 f1 f2 i1 [f3 ...]

• str1 is a non-empty band identifier.

• f1 and f2 are the inclusive DFT start and stop frequencies in Hertz.

• i1 is the number of equally spaced output frequencies from f1 to f2, including both endpoints. A one-point band requires f1=f2; a multi-point band requires f2>f1.

• f3 ... are optional extra output frequencies between f1 and f2 inclusive. They are added to the equally spaced frequencies, and the combined list is sorted from low to high. Repeated values appear only once. This is useful when an S-parameter sweep must also contain specific NTFF frequencies exactly.

Exactly one band is required when eigenmode ports are present. Defining the band once guarantees identical DFT bins at every port.

#eigenmode_port:

Port and outputs

EigenmodePort / #eigenmode_port defines both the modal basis and monitor; EigenmodeExcitation selects what is driven. Each numbered port writes /eigenmode_ports/portN (under /subgrids/<subgrid ID> for an owning HSG grid). frequency, incident and outgoing are modal spectra, not lumped voltage/current histories. Follow ModeIndices and the stored power matrices; do not assign an arbitrary 50 Ohm circuit reference.

With one active channel, per-port S arrays describe responses to that channel. Multiple simultaneous drives produce active reflection data, not a complete S matrix; fully passive virtual ports produce raw spectra without S. Use coefficient_valid_S for coefficient inspection and power_wave_valid_S for propagating-power interpretation. Below-cutoff coefficients must not be interpreted as transmitted power.

Use the modal plots and post-processing examples in Eigenmode ports and S-parameter analysis, not the conventional plot_port tool. See Eigenmode-port and S-parameter output for shapes and files and Eigenmode-port studies and array synthesis for complete matrices and array synthesis. Ports must satisfy the owning grid/backend restrictions in the eigenmode guide.

Ports using automatic tracking also write a versioned mode_tracking group. It contains original candidate mappings, overlaps, residuals, requested and adaptive frequencies, and per-anchor confinement/artifact diagnostics. These quality records do not alter the established anchor_mode_valid or power-wave masks; they explain warnings about anchors that remain in use.

Defines an active or passive modal reference plane. The same command is used for every port; one or more #eigenmode_excitation commands select active port/mode channels.

#eigenmode_port: i1 f1 f2 f3 f4 f5 f6 c1 i2[,i3 ...] str1|f7 [f8 ...] [c2] [degenerate=groups] [mode_polarizations=entries]

• i1 is the unique, one-based port number.

• f1 f2 f3 and f4 f5 f6 are opposite port-plane points in metres. In 3D exactly one finite coordinate pair must match, defining the normal. In 2D the normal must be in-plane and the port must span the complete invariant thickness; use inf for the upper invariant coordinate.

• c1 is + or - and points into the device from the port.

• i2[,i3 ...] is a comma-separated, strictly increasing list of positive one-based modes to monitor, for example 1 or 1,2.

• str1 can be auto. One common automatic anchor list covers the requested band and the significant sampled source spectrum, including its transition regions, and is used by every automatic port.

• Alternatively, f7 [f8 ...] are explicit, strictly increasing modal anchor frequencies. Multiple anchors must cover the required range or the model is rejected with suggested anchors. A single explicit anchor is accepted intentionally as a constant modal basis across the complete band.

• c2 optionally controls field plots: y always writes them and n always suppresses them. If omitted, geometry-only runs write the plots and normal runs do not.

• degenerate optionally declares one group (1,2) or disjoint groups separated by semicolons (1,2;3,4). Every member must be in the monitored mode list. Groups must be exactly or numerically unresolved degenerate; resolved splitting is an error.

• mode_polarizations=y (or x/z) sets the first member of every degenerate pair to that global transverse electric direction. The second follows the positive port-normal axis crossed with the first, independently of propagation sign. A real vector such as 1,1,0 is also accepted and normalized. Omitting this option preserves the current automatic assignment (Python default None). Physical selection requires a 3-D cross-section. Automatic tracking resolves subgroups before applying this single direction within each doublet; legacy tracking uses declared pairs. The deprecated shared-pair form y;x remains supported by both trackers and preserves both directions and signs. The per-mode 1:y;2:x and per-group 1,2:y mappings require legacy tracking. These forms emit a warning and will be removed in gprMax 5.0. Prefer one global direction when its inferred partner matches the intended basis; see Eigenmode tracking guide for migration, phase conventions, and nonorthogonal pairs. Zero, nonfinite, complex, or non-transverse directions are rejected. Values are parsed literally without expression evaluation, after the anchor and plotting arguments.

For a z-normal circular TE11 pair, for example:

#eigenmode_port: 1 0 0 0.02 0.05 0.05 0.02 + 1,2 auto y degenerate=1,2 mode_polarizations=y

Mode 1 then means global y electric polarization and mode 2 means global -x, including in the standard modal plots. Excitation syntax is unchanged. See Eigenmode ports and S-parameter analysis for the runnable example.

For example, these two rectangular-waveguide ports share one DFT band and one automatic anchor list:

#eigenmode_band: wg_band 45e9 65e9 81
#eigenmode_port: 1 0.002 0.001 0.001 0.002 0.007 0.005 + 1 auto
#eigenmode_port: 2 0.011 0.001 0.001 0.011 0.007 0.005 - 1 auto

For independently tracked modes, consecutive anchors are checked using normalized modal-field overlap. If explicit multiple anchors show a severe mismatch, such as at a degeneracy or mode crossing, the run stops and recommends one explicit anchor. With auto, a failure confined to a spectral guard outside the requested band trims that tail only for the affected port and mode and uses its nearest retained modal profile for endpoint extrapolation. A failure in the requested band makes that port and mode warn and use its band-centre anchor. Candidate frequencies remain common, but retained masks and fallbacks are per port and per mode.

Declared degenerate groups use subspace tracking instead. Guard trimming and cutoff exclusion apply to whole groups; a failed group never falls back one member at a time. Physical directions remain aligned at every retained anchor. Resolved splitting, rank loss, or failed in-band group matching raises an error; see Degenerate subspaces and physical alignment for the safeguards.

#eigenmode_field_output:

Export the tracked modal E/H bases prepared during model construction. This also works with a geometry-only run; the arrays are modal profiles, not driven time-domain fields:

#eigenmode_field_output: [filename [port ...]]

filename is a basename without directory components and defaults to port_modes. The file is written as filename.modes.h5 in the run output directory. Optional port numbers are whitespace-separated, unique positive integers; omitting them exports all prepared physical ports. For example, #eigenmode_field_output: bank 1 2 exports ports 1 and 2 to bank.modes.h5.

Repeat the command with distinct filenames for separate exports. Missing or unprepared ports raise an error. Export currently requires a serial 3D main grid; MPI and reduced 2D/subgrid requests are unsupported. The output retains the tracked basis, including any assigned degenerate-mode polarizations.

#eigenmode_excitation:

Port and outputs

EigenmodePort / #eigenmode_port defines both the modal basis and monitor; EigenmodeExcitation selects what is driven. Each numbered port writes /eigenmode_ports/portN (under /subgrids/<subgrid ID> for an owning HSG grid). frequency, incident and outgoing are modal spectra, not lumped voltage/current histories. Follow ModeIndices and the stored power matrices; do not assign an arbitrary 50 Ohm circuit reference.

With one active channel, per-port S arrays describe responses to that channel. Multiple simultaneous drives produce active reflection data, not a complete S matrix; fully passive virtual ports produce raw spectra without S. Use coefficient_valid_S for coefficient inspection and power_wave_valid_S for propagating-power interpretation. Below-cutoff coefficients must not be interpreted as transmitted power.

Use the modal plots and post-processing examples in Eigenmode ports and S-parameter analysis, not the conventional plot_port tool. See Eigenmode-port and S-parameter output for shapes and files and Eigenmode-port studies and array synthesis for complete matrices and array synthesis. Ports must satisfy the owning grid/backend restrictions in the eigenmode guide.

Ports using automatic tracking also write a versioned mode_tracking group. It contains original candidate mappings, overlaps, residuals, requested and adaptive frequencies, and per-anchor confinement/artifact diagnostics. These quality records do not alter the established anchor_mode_valid or power-wave masks; they explain warnings about anchors that remain in use.

Selects one active port and mode after the band and ports have been defined. Repeat the command to drive several distinct port/mode channels. Omit it only when every port is a passive virtual guide; such a model writes raw modal spectra but no S matrix:

#eigenmode_excitation: i1 i2 [str1] [f1] [f2] [f3] [c1]

• i1 is an existing port number.

• i2 is one of that port’s monitored mode indices.

• str1 is auto by default. It constructs a finite, real band-pass pulse whose Gaussian-smoothed lower and upper edges adapt to the requested band, simulation time step, and Nyquist limit. The same waveform drives the source and its modal spectrum analysis. The pulse is placed at the earliest causal time that preserves its significant temporal support, leaving the remaining time window for propagation and ring-down.

• A custom waveform identifier may be supplied instead. Its exact sampled spectrum is checked before any modal solve. Significant DC or Nyquist bins are discarded with a warning. More than one percent spectral power outside the requested band remains an error. Use a band-limited waveform, or select auto to synthesize one automatically for a finite frequency band.

• f1 is an optional non-zero amplitude scale. All simultaneous drives must use the same base waveform.

• f2 is an optional constant phase in degrees. It defaults to zero.

• f3 is an optional true time delay \(\tau\) in seconds. It defaults to zero and applies \(\exp(-j2\pi f\tau)\).

• c1 optionally controls the waveform/DFT plot: y always writes it and n always suppresses it. If omitted, geometry-only runs write the plot and normal runs do not. The flag may follow i2 directly when the default waveform and amplitude are used.

A complete excitation for the ports above is:

#eigenmode_excitation: 1 1 auto y

For a simultaneous two-port state, for example:

#eigenmode_excitation: 1 1 auto 1 0 0 n
#eigenmode_excitation: 2 1 auto 0.5 90 0 n

Exactly one active channel produces an S-parameter column. With two or more active channels, gprMax stores the driven incident/outgoing modal waves and drive metadata but does not create S datasets or an S-parameter CSV. Use an eigenmode study for the independently excited complete S matrix.

A single-frequency band cannot use the automatic finite-band pulse. Supply a matching continuous waveform and one explicit modal anchor instead.

At every FDTD time step the source port and all passive ports project the cell-centred transverse fields onto their requested modes. gprMax writes <output>_sparameters.csv and the corresponding /eigenmode_ports HDF5 groups. Source-port rows contain modal S11 results; other ports contain S21 and modal-conversion results. See Eigenmode ports and S-parameter analysis for the phasor convention, FDFD formulation, power normalization, TF/SF injection, broadband synthesis, and output definitions.

Note

• Hash commands define eigenmode ports on the main grid, where the CPU, CUDA, OpenCL, and Metal solvers are supported. Direct eigenmode ports may also be added to an HSG subgrid through the Python API; they then use the local material slice, spatial step, time step, and CPU or CUDA update cycle. The complete port stencil must remain strictly inside the subgrid working region. See Eigenmode ports and S-parameter analysis.

• Domain-decomposed MPI CPU models are supported. Modal material slices are assembled collectively once, TF/SF corrections are restricted to locally owned Yee components, and modal DFT projections are reduced at finalisation.

#virtual_waveguide:

Replaces the continuation behind an eigenmode-port plane with a finite auxiliary FDTD waveguide terminated by a PML:

#virtual_waveguide: i1 [i2] [i3] [i4] [str1]

• i1 is the one-based number of an existing #eigenmode_port.

• i2 is the total auxiliary-guide length in cells (default 30).

• i3 is the remote PML thickness in cells (default 12).

• i4 is the number of clear cells between an active modal source and the PML (default 6).

• str1 optionally selects a reusable PML profile defined by named #pml_formulation and #pml_cfs commands. Use None to retain the global PML configuration while specifying all preceding positional values.

The guide direction and cross-section are inherited from the referenced port. The material and Yee-component cross-section at the port is repeated through the auxiliary grid. The normal H field and tangential E fields are coupled in both directions at the aperture, so reflected modes enter the guide and are absorbed by its remote PML.

If the referenced port is selected by #eigenmode_excitation, the modal source is placed inside the auxiliary guide. Otherwise the guide is passive. All ports may be passive when every #eigenmode_port has a #virtual_waveguide; raw incident and outgoing modal spectra are then saved, but no S-parameters can be normalized without an active port.

The port plane must be internal, locally uniform along the guide axis, and at least two cells wide in each transverse direction. The minimum guide length is i3 + i4 + 3 cells. In 2D, the two-cell width requirement applies only to the physical transverse axis; the full invariant storage dimension is included automatically. Main-grid virtual waveguides support all 2D TE/TM orientations on the CPU, including passive constant/dispersive SIBC and resistance inf PMC walls. The 3D implementation supports non-dispersive guide cross-sections with the CPU, CUDA, OpenCL, and Metal solvers. Through the Python API, a virtual waveguide may instead be attached to an HSG-subgrid port; it then inherits that subgrid’s fine material slice, spatial and temporal steps, and CPU or CUDA update cycle. Main-grid domain-decomposed MPI CPU models are also supported; HSG with MPI is not. The compact MPI auxiliary guide is replicated, while one aperture-sized collective communicates the three required H sheets after each magnetic halo exchange.

Unlike a direct eigenmode source, a virtual-waveguide source lies outside the main FDTD domain. A closed equivalent-current or KSIR NTFF surface may therefore enclose the antenna and its physical feed aperture without omitting a face. The integration surface must not intersect the virtual aperture.

#rx:

Allows you to introduce output points into the model. These are locations where the values of the electric and magnetic field components over the number of iterations of the model will be saved to file. The syntax of the command is:

#rx: f1 f2 f3 [str1 str2]

• f1 f2 f3 are the coordinates (x,y,z) of the receiver in the model.

• str1 is the identifier of the receiver.

• str2 is a list of outputs with this receiver. It can be any selection from Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, or Iz.

Note

• When the optional parameters str1 and str2 are not given all the electric and magnetic field components will be output with the receiver point.

• On CUDA, OpenCL, and Metal, Ix, Iy, and Iz use the same single-cell Ampere loops as the CPU solver. Their device histories are allocated only for explicitly requested current components.

#rx_array:

Provides a simple method of defining multiple output points in the model. The syntax of the command is:

#rx_array: f1 f2 f3 f4 f5 f6 f7 f8 f9

• f1 f2 f3 are the lower left (x,y,z) coordinates of the output line/rectangle/volume, and f4 f5 f6 are the upper right (x,y,z) coordinates of the output line/rectangle/volume.

• f7 f8 f9 are the increments (x,y,z) which define the number of output points in each direction. f7, f8, or f9 can be set to zero to prevent any output points in a particular direction. Otherwise, the minimum value of f7 is \(\Delta x\), the minimum value of f8 is \(\Delta y\), and the minimum value of f9 is \(\Delta z\).

Automatic voltage-source port output

The output summary is included with #voltage_source:. Equations, current/voltage time lattices and support details are now in Voltage-port equations and timing; dataset definitions are in Voltage-source S11 and impedance output.

#sar:

Requests on-the-fly frequency-domain specific absorption rate (SAR) over one or more semantic geometry tags.

Warning

Only electric absorption contributes to SAR; magnetic absorption is not included. Magnetic properties still affect the FDTD fields. Selected magnetic materials outside PML trigger a warning: with nonzero magnetic loss, the result is not total SAR. See the electric-only absorption limitation for details and exceptions.

#sar: f1 f2 i1 str1 f3 f4|nyquist str2 str3 [str4 ...]

where f1 and f2 are the start and stop frequencies in Hz, i1 is the number of linearly spaced frequency points, str1 is the waveform ID, f3 is the positive target peak-phasor amplitude, f4 is the minimum number of cells per shortest material wavelength, str2 is the output ID, and the remaining strings are geometry tags. For one frequency, f1 and f2 must be equal. The target amplitude uses the source’s native units (for example V/m for a discrete plane wave and V for a voltage source). For example:

#material_density: 1040 brain
#sar: 8e8 1.2e9 41 pulse 1 10 brain_sar brain_region

The positional syntax deliberately requires every parameter. nyquist may replace f4 only as an explicit research override; it bypasses the spatial wavelength criterion but not the temporal Nyquist check. The default Python API criterion is lambda/10, and lambda/8 is requested by supplying 8. Every material in the selected final tagged cells must have a positive mass density. Tagged cells inside boundary or internal PML regions are excluded automatically. Three- and two-dimensional CPU, CUDA, OpenCL, and Metal main-grid models, MPI domain-decomposed CPU models, and three-dimensional CPU or CUDA HSG subgrids are supported. In 2-D, only invariant index zero for TM or index one for TE is sampled; tag mass and absorbed power are reported per unit invariant length. A #sar command inside a subgrid block samples the fine grid at its own timestep and writes /subgrids/<subgrid ID>/sar/<output ID>. Its normalising source may belong to the main grid or the subgrid. See Absorbed power density and specific absorption rate for the formulation and datasets, and Two-dimensional SAR in lossy dielectric cylinders for analytical TMz and TEz validation against homogeneous lossy-cylinder series.

For port-power normalisation, replace the numeric target amplitude by incident_power or accepted_power and insert the target power in watts for 3-D, or watts per metre for 2-D, and port ID before the spectrum limit. No waveform ID is required because the named port supplies its own incident and accepted spectra:

#sar: f1 f2 i1 incident_power f3 str1 f4|nyquist str2 str3 [str4 ...]
#sar: 8e8 1.2e9 41 accepted_power 1 feed 10 brain_sar brain_region

The older form containing an unused waveform ID before the power keyword is still accepted for input-file compatibility. A Hertzian current-moment or plane-wave incident-flux normalisation uses:

#sar: f1 f2 i1 str1 current_moment f3 f4|nyquist str2 str3 [str4 ...]
#sar: f1 f2 i1 str1 incident_flux f3 f4|nyquist str2 str3 [str4 ...]

Here current_moment interprets f3 as A m and is restricted to one active 3-D Hertzian dipole. incident_flux interprets it as W/m² and requires one active discrete plane wave using waveform str1.

The original numeric-amplitude form remains unchanged. By default the output contains local cell SAR, absorbed-power density, and per-tag summaries only. To request spatial mass averages, insert the literal spatial_average, the number of target masses, and those masses in kg immediately before the geometry tags:

#sar: f1 f2 i1 str1 f3 f4|nyquist str2 spatial_average i2 f5 [f6 ...] str3 [str4 ...]
#sar: 8e8 1.2e9 41 pulse 1 10 brain_sar spatial_average 2 0.001 0.01 brain_region

Here i2 must be positive and determines exactly how many following values are interpreted as masses. This explicit positional marker keeps arbitrary geometry-tag names distinguishable from numeric masses. Omitting the marker requests no spatial averaging. The same suffix follows str3 (the output ID) in the port-power form. The spatial_average suffix is currently three-dimensional only; it is rejected in TM and TE models.

#radiometry:

Requests density-independent absorbed-power and radiometric weighting over semantic geometry tags. The waveform, current-moment, plane-wave flux, and port-power forms mirror #sar but do not accept a spatial-average suffix:

#radiometry: f1 f2 i1 str1 f3 f4|nyquist str2 str3 [str4 ...]
#radiometry: f1 f2 i1 str1 current_moment f3 f4|nyquist str2 str3 [str4 ...]
#radiometry: f1 f2 i1 str1 incident_flux f3 f4|nyquist str2 str3 [str4 ...]
#radiometry: f1 f2 i1 accepted_power f3 str1 f4|nyquist str2 str3 [str4 ...]

For example, a unit-flux plane-wave absorption cross section for a tagged layer is requested with:

#radiometry: 5e8 2e9 61 incident incident_flux 1 10 layer_absorption subsurface_layer

The output is written to /radiometry/<output ID> and contains local absorbed-power density, normalised absorption density, and per-tag integrals. It does not require #material_density. See Radiometric absorption weighting for the dimensional meaning of each normalisation.

Warning

Absorbed power and radiometric weights include electric absorption only, not magnetic absorption. Magnetic properties still affect the FDTD fields. Selected magnetic materials outside PML trigger a warning. This limitation applies to every normalisation; see the electric-only absorption limitation for details and exceptions.

#network_port:

Port and outputs

Output is explicit: add #network_port: feed or scene.add(gprMax.NetworkPort(terminal_id='feed')) for terminal feed. The network and its excitation alone do not request port results. Output is /ports/feed, or /subgrids/<subgrid ID>/ports/feed for a CPU or CUDA HSG port.

time, Vgenerator, Vtotal and Inetwork describe the generator, gap voltage (V) and network current (A). Terminal-current spectra remove the background Yee-gap admittance and are distinct from the external-network current. frequency, S11, Zin and Yin have independent validity masks. Reference impedance defaults to 50 Ohms; it does not replace the physical network admittance.

Keep the NetworkPort object to read port.result.total_voltage, port.result.s11, port.result.zin and the other result arrays after solving. Omitting NetworkExcitation leaves a passive measurable load: its received voltage/current remain useful, but source-normalised S11 is undefined. Impedance/admittance are available where their own masks permit.

Use python -m gprMax.toolboxes.Plotting.plot_port model.h5 --port feed --save. See Rational-network port output for the schema, Fixed-topology terminal-source studies for fixed-network drive studies and Ports, antenna power and NTFF for NTFF use.

Requests time- and frequency-domain port quantities for an existing #network_terminal. The terminal ID is also the HDF5 port ID. The syntax is

#network_port: str1 [f1 [f2]]

• str1 is the network-terminal ID.

• f1 is the positive real power-wave reference impedance in ohms; it is 50 Ohms by default.

• f2 is the minimum number of cells per shortest material wavelength; it is 10 by default. The literal nyquist retains the full native spectrum with validity masks for research use.

The output contains terminal voltage and current, incident and reflected power-wave spectra, S11, Zin, and Yin. gprMax removes the background Yee-gap capacitance and conductance from the reported terminal current. A driven network obtains S11 and incident power; an unexcited network remains useful as a passive measured port, for which Zin and Yin are still available when numerically defined.

#src_steps: and #rx_steps:

Provides a simple method to allow you to move the location of all simple sources (#src_steps) or all receivers (#rx_steps) between runs of a model. The syntax of the commands is:

#src_steps: f1 f2 f3
#rx_steps: f1 f2 f3

f1 f2 f3 are increments (x,y,z) to move all simple sources (#hertzian_dipole or #magnetic_dipole) or all receivers (created using either #rx or #rx_array commands).

Note

• #src_steps and #rx_steps are not suitable for moving sources which have associated geometry, e.g. antenna models.

#study:

See Reusable studies for task selection, Python equivalents, complete examples and result processing. This entry defines the CSV interface.

Runs a sequence of source/receiver cases while building the model geometry only once. This is the general counterpart to #src_steps and #rx_steps: positions need not follow a regular increment, individual sources can be activated or scaled, and every output file records the exact case that produced it. The syntax is:

#study: gpr file1

The study type may be gpr for irregular source/receiver acquisition, source for fixed-topology stateful terminal sources, port for a finite-resistance voltage-source S-parameter study, eigenmode for a modal-port S-parameter study, or plane_wave for an angular TFSF/RCS study:

#study: source file1
#study: port file1
#study: eigenmode file1
#study: plane_wave file1

file1 is a CSV table. Its path is resolved relative to the main input file. The required columns are case_id and object_id. The optional columns are active, x_m, y_m, z_m, waveform_id, start_s, stop_s, scale, record, port, and mode. Blank cells mean “use the object’s baseline value”. All three position columns must be supplied together and contain absolute coordinates in metres. port and mode are positive integers used only by an eigenmode study.

The CSV B-scan example provides a complete input file, an editable 60-case schedule, and commands to run, merge and plot the traces. Study positions must be finite: for an invariant coordinate in a 2D model, use its resolved field-plane coordinate rather than inf (for example, z_m=0 in the example’s TMz model).

A plane-wave table uses the deterministic object ID plane_wave_1. Its additional optional columns are theta_deg, phi_deg, psi_deg, axis, m_x, m_y, and m_z. The permitted columns follow the plane-wave command used in the model: angles use theta/phi/psi, vector sources use all three integer mapping columns and psi, and axial sources use axis and psi. waveform_id, start_s, stop_s, and a non-zero scale are available for every form. For example:

case_id,object_id,theta_deg,phi_deg,psi_deg,scale
x_incidence,plane_wave_1,90,0,90,1
y_incidence,plane_wave_1,90,90,90,1

The TFSF box and its background material remain fixed. gprMax rebuilds the auxiliary one-dimensional discrete plane wave and resets/recompiles declarative NTFF accumulators for each row, while retaining the main model geometry. NTFF observation directions also remain fixed: request all angles needed across the cases in the input file and select the appropriate direction from each numbered output file.

A source study manages main-grid #transmission_line, #magnetic_frill_source, and #network_excitation commands. Their deterministic IDs are transmission_line_1, magnetic_frill_source_1, and network_excitation_1 (and so on within each family). Positions, impedances, thin-wire/coax geometry, and rational-network definitions are fixed. The CSV may vary active, waveform_id, start_s, stop_s, and scale. For example:

case_id,object_id,active,waveform_id,start_s,stop_s,scale
full_drive,transmission_line_1,true,pulse,0,4e-9,1
half_drive,transmission_line_1,true,pulse,0,4e-9,0.5
passive,transmission_line_1,false,,,,

Any number of listed sources may be active in a case. An omitted or inactive source keeps its physical terminal and acts as a passive termination with zero generator drive. gprMax resets the terminal’s internal field, recurrence, history, and derived port state before every case. Declarative NTFF outputs are also rebuilt with pristine accumulators. source studies are therefore suitable for comparing multiple-feed antenna excitations, but do not calculate a complete S matrix; use port or eigenmode for that purpose. The normal CPU, CUDA, OpenCL, and Metal implementations are available. Stateful sources inside subgrids, MPI, and task farming are not currently supported by this study type.

For example:

case_id,object_id,active,x_m,y_m,z_m,waveform_id,start_s,stop_s,scale,record
trace_1,hertzian_dipole_1,true,0.100,0.050,0.030,,,,1,
trace_1,rx_1,,0.140,0.050,0.030,,,,,true
trace_2,hertzian_dipole_1,true,0.102,0.052,0.030,,,,0.8,
trace_2,rx_1,,0.145,0.052,0.030,,,,,true

Objects receive deterministic IDs from their order of appearance within each object family: hertzian_dipole_1, hertzian_dipole_2, magnetic_dipole_1, voltage_source_1, and rx_1. An explicit #rx identifier is also accepted as an alias. A source listed in a case is active by default; a source omitted from that case, or listed with active=false, is inactive. A receiver omitted from a case remains at its baseline position and is still recorded. record=false is reserved for future selective-output support and is currently rejected rather than silently ignored.

For a port study, every #voltage_source must have finite, non-zero resistance and a unique automatic port ID. The CSV must contain one case for every voltage source and drive exactly one source in each case. Omitted sources retain their fixed source resistance but receive a zero generator waveform, so they behave as passive matched terminations. For example, a two-port schedule is:

case_id,object_id,active,scale
drive_port1,voltage_source_1,true,1
drive_port2,voltage_source_2,true,1

Source position and resistance cannot vary in a port study because they are part of the built electric-edge material. Hard sources are not accepted: zero drive on a hard source enforces zero field and is not a matched termination. Each case output stores its source-plane S-matrix column. After the cases finish, <output>_study.h5 stores the complete source-plane and gap-corrected matrices using S[frequency, output_port, input_port]. The correction removes all numerical gap capacitances/conductances through the full admittance matrix, including the coupled off-diagonal terms.

For an eigenmode study, every case contains the single deterministic object eigenmode_excitation_1 and selects one declared port/mode channel. Every mode on every #eigenmode_port must be selected exactly once:

case_id,object_id,port,mode
p1m1,eigenmode_excitation_1,1,1
p2m1,eigenmode_excitation_1,2,1

All port modal anchors are solved during the first geometry build, including source-grade spectral-guard anchors. Later cases reuse those bases, reset the modal DFT and any #virtual_waveguide fields/PML history, and switch the active channel without another FDFD solve. The aggregate file stores S[frequency, output_channel, input_channel] together with channel_ports and channel_modes. It retains the full measured incident and outgoing matrices from all cases and obtains \(S\) from the conditioned system \(B=SA\); it does not assume that nominally passive ports have exactly zero incident wave. The same de-embedding is applied to any embedded far-field basis selected by an array codebook.

The number of CSV cases determines the number of model runs, so -n is not required. -i N restarts at case N and retains absolute output numbering. The original geometry, materials, PMLs, and grid allocation are reused, but field arrays and receiver histories are reset before every case. Each output contains a /study group with the case ID, resolved parameters, and a copy of the CSV source.

Note

The five study families have different source and subgrid restrictions; see Execution and reuse. source studies support transmission lines, magnetic frills and rational-network excitations; plane_wave studies support a main-grid discrete-plane-wave template. MPI domain decomposition and task farming are not enabled for any study family.

#array_codebook:

Associates a versioned JSON array-state codebook with an eigenmode study:

#study: eigenmode modal_cases.csv
#array_codebook: array_states.json

The command has one positional argument, resolved relative to the main input file. It is valid only with #study: eigenmode. The study CSV still defines the independent one-active-channel solves needed to obtain the complete S matrix. The JSON file defines any coherent states to evaluate after that matrix has been assembled and, optionally, the far-field requests whose linear complex fields should be retained as an embedded basis:

{
  "schema": "gprMax-array-codebook-v1",
  "embedded_far_fields": [
    {"transform_id": "antenna_band", "output_id": "pattern"}
  ],
  "states": [
    {
      "id": "broadside",
      "drives": [
        {"port": 1, "mode": 1, "power_w": 1.0},
        {"port": 2, "mode": 1, "power_w": 1.0,
         "phase_deg": 0.0, "delay_s": 0.0}
      ]
    }
  ]
}

power_w is incident modal power, so the complex power-wave magnitude is its square root. phase_deg is a frequency-independent phase shifter and delay_s is a true time delay. Both default to zero and power_w defaults to one. Omit a channel from a state to leave it passive. Every selected far field must identify one existing frequency-domain KSIR or equivalent-current far-field request. Retaining fields is opt-in because the requested directions, a full-sphere radiation quadrature, and every modal channel are stored.

#snapshot:

Allows you to obtain information about the electromagnetic fields within a volume of the model at a given time instant. By default, the files use the open source Visualization ToolKit (VTK) VTKHDF format and require ParaView 5.12 or newer for viewing. The syntax of this command is:

#snapshot: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 file1

or

#snapshot: f1 f2 f3 f4 f5 f6 f7 f8 f9 i1 file1

• f1 f2 f3 are the lower left (x,y,z) coordinates of the volume of the snapshot in metres.

• f4 f5 f6 are the upper right (x,y,z) coordinates of the volume of the snapshot in metres.

• f7 f8 f9 are the spatial discretisation of the snapshot in metres.

• f10 or i1 are the time in seconds (float) or the zero-based iteration number (integer) which denote the electric-field time level at which the snapshot will be taken.

• file1 is the name of the file where the snapshot will be stored. Snapshot files are automatically stored in a directory with the name of the input file appended with ‘_snaps’. For multiple model runs each model run will have its own directory, i.e. ‘_snaps1’, ‘snaps2’ etc…

For example to save a snapshot of the electromagnetic fields in the model at a simulated time of 3 nanoseconds use: #snapshot: 0 0 0 1 1 1 0.1 0.1 0.1 3e-9 snap1

Snapshot timing follows the Yee leapfrog convention. A requested real time is mapped to the nearest full time step \(n\Delta t\) (a halfway value is mapped to the earlier step), and an integer selects that same zero-based time-step index \(n\) directly. The electric fields are therefore at \(n\Delta t\); the magnetic fields are the native staggered values at \((n-1/2)\Delta t\). No temporal averaging is applied. Valid indices run from zero to one less than the model’s number of iterations. The existing spatial collocation of the Yee components is independent of this timing convention.

Spatial samples lie at the centres of the regular output cells described by the file’s origin and spacing. For spacing larger than one native cell, each component is linearly interpolated from its surrounding native Yee samples; this is not a volume average. Stride-one arithmetic is unchanged. A non-dividing interior extent retains ceil((p2-p1)/dl) output cells, including its final regular cell. If that cell’s centre or interpolation stencil falls outside the physical Yee grid, model construction raises an error: p2 and dl are no longer silently clipped or enlarged. Reduce the requested extent or spacing in that case. In MPI, all ranks participate in sparse native-sample exchange only at the requested snapshot iterations, including ranks that own no output cells.

In a reduced 2-D model, the invariant-axis range is collapsed to the genuine field plane: index zero for TM and index one for TE. The two-cell TE thickness is required by Yee staggering and is not exported as two physical snapshot planes.

Tip

A series of snapshots can be more easily defined using a loop and our Python API, see Surface-impedance reproducibility metadata.

The Python API can also add a snapshot to an HSG subgrid. Such a snapshot uses the subgrid’s finer spatial discretisation and time step, while its file origin remains in the global model coordinate system. The hash command above always defines a main-grid snapshot.

Near-to-far-field transformation commands

The NTFF commands separate the integration surface from the formulation and its output points. A closed surface can therefore be reused by KSIR and equivalent-current transforms, and by many output directions, without repeating the six surface coordinates. An open surface is specific to the frequency-domain Huygens/equivalent-current transform, which can use any user-selected nonempty subset of the six faces. All optional parameters use the traditional positional gprMax syntax; name=value tokens are not used. The derivations, field normalisation, Yee placement, and comparison of the three implemented transforms are given in Near-to-far-field transformations.

The following conventions apply to every NTFF command:

• coordinates and radii are in metres, frequencies are in Hz, and angles are in degrees;

• requested frequencies must not exceed the temporal Nyquist limit \(1/(2\Delta t)\) for the model time step;

• \(\theta\) is the polar angle measured from +z and \(\phi\) is the azimuth measured from +x towards +y;

• a spherical coordinate is relative to the centre of its integration surface. The Python API can instead give a custom surface origin;

• Cartesian outputs are Ex Ey Ez Hx Hy Hz. Spherical outputs are Er Etheta Ephi Hr Htheta Hphi. Far-field derived outputs are described under #ksir_far_field below;

• every exact time- or frequency-domain point must be strictly outside the completed integration surface. A point may be outside the FDTD model domain;

• the sampled surface and exterior must be one homogeneous, lossless, non-dispersive material for KSIR, conventional #ntff_frequency, and the transient transform. A #ntff_layered_frequency surface may instead cross the declared planar interfaces; its two semi-infinite observation materials must be lossless;

• surface samples must remain outside the PML and clear of the TFSF correction stencil. A closed surface must enclose the radiating source or the complete TFSF box and scatterer. An open Huygens surface instead permits an impressed source outside only through one of its omitted faces;

• the implementation requires a three-dimensional model and does not support geometry-fixed reuse. NTFF commands are main-grid objects, but their closed surface may contain complete subgrids. A surface that overlaps a subgrid must strictly enclose its HSG outer coupling surface; it cannot touch or cut the coupling region. Both time- and frequency-domain KSIR collection and frequency-domain equivalent-current collection are available with the serial CPU, CUDA, OpenCL, Metal, and MPI domain-decomposition solvers. Equivalent-current transient far fields are available with the same solvers. MPI uses the CPU field-update backend;

• CPU collection uses the Cython/OpenMP implementation. Accelerator surface state and time-domain output storage remain on the device during FDTD iterations and are transferred to the host once, after the solve. CUDA and OpenCL are hardware-qualified on the development server. Metal has source-generation and dispatch coverage, plus a homogeneous, frequency-domain KSIR integration test for execution on Apple hardware. This test does not qualify symmetry completion, time-domain transforms, equivalent-current transforms, or layered media on Metal. Hosted CI does not run Apple GPU tests. With MPI, every surface patch is sampled by the rank which owns its inside Yee sample; the neighbouring outside sample is read from the normal one-cell halo. There is therefore no additional NTFF communication inside the FDTD iteration. Compact time histories are reduced and frequency-domain surface phasors are assembled on the coordinator after time stepping, before the normal HDF5 output is written.

Directivity, gain, efficiency, and port normalisation are post-processing operations after the FDTD solve. Angular KSIR and equivalent-current evaluation uses Cython/OpenMP kernels on the host for CPU and accelerator simulations alike; it does not add a new per-iteration GPU operation or transfer fields back to the CPU during time stepping.

A minimal dipole workflow can reuse one surface for an exact KSIR time-domain point and an equivalent-current frequency-domain radiation pattern:

#domain: 0.1 0.1 0.1
#dx_dy_dz: 0.002 0.002 0.002
#time_window: 10e-9

#waveform: gaussiandot 1 1e9 pulse
#hertzian_dipole: z 0.05 0.05 0.05 pulse

#ntff_surface: 0.03 0.03 0.03 0.07 0.07 0.07 radiation_surface
#ksir_time_rx: 0.12 0.05 0.05 radiation_surface transient Ez first_arrival
#ntff_frequency: radiation_surface antenna_band 0.8e9 1.0e9 1.2e9 hann
#ntff_far_field_array: 0 180 5 0 360 5 antenna_band pattern Etheta Ephi radiation_intensity

The observation point may lie outside the FDTD domain because KSIR evaluates the homogeneous exterior analytically. The integration surface itself must be inside the non-PML FDTD region and enclose the source.

KSIR independently reconstructs each Cartesian field component and is the only formulation that provides finite-distance fields. The conventional equivalent-current formulation combines tangential electric and magnetic fields and provides far-zone fields only. Its commands use the ntff prefix. Both frequency-domain formulations support the same far-field, antenna, and RCS outputs, so they can be requested together for comparison.

#ntff_surface:

Defines a reusable Yee-aligned cuboidal integration surface:

#ntff_surface: x1 y1 z1 x2 y2 z2 surface_id [omit_face1 ... omit_face5]

x1 y1 z1 and x2 y2 z2 are the lower and upper logical corners. surface_id must be unique and must not contain /. Zero to five trailing face names may be supplied from x0, xmax, y0, ymax, z0, and zmax; duplicates are not allowed. Listed faces are omitted from the frequency-domain Huygens/equivalent-current integral. At least one face must remain active. An impressed source outside the Huygens volume must enter through an omitted face. A feed crossing an opening should remain uniform to the corresponding PML, with its impressed source plane outside the volume. Every sampled face must lie in the homogeneous exterior. This open-surface form is rejected by every KSIR/Ramahi command and by the transient equivalent-current transform.

An arbitrary open surface is not a mathematically exact closure. It is useful when the omitted contribution is intentionally excluded, but its far field can depend on the selected faces. Check convergence by moving and enlarging the sampled faces. The edge rows shared with an omitted face are excluded from the homogeneous-material check, allowing an active side face to terminate on a PEC backplane.

Warning

The NTFF integration surface is not closed. Equivalent-current NTFF normally assumes a closed Huygens surface. This option is intended for configurations where the omitted face is associated with an eigenmode port or other modelling scenarios that require an open Huygens surface. Results may be incomplete or inaccurate if the omitted field contribution is not represented correctly or is significant for your calculations.

Without omitted face names, the surface is physically closed unless one or more faces coincide exactly with a declared #symmetry_boundary. Coincident PEC/PMC faces are then omitted from direct sampling and completed automatically using the exact image parity, reflected normal, edge quadrature, and propagation distance for every component.

For example:

#ntff_surface: 0.034 0.034 0.034 0.066 0.066 0.066 radiation_surface

An antenna fed from the negative x direction can omit its feed face:

#ntff_surface: 0.034 0.020 0.020 0.080 0.080 0.080 radiation_surface x0

A leaky-wave antenna can omit both waveguide end faces so that a passive eigenmode receiver beyond xmax still measures S21:

#ntff_surface: 0.034 0.020 0.020 0.080 0.080 0.080 radiation_surface x0 xmax

A surface whose lower side terminates on a PEC backplane can omit z0:

#ntff_surface: 0.020 0.020 0.030 0.080 0.080 0.080 radiation_surface z0

#ksir_frequency:

Declares a streaming frequency transform for a previously defined surface:

#ksir_frequency: surface_id transform_id f1 [f2 ...] [window]

transform_id must be globally unique. At least one non-negative frequency is required. The optional final window is rectangular (the default) or hann. Frequencies above the temporal Nyquist limit are rejected rather than being silently aliased. Frequencies are accumulated directly during time stepping; field histories are not retained. Surface phasors are saved under the transform’s HDF5 group so that it remains useful even when it has no receiver or far-field command.

The engineering convention is used throughout: phasors have time dependence exp(+j*omega*t), the forward transform kernel is exp(-j*omega*t), and the outgoing Green function contains exp(-j*k*R).

A frequency transform still requires a sufficiently long FDTD time window. The hann window reduces leakage from a truncated non-zero tail, but it is not a replacement for allowing the physical surface fields to decay. A rectangular window retains the unmodified engineering phasor and is required by gain normalisation, but is more sensitive to end-of-record truncation.

#ksir_frequency: radiation_surface antenna_band 0.8e9 1.0e9 1.2e9 hann

#ntff_frequency:

Declares a streaming frequency transform using the conventional equivalent-current method of Luebbers et al. [LUE1991]:

#ntff_frequency: surface_id transform_id f1 [f2 ...] [window]

The frequency, window, engineering convention, Nyquist check, and globally unique transform-ID rules are identical to #ksir_frequency.

When #ntff_antenna_ports associates an eigenmode port with this transform, every transform frequency must be present in the final #eigenmode_band DFT grid. The transform frequencies may be a strict subset. If a requested NTFF frequency is not one of the uniformly spaced modal bins, append it after the point count; gprMax sorts and deduplicates the combined grid. For example:

#eigenmode_band: antenna_band 8e9 12e9 101 8.5e9 9.5e9 10.5e9 11.5e9
#ntff_frequency: radiation_surface antenna_pattern 8e9 8.5e9 9e9 9.5e9 10e9 10.5e9 11e9 11.5e9 12e9 rectangular

The extra modal bins affect the S-parameter DFT grid; they do not replace the auto or explicit modal anchor policy on #eigenmode_port.

The tangential Yee fields are arithmetically collocated on common cell-face centres and form

\[\mathbf J_s=\hat{\mathbf n}\times\mathbf H, \qquad \mathbf M_s=-\hat{\mathbf n}\times\mathbf E.\]

If

\[\mathbf N=\oint_S\mathbf J_s e^{+jk\hat{\mathbf r}\cdot(\mathbf r'-\mathbf r_0)}\,\mathrm dS', \qquad \mathbf L=\oint_S\mathbf M_s e^{+jk\hat{\mathbf r}\cdot(\mathbf r'-\mathbf r_0)}\,\mathrm dS',\]

the stored range-normalized electric far field is

\[\mathbf F_E =-\frac{jk}{4\pi} \left[\eta\left(\mathbf N- \hat{\mathbf r}(\hat{\mathbf r}\cdot\mathbf N)\right) -\hat{\mathbf r}\times\mathbf L\right].\]

Here \(\mathbf r_0\) is the surface phase origin. This transform cannot be used by the finite-distance #ksir_frequency_rx commands. It is consumed by #ntff_far_field, #ntff_far_field_array, and optionally #ntff_antenna_ports. It accepts either all six physical integration faces or any nonempty subset selected by trailing omitted face names on #ntff_surface. Symmetry-completed equivalent-current surfaces are not yet enabled.

#ntff_frequency: radiation_surface current_band 0.8e9 1.0e9 1.2e9 hann

#ntff_layered_background:

Declares the material order and physical interfaces of a planar background:

#ntff_layered_background: background_id axis material0 [interface0 material1 ...]

axis is x, y, or z. Materials run from the positive-axis semi-infinite exterior to the negative-axis semi-infinite exterior. Interface coordinates are in metres and must therefore be strictly descending. There is exactly one fewer interface than materials. For example, air above a finite dry-soil layer and a wet-soil lower half-space is

#ntff_layered_background: ground z free_space 0.0 dry_soil -0.1 wet_soil

This command describes the analytical propagation background; it does not build FDTD geometry. The corresponding boxes/layers must be built separately at the same coordinates. Repeating a material on both sides of an interface is permitted and is useful for homogeneous-reduction checks.

The built-in pec identifier may be used once as the first or last material. It then represents a perfect short-circuit termination at its adjacent interface, rather than a semi-infinite propagation material. For a 2 mm substrate below free space and backed by a PEC plane, for example,

#ntff_layered_background: grounded z free_space 0.0 substrate -0.002 pec

The FDTD substrate and PEC geometry must still be built separately. Far-field directions through the PEC are invalid. Directivity and radiated power use the open hemisphere with conventional full-sphere directivity normalisation; exterior_power, exterior_efficiency, and exterior_maximum are not defined for a terminated stack.

#ntff_layered_frequency:

Declares a streaming equivalent-current transform which uses a previously defined planar background:

#ntff_layered_frequency: surface_id transform_id background_id f1 [f2 ...] [window]

The frequency and window rules are identical to #ntff_frequency. The surface may cross one or more declared layer interfaces. Internal materials may be electrically conductive or Debye, Lorentz, Drude, or inclusive multipole dispersive materials; their complex properties are evaluated at each requested frequency. The positive- and negative-axis exterior materials must be lossless. A terminal PEC is supported as described above; internal PEC layers and PMC layers are not supported.

The transform is consumed by the existing #ntff_far_field, #ntff_far_field_array, and #ntff_antenna_ports commands. It returns far-zone quantities only. For directions in the positive-axis hemisphere, range normalisation and radiation intensity use the positive exterior wavenumber and impedance; the negative exterior is used in the opposite hemisphere. Exact grazing directions (zero direction cosine along axis) are rejected. Avoid requesting such a point explicitly; internal directivity quadrature omits it automatically.

Far-field and antenna metrics support different lossless exterior materials. Coherent array-codebook synthesis is currently limited to stacks whose two exterior impedances are equal and frequency independent.

Three grouped outputs are specific to a planar-layered transform:

• exterior_power stores radiated power, radiated fraction, and the positive-to-negative-axis integrated power ratio;

• exterior_maximum stores each exterior’s peak radiation intensity, conventional full-sphere-normalised directivity, and peak direction;

• exterior_efficiency stores accepted- and incident-power-normalised coupling efficiencies and maximum gain for each exterior. It requires #ntff_antenna_ports and a rectangular transform window.

The first two outputs do not require an antenna-port association. The region names follow the declared stack axis rather than assuming that +z is air or -z is ground.

#ntff_layered_frequency: radiation_surface ground_band ground 0.8e9 1.0e9 1.2e9 rectangular
#ntff_antenna_ports: ground_band feed
#ntff_far_field: 30 0 ground_band upper_cut Etheta Ephi directivity_dbi exterior_power exterior_efficiency exterior_maximum
#ntff_far_field: 150 0 ground_band lower_cut Etheta Ephi directivity_dbi

#ntff_layered_time:

Declares a direct time-domain equivalent-current transform for a previously defined lossless planar background [CAP2007]:

#ntff_layered_time: surface_id transform_id background_id [impulse_tolerance [max_impulses]]

The closed surface may cross declared interfaces. Every stack material must be positive, lossless, and nondispersive. Requested directions must be non-grazing and propagating in every layer; use #ntff_layered_frequency for conductive, dispersive, evanescent, or total- internal-reflection cases. impulse_tolerance is the relative per-path cutoff for the multiple-reflection impulse train and defaults to 1e-10. max_impulses is a safety limit on processed, time-coalesced path events and defaults to 100000. Reducing the tolerance or observing close to a critical angle may require a larger value and a longer FDTD time window.

For a PEC-terminated stack, the surface may either remain closed above the termination or place one boundary exactly on the PEC coordinate and omit only that face with #ntff_surface. For example, a negative-z termination uses the trailing z0 option. The grounded Green function supplies the image cancellation at that plane. No other face may be omitted from a direct time-domain transform.

The transform is available with the CPU, MPI CPU, CUDA, OpenCL, and Metal solvers. CPU and MPI surface gathering and sparse impulse deposition use Cython/OpenMP; accelerator backends keep the sampled currents, TE/TM impulse metadata, and accumulated far-field traces on the device until finalisation.

#ntff_layered_time: radiation_surface ground_transient ground 1e-10 100000

#ksir_antenna_ports:

Associates a complete set of physical antenna ports with a frequency transform so that accepted power, gain, realized gain, and efficiency can be calculated:

#ksir_antenna_ports: transform_id port_id1 [port_id2 ...]

For a voltage source, port_id is the source’s optional ID, or its automatically assigned portN ID. Transmission-line and magnetic-frill sources provide automatic port IDs tl1, tl2, … and frill1, frill2, … respectively, in source creation order. Rational-network ports use their terminal ID and require an explicit #network_port. See Ports, antenna power and NTFF for passive-load power accounting. The association is not required for electric or magnetic far fields, radiation intensity, RCS, or directivity. It is required when a far-field command asks for gain or efficiency.

Every Ramahi/KSIR command requires a closed six-face #ntff_surface; this rule is independent of whether an eigenmode source is present. Traditionally a real guide extends through the domain PML and the crossed Huygens face is omitted, so that open surface cannot be reused for KSIR. With #virtual_waveguide the matched continuation and source are outside the main domain, permitting a closed KSIR or equivalent-current surface around the physical antenna.

A port on a subgrid is named as subgrid_id/port_id. For example, fine_grid/feed identifies a voltage source called feed on subgrid fine_grid; automatic source ports use forms such as fine_grid/tl1 and fine_grid/frill1. Main-grid IDs remain unqualified. Each subgrid port is post-processed using its owning grid’s finer spatial and temporal steps.

The listed set must include every physical voltage, transmission-line, magnetic-frill and explicitly monitored rational-network port in the model. Every excited network terminal must have a #network_port. This requirement makes the net accepted power unambiguous in coupled multiport antennas. A source whose waveform amplitude is zero is still a terminated physical port: list it normally. It has zero incident power, but coupling from driven ports can make its accepted power negative because it delivers coupled energy into its termination. gprMax sums all signed port powers when normalising antenna gain.

Gain normalisation currently requires the transform to use the rectangular window. Active Hertzian electric or magnetic dipoles and plane-wave sources cannot be mixed with a port-normalised antenna result, because their input power is not represented by this port set. The normal per-port wavelength-sampling limit also applies to gain validity. For a voltage-source or magnetic-frill port, an explicit nyquist research override on the source retains the full temporal band, including spatially under-resolved values, as it does for S11 and impedance.

For example, a two-element array with one driven and one terminated element uses:

#waveform: ricker 1 1e9 driven
#waveform: ricker 0 1e9 terminated
#voltage_source: z 0.045 0.050 0.050 50 driven 0 10e-9 element1 10
#voltage_source: z 0.055 0.050 0.050 50 terminated 0 10e-9 element2 10
#ksir_frequency: radiation_surface antenna_band 0.8e9 1.0e9 1.2e9 rectangular
#ksir_antenna_ports: antenna_band element1 element2
#ksir_far_field_array: 0 180 2 0 360 2 antenna_band pattern gain realized_gain

#ntff_antenna_ports:

Associates antenna ports with an equivalent-current frequency transform:

#ntff_antenna_ports: transform_id port_id1 [port_id2 ...]

Its port-set, rectangular-window, multiport-power, subgrid, and validity rules are the same as for #ksir_antenna_ports. The separate command name prevents accidentally associating a port set with a transform from the other formulation.

This is the antenna-port association to use with an eigenmode excitation. Every #eigenmode_port uses portN, where N is its explicit port number. The listed set must include every physical conventional and modal port in the model. A passive eigenmode port has zero generator incident power and contributes signed net modal power to the accepted-power balance.

For every associated eigenmode port, every transform frequency must be present in that port’s direct-DFT bins. The eigenmode grid may contain additional bins, so S-parameters may use a denser grid than the NTFF. Degenerate modes and mode crossings should use a single-frequency modal solve rather than broadband profile interpolation.

#ksir_time_rx: and #ksir_time_rx_spherical:

Request exact physical time-domain fields at one Cartesian or spherical point:

#ksir_time_rx: x y z surface_id [rx_id [output1 output2 ... [time_origin]]]
#ksir_time_rx_spherical: r theta phi surface_id [rx_id [output1 output2 ... [time_origin]]]

The Cartesian command defaults to all six Cartesian components. The spherical command defaults to all six spherical components. rx_id is optional and is generated as rx1, rx2, and so on when omitted. time_origin is the final token and is either:

• simulation (default): retain time from the start of the FDTD run; or

• first_arrival: omit the guaranteed retarded-time zero prefix separately for each point while recording its absolute physical origin.

Optional parameters are positional. Therefore an ID must be supplied before component names, and both the ID and any desired components must precede time_origin. For example:

#ksir_time_rx: 0.074 0.05 0.051 radiation_surface outside Ez Hy first_arrival
#ksir_time_rx_spherical: 0.25 90 0 radiation_surface principal Etheta Ephi simulation

The spherical radius is explicit because these commands return the actual finite-distance field, including all 1/R and 1/R^2 terms. It is not a normalization constant.

The output file stores both fully_supported_lengths and valid_lengths. Use the former by default: it stops before retarded propagation causes any surface patch to run beyond the available FDTD history. The latter exposes the longer partial tail for research use. gprMax also records a terminal decay ratio and warns when the field has not decayed adequately within the fully supported interval; in that case the simulation time window should be increased.

#ksir_time_rx_array:

Defines a Cartesian line, plane, or volume of exact time-domain points:

#ksir_time_rx_array: x1 y1 z1 x2 y2 z2 dx dy dz surface_id [rx_id [output1 output2 ... [time_origin]]]

The bounds are inclusive and must contain an integer number of increments. An increment can be zero only on an axis whose lower and upper coordinates are equal. All points share one output ID and are stored as the first dimension of each field dataset.

#ntff_time_far_field: and #ntff_time_far_field_array:

Request transient far-zone fields using homogeneous time-domain equivalent-current NTFF:

#ntff_time_far_field: theta phi surface_id [output_id [output1 output2 ...]]
#ntff_time_far_field_array: theta1 theta2 dtheta phi1 phi2 dphi surface_id [output_id [output1 output2 ...]]

The default outputs are Etheta Ephi. Any Cartesian or spherical electric or magnetic component may be requested. These are range-normalized far-zone waveforms, not finite-distance receivers, so the commands deliberately have no radius. If \(\tau=t-r/c_b\), the electric result is

\[\mathbf F_E(\hat{\mathbf r},\tau) =-\frac{1}{4\pi c_b}\oint_S \left[\eta\,\dot{\mathbf J}_{s,t} -\hat{\mathbf r}\times\dot{\mathbf M}_s\right] \left(\tau+\frac{\hat{\mathbf r}\cdot (\mathbf r'-\mathbf r_0)}{c_b}\right)\,\mathrm dS'.\]

Here \(\mathbf J_{s,t}\) is the component of \(\mathbf J_s=\hat{\mathbf n}\times\mathbf H\) transverse to the observation direction and \(\mathbf M_s=-\hat{\mathbf n}\times\mathbf E\). The electric and magnetic current differences retain their natural half-step offset: the derivative of \(\mathbf M_s^n\) is placed at \((n-1/2)\Delta t\), and the derivative of \(\mathbf J_s^{n+1/2}\) at \(n\Delta t\). No extra interpolation is used to force them onto one FDTD time level. This is the modification to the original Luebbers method [LUE1991] introduced by Giannopoulos et al. [GIAFF1997]. Linear interpolation is used only for fractional propagation delays.

The stored reduced-time axis excludes both the range-dependent leading-zero delay and all final bins that are not supported by every integration patch. This prevents an incomplete retarded-time tail being mistaken for a physical late-time response. Increase #time_window if the stored terminal-decay test fails.

The current implementation supports the CPU, CUDA, OpenCL, and Metal solvers, requires a homogeneous lossless background and all six physical faces, and supports 3-D serial models. KSIR remains available for finite-distance time-domain points and for symmetry-completed surfaces.

#ntff_time_far_field_array: 0 180 2 0 360 2 radiation_surface transient Etheta Ephi

#ntff_layered_time_far_field: and #ntff_layered_time_far_field_array:

Request range-normalized transient far fields from a #ntff_layered_time transform:

#ntff_layered_time_far_field: theta phi transform_id [output_id [output1 output2 ...]]
#ntff_layered_time_far_field_array: theta1 theta2 dtheta phi1 phi2 dphi transform_id [output_id [output1 output2 ...]]

The defaults and available Cartesian/spherical field components are the same as for #ntff_time_far_field. The reduced time is \(t-r/c_o\), where \(c_o\) is the wave speed of the positive- or negative-axis observation half-space selected by each direction. The corresponding exterior impedance is used to derive the magnetic field. Only the interval supported by the complete FDTD history and every retained multiple-reflection impulse is written.

#ntff_layered_time_far_field_array: 5 175 5 0 355 5 ground_transient transient_pattern Etheta Ephi

#ksir_frequency_rx: and #ksir_frequency_rx_spherical:

Request the exact finite-distance physical phasor at one point using a previously declared transform:

#ksir_frequency_rx: x y z transform_id [rx_id [output1 output2 ...]]
#ksir_frequency_rx_spherical: r theta phi transform_id [rx_id [output1 output2 ...]]

The Cartesian and spherical component defaults match their time-domain counterparts. These results retain the full outgoing Green function and are not range normalized. Their arrays have shape (nfrequencies, npoints); npoints is one for these two commands.

#ksir_frequency_rx: 0.074 0.05 0.051 antenna_band near_phasor Ez
#ksir_frequency_rx_spherical: 0.25 90 0 antenna_band spherical_phasor Etheta Ephi

#ksir_frequency_rx_array:

Defines a Cartesian line, plane, or volume of exact frequency-domain points:

#ksir_frequency_rx_array: x1 y1 z1 x2 y2 z2 dx dy dz transform_id [rx_id [output1 output2 ...]]

The inclusive bounds and zero-increment rule are the same as for #ksir_time_rx_array.

#ksir_far_field:

Requests a range-normalized far field in one spherical direction:

#ksir_far_field: theta phi transform_id [output_id [output1 output2 ...]]

The default outputs are Etheta Ephi. Cartesian and spherical electric or magnetic components may be requested. The derived outputs are radiation_intensity, directivity, directivity_dbi, gain, gain_dbi, realized_gain, realized_gain_dbi, radiation_efficiency, total_efficiency, and rcs. Linear and dBi forms are separate outputs; gprMax does not silently convert one into the other.

For the range-normalized electric field, radiation intensity is

\[U(\theta,\phi,f) = \frac{|F_\theta|^2+|F_\phi|^2}{2\eta},\]

where \(\eta\) is the wave impedance of the homogeneous material around the NTFF surface. When directivity or either efficiency is requested, gprMax also evaluates a temporary full sphere using Gauss–Legendre quadrature in \(\cos\theta\) and periodic quadrature in \(\phi\):

\[P_\mathrm{rad}(f) = \int_{4\pi}U(\theta,\phi,f)\,\mathrm{d}\Omega, \qquad D(\theta,\phi,f) = \frac{4\pi U(\theta,\phi,f)}{P_\mathrm{rad}(f)}.\]

The quadrature order is selected from the largest requested value of \(ka\), where \(a\) is the bounding radius of the completed integration surface. The temporary full-sphere fields are processed in blocks and are not stored; only the radiated power, estimated pattern maximum and its direction, and quadrature metadata are retained. The maximum estimate is additionally refined with the directions explicitly requested for that output, so a fine user grid cannot report a larger directivity than the stored maximum. Therefore a user may request only a principal-plane cut and still obtain correctly full-sphere-normalised directivity.

For an associated antenna port set, the exact-frequency terminal spectra give

\[P_\mathrm{acc}(f) = \sum_p \frac{1}{2}\Re\{V_p I_p^*\}, \qquad P_\mathrm{inc}(f) = \sum_p \frac{|V_p^+|^2}{2Z_{0p}}.\]

The requested gain quantities are

\[G = \frac{4\pi U}{P_\mathrm{acc}}, \qquad G_\mathrm{realized} = \frac{4\pi U}{P_\mathrm{inc}},\]

and the scalar efficiencies stored for each frequency are

\[\eta_\mathrm{rad}=\frac{P_\mathrm{rad}}{P_\mathrm{acc}}, \qquad \eta_\mathrm{total}=\frac{P_\mathrm{rad}}{P_\mathrm{inc}}.\]

All surface and terminal spectra use the same engineering DFT and transform scale, which cancels in these dimensionless ratios. Frequencies below -40 dB of the peak total incident spectrum, invalid terminal samples, or non-positive normalising powers are written as NaN and marked invalid in the HDF5 port metadata. A radiation efficiency materially above unity emits a warning and normally indicates an insufficient time window, mesh error, integration-surface error, or inconsistent port definition.

rcs requests bistatic radar cross section and requires a TFSF plane-wave source. The NTFF surface must strictly enclose the TFSF box and be clear of its field-correction stencil. With hash commands, exactly one plane wave must be present and it is associated automatically. Use one plane wave per simulation for an unambiguous RCS result. RCS and port-normalised gain belong to different excitation workflows and cannot be combined in one result.

Unlike the exact spherical receiver commands, #ksir_far_field has no radius. Each field component is the range-normalized quantity

\[F_\mathrm{s}(\theta,\phi,f) = r\,\exp(+jkr)\,E_\mathrm{s}(r,\theta,\phi,f),\]

in the far-zone limit, where the subscript s denotes the scattered field. The RCS is

\[\sigma(\theta,\phi,f) = 4\pi \frac{|F_{\mathrm{s},\theta}|^2+|F_{\mathrm{s},\phi}|^2} {|E_{\mathrm{inc},x}|^2+|E_{\mathrm{inc},y}|^2 +|E_{\mathrm{inc},z}|^2}.\]

The incident spectrum is not inferred from the nominal waveform amplitude. gprMax samples the actual numerically propagated field of the associated discrete plane wave and transforms it using the same frequencies and time window as the NTFF surface data. Plane-wave start and stop times are therefore included automatically. Frequencies at which the incident spectrum is zero produce NaN; results where it is very small can be poorly conditioned and should not be used.

rcs is stored as a real, linear quantity in square metres, not in dBsm. It can be converted using

\[\sigma_\mathrm{dBsm} = 10\log_{10}\!\left(\frac{\sigma}{1\,\mathrm{m}^2}\right).\]

The requested \((\theta,\phi)\) specifies the observation direction. For monostatic RCS it must be opposite to the incident propagation direction. For example, a plane wave propagating along +x is observed in backscatter at theta=90 and phi=180 degrees. The far-field normalization and engineering phase convention are also written as HDF5 attributes.

#ksir_far_field: 90 180 antenna_band backscatter Etheta Ephi rcs

#ksir_far_field_array:

Requests the Cartesian product of inclusive theta and phi ranges:

#ksir_far_field_array: theta1 theta2 dtheta phi1 phi2 dphi transform_id [output_id [output1 output2 ...]]

Each range must contain an integer number of positive increments. For example, the following requests a five-degree full-sphere pattern:

#ksir_far_field_array: 0 180 5 0 360 5 antenna_band pattern Etheta Ephi radiation_intensity

#ntff_far_field: and #ntff_far_field_array:

Request conventional equivalent-current frequency-domain far fields:

#ntff_far_field: theta phi transform_id [output_id [output1 output2 ...]]
#ntff_far_field_array: theta1 theta2 dtheta phi1 phi2 dphi transform_id [output_id [output1 output2 ...]]

transform_id must refer to #ntff_frequency or #ntff_layered_frequency. Angles, default field components, range and increment rules, range normalization, derived radiation quantities, RCS, and antenna-port normalization are identical to the corresponding #ksir_far_field commands. Because the underlying surface integral is independent, requesting both formulations on the same #ntff_surface provides a useful numerical cross-check.

Only a planar-layered transform accepts exterior_power, exterior_efficiency, and exterior_maximum. Their definitions and port requirements are given under #ntff_layered_frequency above.

#ntff_far_field_array: 0 180 5 0 360 5 current_band current_pattern Etheta Ephi directivity_dbi

Symmetry-completed surface example

A symmetry plane can coincide with an integration-surface face. No closure option is entered on the KSIR command:

#symmetry_boundary: x0 pmc
#ntff_surface: 0 0.034 0.034 0.026 0.066 0.066 half_surface
#ksir_time_rx: 0.04 0.05 0.051 half_surface half_fields Ez

The physical x0 face is not sampled. The other five faces and their reflected images form the completed closed surface. Observation points must be outside this completed physical-plus-image surface, not merely outside the simulated half. This KSIR workflow is supported by the local CPU, CUDA, OpenCL, and Metal solvers for nondispersive models. Equivalent-current transforms do not yet support symmetry image completion; physical faces can instead be omitted explicitly for an open frequency-domain Huygens surface. OpenCL has end-to-end qualification on the development server. Metal has source-generation and dispatch coverage for symmetry completion, but this workflow still requires execution qualification on suitable Apple hardware. The homogeneous frequency-domain KSIR integration test does not exercise symmetry completion or this time-domain example. Hosted CI does not run Apple GPU tests.

PML commands

The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain.

#pml_cells:

Allows you to control the number of cells (thickness) of PML that are used on the six sides of the model domain. The PML is defined within the model domain, i.e. it is not added to the domain size. The syntax of the command is:

#pml_cells: i1 [i2 i3 i4 i5 i6]

• i1 is the number of cells of PML to use on all sides of the model domain (can be set to zero to completely switch off the PML), or i1 is the number of cells of PML to use on the side of the model domain nearest the origin of the x-axis (x0).

• i2 is the number of cells of PML to use on the side of the model domain nearest the origin of the y-axis (y0).

• i3 is the number of cells of PML to use on the side of the model domain nearest the origin of the z-axis (z0).

• i4 is the number of cells of PML to use on the side of the model domain furthest from the origin of the x-axis (xmax).

• i5 is the number of cells of PML to use on the side of the model domain furthest from the origin of the y-axis (ymax).

• i6 is the number of cells of PML to use on the side of the model domain furthest from the origin of the z-axis (zmax).

• i1 i2 i3 i4 i5 i6 may be set to zero to turn off the PML on a specific side of the model domain.

For example to use a PML with 20 cells (thicker than the default 10 cells) on only the z-axis sides of the domain use:

#pml_cells: 10 10 20 10 10 20

#pml_formulation:

Allows you to alter the formulation used for the PML. The current options are to use the Higher Order RIPML (HORIPML) - https://doi.org/10.1109/TAP.2011.2180344, or Multipole RIPML (MRIPML) - https://doi.org/10.1109/TAP.2018.2823864. The syntax of the command is:

#pml_formulation: str1 [str2]

• str1 can be either HORIPML or MRIPML.

• str2 is an optional reusable PML profile ID. Without an ID the formulation is applied globally. With an ID it is stored for use by a #pml_slab and does not change the domain-boundary PMLs.

For example to use the Multipole RIPML:

#pml_formulation: MRIPML

#pml_cfs:

Allows you (advanced) control of the parameters that are used to build each order of the PML. Up to a second order PML can currently be specified, i.e. by using two #pml_cfs commands. The syntax of the command is:

#pml_cfs: str1 str2 f1 f2 str3 str4 f3 f4 str5 str6 f5 f6 [str7]

• str1 is the type of scaling to use for the CFS \(\alpha\) parameter. It can be constant, linear, quadratic, cubic, quartic, quintic and sextic.

• str2 is the direction of the scaling to use for the CFS \(\alpha\) parameter. It can be forward or reverse.

• f1 f2 are the minimum and maximum values for the CFS \(\alpha\) parameter.

• str3 is the type of scaling to use for the CFS \(\kappa\) parameter. It can be constant, linear, quadratic, cubic, quartic, quintic and sextic.

• str4 is the direction of the scaling to use for the CFS \(\kappa\) parameter. It can be forward or reverse.

• f3 f4 are the minimum and maximum values for the CFS \(\kappa\) parameter. The minimum value for the CFS \(\kappa\) parameter is one.

• str5 is the type of scaling to use for the CFS \(\sigma\) parameter. It can be constant, linear, quadratic, cubic, quartic, quintic and sextic.

• str6 is the direction of the scaling to use for the CFS \(\sigma\) parameter. It can be forward or reverse.

• f5 f6 are the minimum and maximum values for the CFS \(\sigma\) parameter.

• str7 is an optional reusable PML profile ID. It must match the ID on a named #pml_formulation and may then be selected by #pml_slab.

The CFS values (which are internally specified) used for the default standard first order PML are: #pml_cfs: constant forward 0 0 constant forward 1 1 quartic forward 0 None. Specifying ‘None’ for the maximum value of \(\sigma\) forces gprMax to calculate it internally based on the relative permittivity and permeability of the underlying materials in the model.

The parameters will be applied to all slabs of the PML that are switched on. When a profile ID is supplied, they are applied only to slabs that select that profile. A profile can contain one or two CFS terms. A named formulation with no named #pml_cfs uses the default first-order CFS parameters.

Do not duplicate the unshifted first-order defaults to construct a second-order HORIPML. The product can have a negative real stretch and amplify evanescent fields even far below the CFL time-step limit. See Selecting a stable second-order HORIPML profile for the classical/CFS shift condition and a tested repair with matching spatial grading.

Do not duplicate the unshifted first-order defaults to construct a second-order HORIPML. The product can have a negative real stretch and amplify evanescent fields even far below the CFL time-step limit. See Selecting a stable second-order HORIPML profile for the classical/CFS shift condition and a tested repair with matching spatial grading.

Tip

forward direction implies minimum parameter value at the inner boundary of the PML and maximum parameter value at the edge of computational domain, reverse is the opposite.

#pml_slab:

Places an experimental, one-axis RIPML correction in an axis-aligned region of the main 3D grid. This can form a local matched load inside a PEC guiding structure, or replace one domain-boundary PML with an independently configured profile. The syntax is:

#pml_slab: f1 f2 f3 f4 f5 f6 str1 [str2 [str3]]

• f1 f2 f3 and f4 f5 f6 are the lower and upper corners of the slab.

• str1 is the maximum-stretch face: x0, y0, z0, xmax, ymax, or zmax. The opposite face is the zero-stretch entrance.

• str2 is an optional profile ID defined by #pml_formulation and, optionally, #pml_cfs. If omitted, the global PML configuration is used.

• str3 optionally controls automatic PEC enclosure generation: y (default) or n. To specify str3 without a profile, use None for str2.

For example, this replaces a disabled x0 boundary PML with a 12-cell slab using a locally defined MRIPML recipe:

#pml_cells: 0 10 10 10 10 10
#pml_formulation: MRIPML port_load
#pml_slab: 0 0 0 0.012 0.080 0.060 x0 port_load

For an internal absorber, gprMax normally creates five PEC plates after all user geometry has been processed: four transverse walls and a backing plate on the maximum-stretch face. Plates coincident with model-domain faces are not needed and are omitted. The zero-stretch entrance remains open. Existing geometry on an automatically enclosed face is therefore made PEC.

Advanced experiments may disable this behaviour while retaining the global PML profile:

#pml_slab: 0.010 0.020 0.020 0.022 0.060 0.050 x0 None n

With automatic plates disabled, gprMax checks the final Yee edges and warns about every exposed transverse or maximum-stretch face, but permits the model to run for research and custom-enclosure experiments. Incomplete enclosures have no stability guarantee: long-duration testing has demonstrated exponential late-time growth for exposed transverse and open maximum-stretch terminations. The material cross-section must remain constant along the absorption direction; inconsistently filled slabs, same-axis slab overlaps, and overlaps with a native PML remain errors because they invalidate the current formulation.

An automatically generated internal identifier is reported in the log. The hash command defines a slab on the main 3D grid. With the Python API a slab may instead be added to an HSG subgrid, provided the complete slab lies inside its working region and does not overlap the HSG coupling or auxiliary-PML regions. Subgrid slabs use the CPU or CUDA solver and the subgrid’s local spatial and temporal discretisation. Main-grid domain-decomposed MPI CPU models are also supported. A slab may cross any number of rank boundaries. Its CFS grading is evaluated over the complete global thickness and sliced consistently between ranks, while each rank allocates only the PML history arrays for its local part. Consequently, the timestep loop needs no slab-specific communication in addition to the normal electric- and magnetic-field halo exchanges. Automatic PEC enclosure, material-extrusion validation, custom profiles, and boundary-replacement slabs have the same behaviour as in serial models. Orthogonal PML slabs may overlap, as domain PMLs do at edges and corners.

#symmetry_boundary:

Sets a PEC or PMC symmetry boundary on one face of the model domain, replacing the PML on that face. The command may be used more than once to set different faces. The syntax is:

#symmetry_boundary: str1 str2

• str1 is x0, y0, z0, xmax, ymax, or zmax.

• str2 is pec or pmc.

For example:

#symmetry_boundary: x0 pec
#symmetry_boundary: ymax pmc

Note

• The PML thickness on a symmetry face is set to zero automatically.

• PEC and PMC boundaries, including PMC boundaries in models containing dispersive materials, are supported by the CPU, CUDA, OpenCL, Metal, and domain-decomposed MPI CPU solvers.

• In MPI models, each face is constructed and updated only by ranks that touch that global domain face. Domain-edge corrections are likewise restricted to ranks touching both physical faces; internal rank seams are not treated as symmetry edges.

• Symmetry boundaries are not currently supported in 2D mode or on a subgrid. They may be used on the main grid of a model that contains subgrids.