diff --git a/docs/source/features.rst b/docs/source/features.rst index 51e05d59..e1dc7b7f 100644 --- a/docs/source/features.rst +++ b/docs/source/features.rst @@ -16,7 +16,7 @@ The code has been completely re-written in Python/Cython. In the process a lot o * Scripting in the input file. * Built-in library of antenna models * Anisotropic material modelling. -* Dispersive material modelling using multiple pole Debye, Lorenz or Drude formulations. +* Dispersive material modelling using multiple pole Debye, Lorentz or Drude formulations. * Building heterogeneous objects using fractal distributions. * Building objects with rough surfaces. * Modelling soils with realistic dielectric and geometric properties. @@ -28,7 +28,7 @@ New commands * ``#python`` and ``#end_python`` are used to define blocks of the input file where Python code will be executed. This allows the user to use scripting directly in the input file. * ``#material`` replaces ``#medium`` with a new syntax. * ``#add_dispersion_debye`` is used to add Debye dispersive properties to a ``#material``. -* ``#add_dispersion_lorenz`` is used to add Lorenz dispersive properties to a ``#material``. +* ``#add_dispersion_lorentz`` is used to add Lorentz dispersive properties to a ``#material``. * ``#add_dispersion_drude`` is used to add Drude dispersive properties to a ``#material``. * ``#soil_peplinski`` is a soil mixing model that can be used with ``#fractal_box`` to generate soil(s) with more realistic dielectric and geometric properties. * ``#cylindrical_sector`` (like a slice of pie shape) is a new object building command. @@ -95,7 +95,7 @@ The input file has now been made scriptable by permitting blocks of Python code Dispersive media ---------------- -gprMax has always included the ability to represent dispersive materials using a single-pole Debye model. Many materials can be adequately represented using this approach for the typical frequency ranges associated with GPR. However, multi-pole Debye, Drude and Lorenz functions are often used to simulate the electric susceptibility of materials such as: water [PIE2009]_, human tissue [IRE2013]_, cold plasma [LI2013]_, gold [VIA2005]_, and soils [BER1998]_, [GIAK2012]_, [TEI1998]_. Electric susceptibility relates the polarization density to the electric field, and includes both the real and imaginary parts of the complex electric permittivity variation. In the new version of gprMax a recursive convolution based method is used to express dispersive properties as apparent current density sources [GIA2014]_. A major advantage of this implementation is that it creates an inclusive susceptibility function that holds, as special cases, Debye, Drude and Lorenz materials. For further details see the :ref:`material commands section `. +gprMax has always included the ability to represent dispersive materials using a single-pole Debye model. Many materials can be adequately represented using this approach for the typical frequency ranges associated with GPR. However, multi-pole Debye, Drude and Lorentz functions are often used to simulate the electric susceptibility of materials such as: water [PIE2009]_, human tissue [IRE2013]_, cold plasma [LI2013]_, gold [VIA2005]_, and soils [BER1998]_, [GIAK2012]_, [TEI1998]_. Electric susceptibility relates the polarization density to the electric field, and includes both the real and imaginary parts of the complex electric permittivity variation. In the new version of gprMax a recursive convolution based method is used to express dispersive properties as apparent current density sources [GIA2014]_. A major advantage of this implementation is that it creates an inclusive susceptibility function that holds, as special cases, Debye, Drude and Lorentz materials. For further details see the :ref:`material commands section `. Realistic soils, heterogeneous objects and rough surfaces --------------------------------------------------------- diff --git a/docs/source/input.rst b/docs/source/input.rst index 18f4e996..57e9b158 100644 --- a/docs/source/input.rst +++ b/docs/source/input.rst @@ -314,10 +314,10 @@ For example to create a model of water with a single Debye pole, :math:`\epsilon * Temporal values associated with pole frequencies and relaxation times should always be greater than the time step :math:`\Delta t` used in the model. -#add_dispersion_lorenz: +#add_dispersion_lorentz: ----------------------- -Allows you to add dispersive properties to an already defined ``#material`` based on a multiple pole Lorenz formulation (see :ref:`capabilities` section). For example, the susceptability function for a single-pole Lorentz material is given by: +Allows you to add dispersive properties to an already defined ``#material`` based on a multiple pole Lorentz formulation (see :ref:`capabilities` section). For example, the susceptability function for a single-pole Lorentz material is given by: .. math:: @@ -335,21 +335,21 @@ The syntax of the command is: .. code-block:: none - #add_dispersion_lorenz: i1 f1 f2 f3 f4 f5 f6 ... str1 + #add_dispersion_lorentz: i1 f1 f2 f3 f4 f5 f6 ... str1 -* ``i1`` is the number of Lorenz poles. -* ``f1`` is the difference between the zero-frequency relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_{rp1} = \epsilon_{rsp1} - \epsilon_{r \infty p1}` , for the first Lorenz pole. -* ``f2`` is the frequency (Hertz), :math:`\omega_{p1}`, for the first Lorenz pole. -* ``f3`` is the damping coefficient (Hertz), :math:`\delta_{p1}`, for the first Lorenz pole. -* ``f4`` is the difference between the zero-frequency relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_{rp2} = \epsilon_{rsp2} - \epsilon_{r \infty p2}` , for the second Lorenz pole. -* ``f5`` is the frequency (Hertz), :math:`\omega_{p2}`, for the second Lorenz pole. -* ``f6`` is the damping coefficient (Hertz), :math:`\delta_{p2}`, for the second Lorenz pole. +* ``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. :math:`\Delta \epsilon_{rp1} = \epsilon_{rsp1} - \epsilon_{r \infty p1}` , for the first Lorentz pole. +* ``f2`` is the frequency (Hertz), :math:`\omega_{p1}`, for the first Lorentz pole. +* ``f3`` is the damping coefficient (Hertz), :math:`\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. :math:`\Delta \epsilon_{rp2} = \epsilon_{rsp2} - \epsilon_{r \infty p2}` , for the second Lorentz pole. +* ``f5`` is the frequency (Hertz), :math:`\omega_{p2}`, for the second Lorentz pole. +* ``f6`` is the damping coefficient (Hertz), :math:`\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 :math:`\Delta \epsilon_{rp}`, :math:`\omega_p` and :math:`\delta_p` for as many Lorenz poles as you have specified with ``i1``. + * You can continue to add triplets of values for :math:`\Delta \epsilon_{rp}`, :math:`\omega_p` and :math:`\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. :math:`\epsilon_{r \infty}`. * Temporal values associated with pole frequencies and relaxation times should always be greater than the time step :math:`\Delta t` used in the model. diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index 0cabef72..c40c251b 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -130,7 +130,7 @@ def check_cmd_names(processedlines): singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#time_step_limit_type', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps'], 'None') # Commands that there can be multiple instances of in a model - these will be lists within the dictionary - multiplecmds = {key: [] for key in ['#geometry_view', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorenz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#rx', '#rx_box', '#snapshot', '#pml_cfs']} + multiplecmds = {key: [] for key in ['#geometry_view', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#rx', '#rx_box', '#snapshot', '#pml_cfs']} # Geometry object building commands that there can be multiple instances of in a model - these will be lists within the dictionary geometrycmds = ['#edge', '#plate', '#triangle', '#box', '#sphere', '#cylinder', '#cylindrical_sector', '#fractal_box', '#add_surface_roughness', '#add_surface_water', '#add_grass'] diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 52feab05..8e19425b 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -432,7 +432,7 @@ def process_multicmds(multicmds, G): if G.messages: print('Debye-type disperion added to {} with delta_epsr={}, and tau={} secs created.'.format(material.ID, ','.join('%4.2f' % deltaer for deltaer in material.deltaer), ','.join('%4.3e' % tau for tau in material.tau))) - cmdname = '#add_dispersion_lorenz' + cmdname = '#add_dispersion_lorentz' if multicmds[cmdname] != 'None': for cmdinstance in multicmds[cmdname]: tmp = cmdinstance.split() @@ -452,7 +452,7 @@ def process_multicmds(multicmds, G): raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' material(s) {} do not exist'.format(notfound)) for material in materials: - material.type = 'lorenz' + material.type = 'lorentz' material.poles = poles material.average = False for pole in range(1, 3 * poles, 3): @@ -461,12 +461,12 @@ def process_multicmds(multicmds, G): material.tau.append(float(tmp[pole + 1])) material.alpha.append(float(tmp[pole + 2])) else: - raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the permittivity difference and relaxation times, and relaxation times that are greater than the time step for the model.') + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the permittivity difference and frequencies, and associated times that are greater than the time step for the model.') if material.poles > Material.maxpoles: Material.maxpoles = material.poles if G.messages: - print('Lorenz-type disperion added to {} with delta_epsr={}, tau={} secs, and alpha={} created.'.format(material.ID, ','.join('%4.2f' % deltaer for deltaer in material.deltaer), ','.join('%4.3e' % tau for tau in material.tau), ','.join('%4.3e' % alpha for alpha in material.alpha))) + print('Lorentz-type disperion added to {} with delta_epsr={}, omega={} secs, and gamma={} created.'.format(material.ID, ','.join('%4.2f' % deltaer for deltaer in material.deltaer), ','.join('%4.3e' % tau for tau in material.tau), ','.join('%4.3e' % alpha for alpha in material.alpha))) cmdname = '#add_dispersion_drude' @@ -497,12 +497,12 @@ def process_multicmds(multicmds, G): material.tau.append(float(tmp[pole ])) material.alpha.append(float(tmp[pole + 1])) else: - raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the relaxation times, and relaxation times that are greater than the time step for the model.') + raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires positive values for the frequencies, and associated times that are greater than the time step for the model.') if material.poles > Material.maxpoles: Material.maxpoles = material.poles if G.messages: - print('Drude-type disperion added to {} with tau1={} secs, and tau2={} secs created.'.format(material.ID, ','.join('%4.3e' % tau for tau in material.tau), ','.join('%4.3e' % alpha for alpha in material.alpha))) + print('Drude-type disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ','.join('%4.3e' % tau for tau in material.tau), ','.join('%4.3e' % alpha for alpha in material.alpha))) cmdname = '#soil_peplinski' diff --git a/gprMax/materials.py b/gprMax/materials.py index 6f332098..72408563 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -28,7 +28,7 @@ class Material(): maxpoles = 0 # Types of material - types = ['standard', 'debye', 'lorenz', 'drude'] + types = ['standard', 'debye', 'lorentz', 'drude'] # Properties of water from: http://dx.doi.org/10.1109/TGRS.2006.873208 waterer = 80.1 @@ -102,12 +102,16 @@ class Material(): if self.type == 'debye': self.w[x] = self.deltaer[x] / self.tau[x] self.q[x] = -1 / self.tau[x] - elif self.type == 'lorenz': - wp2 = (2 * np.pi * (1 / self.tau[x])) * (2 * np.pi * (1 / self.tau[x])) + elif self.type == 'lorentz': + # tau for Lorentz materials are pole frequencies + # alpha for Lorentz materials are the damping coefficients + wp2 = (2 * np.pi * self.tau[x]) * (2 * np.pi * (1 / self.tau[x])) self.w[x] = -(wp2 * self.deltaer[x]) * j / np.sqrt(wp2 - (self.alpha[x] * self.alpha[x])) self.q[x] = -self.alpha[x] + np.sqrt(wp2 - (self.alpha[x] * self.alpha[x])) * j elif self.type == 'drude': - wp2 = (2 * np.pi * (1 / self.tau[x])) * (2 * np.pi * (1 / self.tau[x])) + # tau for Drude materials are pole frequencies + # alpha for Drude materials are the inverse of relaxation times + wp2 = (2 * np.pi * self.tau[x]) * (2 * np.pi * self.tau[x]) self.se += wp2 / self.alpha[x] self.w[x] = - (wp2 / self.alpha[x]) self.q[x] = - self.alpha[x] diff --git a/tools/inputfile_new2old.py b/tools/inputfile_new2old.py index 44a4318b..10d94180 100644 --- a/tools/inputfile_new2old.py +++ b/tools/inputfile_new2old.py @@ -43,7 +43,7 @@ newfile += '_oldstyle' print("Attempting to convert inputfile '{}' to use old syntax...\n".format(inputfile)) -newcommands = ['#add_dispersion_lorenz', '#add_dispersion_drude', '#fractal_box', '#soil_peplinski', '#python', '#end_python', '#add_surface_roughness', '#add_surface_water', '#add_grass', '#magnetic_dipole', '#pml_cfs', '#cylindrical_sector', '#time_step_limit_type'] +newcommands = ['#add_dispersion_lorentz', '#add_dispersion_drude', '#fractal_box', '#soil_peplinski', '#python', '#end_python', '#add_surface_roughness', '#add_surface_water', '#add_grass', '#magnetic_dipole', '#pml_cfs', '#cylindrical_sector', '#time_step_limit_type'] materials = {} debyes = [] waveforms = []