diff --git a/examples/antenna_like_GSSI_1500.vti b/examples/antenna_like_GSSI_1500.vti
new file mode 100644
index 00000000..56e7a075
Binary files /dev/null and b/examples/antenna_like_GSSI_1500.vti differ
diff --git a/examples/antenna_like_GSSI_1500_fs.h5 b/examples/antenna_like_GSSI_1500_fs.h5
new file mode 100644
index 00000000..0012b8b2
Binary files /dev/null and b/examples/antenna_like_GSSI_1500_fs.h5 differ
diff --git a/examples/antenna_like_GSSI_1500_fs.py b/examples/antenna_like_GSSI_1500_fs.py
index 7ae28dc2..0ee0453b 100644
--- a/examples/antenna_like_GSSI_1500_fs.py
+++ b/examples/antenna_like_GSSI_1500_fs.py
@@ -41,7 +41,7 @@ gv2 = gprMax.GeometryView(p1=(ant_pos[0] - 0.170/2, ant_pos[1] - 0.108/2, ant_po
p2=(ant_pos[0] + 0.170/2, ant_pos[1] + 0.108/2, ant_pos[2] + 0.010),
dl=(dl, dl, dl), filename='antenna_like_GSSI_1500_pcb',
output_type='f')
-#scene.add(gv1)
+scene.add(gv1)
#scene.add(gv2)
# Run model
diff --git a/examples/heterogeneous_soil_v2.h5 b/examples/heterogeneous_soil_v2.h5
new file mode 100644
index 00000000..a5ea069d
Binary files /dev/null and b/examples/heterogeneous_soil_v2.h5 differ
diff --git a/examples/heterogeneous_soil_v2.in b/examples/heterogeneous_soil_v2.in
new file mode 100644
index 00000000..d9cdc811
--- /dev/null
+++ b/examples/heterogeneous_soil_v2.in
@@ -0,0 +1,34 @@
+#title: Heterogeneous soil using a stochastic distribution of dielectric properties given by a mixing model from Peplinski
+#domain: 0.15 0.15 0.15
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 1e-9
+
+#waveform: ricker 1 1.5e9 my_ricker
+#hertzian_dipole: z 0.045 0.075 0.07 my_ricker
+#rx: 0.105 0.075 0.07
+
+#material: 3 0 1 0 sand
+#material: 4 0.1 1 0 sand2
+#material: 6 0.001 1 0 asphalt
+#material: 4.9 0.001 1 0 my_water
+#add_dispersion_debye: 1 75.2 9.231e-12 my_water
+
+ellipsoid: 0.045 0.045 0.045 0.03 0.03 0.03 sand
+sphere: 0.075 0.075 0.075 0.03 sand2
+box: 0.029 0.029 0.029 0.03 0.03 0.03 sand
+box: 0.08 0.08 0.08 0.081 0.081 0.081 sand
+cylinder: 0.03 0.03 0.03 0.08 0.08 0.08 0.03 sand
+cone: 0.02 0.075 0.075 0.08 0.075 0.075 0.03 0.0 sand
+cone: 0.075 0.02 0.075 0.075 0.08 0.075 0.03 0.0 sand
+cone: 0.075 0.075 0.02 0.075 0.075 0.08 0.03 0.0 sand
+cone: 0.03 0.03 0.03 0.08 0.08 0.08 0.028 0.005 free_space
+
+
+
+#material_list: sand pec sand2 asphalt my_water my_list
+#material_range: 2 6 0.1 0.001 1.0 1.0 0.0 0.0 test
+soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil
+#fractal_box: 0 0 0 0.15 0.15 0.1 1.5 1 1 1 5 test my_soil_box
+add_surface_roughness: 0 0 0.070 0.15 0.15 0.070 1.5 1 1 0.065 0.080 my_soil_box
+
+#geometry_view: 0 0 0 0.15 0.15 0.15 0.001 0.001 0.001 heterogeneous_soil_v2 n
diff --git a/examples/heterogeneous_soil_v2.vti b/examples/heterogeneous_soil_v2.vti
new file mode 100644
index 00000000..ad297809
Binary files /dev/null and b/examples/heterogeneous_soil_v2.vti differ
diff --git a/gprMax/__init__.py b/gprMax/__init__.py
index 027f21bd..b5ccc50b 100644
--- a/gprMax/__init__.py
+++ b/gprMax/__init__.py
@@ -26,7 +26,7 @@ from .cmds_multiuse import (PMLCFS, AddDebyeDispersion, AddDrudeDispersion,
AddLorentzDispersion, GeometryObjectsWrite,
GeometryView, HertzianDipole, MagneticDipole,
Material, Rx, RxArray, Snapshot, SoilPeplinski,
- TransmissionLine, VoltageSource, Waveform)
+ TransmissionLine, VoltageSource, Waveform, MaterialRange, MaterialList)
from .cmds_singleuse import (Discretisation, Domain, ExcitationFile,
OMPThreads, PMLProps, RxSteps, SrcSteps,
TimeStepStabilityFactor, TimeWindow, Title)
diff --git a/gprMax/cmds_geometry/box.py b/gprMax/cmds_geometry/box.py
index 40311f21..aa509bf7 100644
--- a/gprMax/cmds_geometry/box.py
+++ b/gprMax/cmds_geometry/box.py
@@ -132,8 +132,8 @@ class Box(UserObjectGeometry):
materials[2].se), axis=0)
m.mr = np.mean((materials[0].mr, materials[1].mr,
materials[2].mr), axis=0)
- m.sm = np.mean((materials[0].mr, materials[1].mr,
- materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
# Append the new material object to the materials list
grid.materials.append(m)
diff --git a/gprMax/cmds_geometry/cone.py b/gprMax/cmds_geometry/cone.py
new file mode 100644
index 00000000..34e5ec13
--- /dev/null
+++ b/gprMax/cmds_geometry/cone.py
@@ -0,0 +1,147 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax. If not, see .
+
+import logging
+
+import numpy as np
+
+from ..cython.geometry_primitives import build_cone
+from ..materials import Material
+from .cmds_geometry import UserObjectGeometry
+
+logger = logging.getLogger(__name__)
+
+
+class Cone(UserObjectGeometry):
+ """Introduces a circular cone into the model. The difference with the cylinder is that the faces of the cone
+ can have different radii and one of them can be zero.
+
+ Attributes:
+ p1: list of the coordinates (x,y,z) of the centre of the first face
+ of the cone.
+ p2: list of the coordinates (x,y,z) of the centre of the second face
+ of the cone.
+ r1: float of the radius of the first face of the cone.
+ r2: float of the radius of the second face of the cone.
+ material_id: string for the material identifier that must correspond
+ to material that has already been defined.
+ material_ids: list of material identifiers in the x, y, z directions.
+ averaging: string (y or n) used to switch on and off dielectric smoothing.
+ """
+
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
+ self.hash = '#cylinder'
+
+ def create(self, grid, uip):
+ try:
+ p1 = self.kwargs['p1']
+ p2 = self.kwargs['p2']
+ r1 = self.kwargs['r1']
+ r2 = self.kwargs['r2']
+ except KeyError:
+ logger.exception(self.__str__() + ' please specify 2 points and two radii')
+ raise
+
+ # Check averaging
+ try:
+ # Try user-specified averaging
+ averagecylinder = self.kwargs['averaging']
+ except KeyError:
+ # Otherwise go with the grid default
+ averagecylinder = grid.averagevolumeobjects
+
+ # Check materials have been specified
+ # Isotropic case
+ try:
+ materialsrequested = [self.kwargs['material_id']]
+ except KeyError:
+ # Anisotropic case
+ try:
+ materialsrequested = self.kwargs['material_ids']
+ except KeyError:
+ logger.exception(self.__str__() + ' no materials have been specified')
+ raise
+
+ p3 = uip.round_to_grid_static_point(p1)
+ p4 = uip.round_to_grid_static_point(p2)
+
+ x1, y1, z1 = uip.round_to_grid(p1)
+ x2, y2, z2 = uip.round_to_grid(p2)
+
+ if r1 < 0:
+ logger.exception(self.__str__() + f' the radius of the first face {r1:g} should be a positive value.')
+ raise ValueError
+
+ if r2 < 0:
+ logger.exception(self.__str__() + f' the radius of the second face {r2:g} should be a positive value.')
+ raise ValueError
+
+ if r1 == 0 and r2 == 0:
+ logger.exception(self.__str__() + f' not both radii can be zero.')
+ raise ValueError
+
+ # Look up requested materials in existing list of material instances
+ materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+ if len(materials) != len(materialsrequested):
+ notfound = [x for x in materialsrequested if x not in materials]
+ logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+ raise ValueError
+
+ # Isotropic case
+ if len(materials) == 1:
+ averaging = materials[0].averagable and averagecylinder
+ numID = numIDx = numIDy = numIDz = materials[0].numID
+
+ # Uniaxial anisotropic case
+ elif len(materials) == 3:
+ averaging = False
+ numIDx = materials[0].numID
+ numIDy = materials[1].numID
+ numIDz = materials[2].numID
+ requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+ averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+ if averagedmaterial:
+ numID = averagedmaterial.numID
+ else:
+ numID = len(grid.materials)
+ m = Material(numID, requiredID)
+ m.type = 'dielectric-smoothed'
+ # Create dielectric-smoothed constituents for material
+ m.er = np.mean((materials[0].er, materials[1].er,
+ materials[2].er), axis=0)
+ m.se = np.mean((materials[0].se, materials[1].se,
+ materials[2].se), axis=0)
+ m.mr = np.mean((materials[0].mr, materials[1].mr,
+ materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
+
+ # Append the new material object to the materials list
+ grid.materials.append(m)
+
+ build_cone(x1, y1, z1, x2, y2, z2, r1, r2, grid.dx, grid.dy, grid.dz,
+ numID, numIDx, numIDy, numIDz, averaging, grid.solid,
+ grid.rigidE, grid.rigidH, grid.ID)
+
+ dielectricsmoothing = 'on' if averaging else 'off'
+ logger.info(self.grid_name(grid) + f"Cone with face centres {p3[0]:g}m, " +
+ f"{p3[1]:g}m, {p3[2]:g}m and {p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m, " +
+ f"with radii {r1:g}m and {r2:g}, of material(s) {', '.join(materialsrequested)} " +
+ f"created, dielectric smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/cylinder.py b/gprMax/cmds_geometry/cylinder.py
index 1dbcbbc7..d70b824a 100644
--- a/gprMax/cmds_geometry/cylinder.py
+++ b/gprMax/cmds_geometry/cylinder.py
@@ -119,8 +119,8 @@ class Cylinder(UserObjectGeometry):
materials[2].se), axis=0)
m.mr = np.mean((materials[0].mr, materials[1].mr,
materials[2].mr), axis=0)
- m.sm = np.mean((materials[0].mr, materials[1].mr,
- materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
# Append the new material object to the materials list
grid.materials.append(m)
diff --git a/gprMax/cmds_geometry/cylindrical_sector.py b/gprMax/cmds_geometry/cylindrical_sector.py
index 6262bd8b..42005250 100644
--- a/gprMax/cmds_geometry/cylindrical_sector.py
+++ b/gprMax/cmds_geometry/cylindrical_sector.py
@@ -143,8 +143,8 @@ class CylindricalSector(UserObjectGeometry):
materials[2].se), axis=0)
m.mr = np.mean((materials[0].mr, materials[1].mr,
materials[2].mr), axis=0)
- m.sm = np.mean((materials[0].mr, materials[1].mr,
- materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
# Append the new material object to the materials list
grid.materials.append(m)
diff --git a/gprMax/cmds_geometry/ellipsoid.py b/gprMax/cmds_geometry/ellipsoid.py
new file mode 100644
index 00000000..115ee10c
--- /dev/null
+++ b/gprMax/cmds_geometry/ellipsoid.py
@@ -0,0 +1,132 @@
+# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom
+# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
+#
+# This file is part of gprMax.
+#
+# gprMax is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# gprMax is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with gprMax. If not, see .
+
+import logging
+
+import numpy as np
+
+from ..cython.geometry_primitives import build_ellipsoid
+from ..materials import Material
+from .cmds_geometry import UserObjectGeometry
+
+logger = logging.getLogger(__name__)
+
+
+class Ellipsoid(UserObjectGeometry):
+ """Introduces an ellipsoidal object with specific parameters into the model.
+
+ Attributes:
+ p1: list of the coordinates (x,y,z) of the centre of the ellipsoid.
+ xr: float of x semiaxis of the ellipsoid.
+ xy: float of y semiaxis of the ellipsoid.
+ xz: float of z semiaxis of the ellipsoid.
+ material_id: string for the material identifier that must correspond
+ to material that has already been defined.
+ material_ids: list of material identifiers in the x, y, z directions.
+ averaging: string (y or n) used to switch on and off dielectric smoothing.
+ """
+
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
+ self.hash = '#ellipsoid'
+
+ def create(self, grid, uip):
+ try:
+ p1 = self.kwargs['p1']
+ xr = self.kwargs['xr']
+ yr = self.kwargs['yr']
+ zr = self.kwargs['zr']
+
+ except KeyError:
+ logger.exception(self.__str__() + ' please specify a point and the three semiaxes.')
+ raise
+
+ # Check averaging
+ try:
+ # Try user-specified averaging
+ averageellipsoid = self.kwargs['averaging']
+ except KeyError:
+ # Otherwise go with the grid default
+ averageellipsoid = grid.averagevolumeobjects
+
+ # Check materials have been specified
+ # Isotropic case
+ try:
+ materialsrequested = [self.kwargs['material_id']]
+ except KeyError:
+ # Anisotropic case
+ try:
+ materialsrequested = self.kwargs['material_ids']
+ except KeyError:
+ logger.exception(self.__str__() + ' no materials have been specified')
+ raise
+
+ # Centre of sphere
+ p2 = uip.round_to_grid_static_point(p1)
+ #xc, yc, zc = uip.round_to_grid(p1)
+ xc, yc, zc = uip.discretise_point(p1)
+
+ # Look up requested materials in existing list of material instances
+ materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
+
+ if len(materials) != len(materialsrequested):
+ notfound = [x for x in materialsrequested if x not in materials]
+ logger.exception(self.__str__() + f' material(s) {notfound} do not exist')
+ raise ValueError
+
+ # Isotropic case
+ if len(materials) == 1:
+ averaging = materials[0].averagable and averageellipsoid
+ numID = numIDx = numIDy = numIDz = materials[0].numID
+
+ # Uniaxial anisotropic case
+ elif len(materials) == 3:
+ averaging = False
+ numIDx = materials[0].numID
+ numIDy = materials[1].numID
+ numIDz = materials[2].numID
+ requiredID = materials[0].ID + '+' + materials[1].ID + '+' + materials[2].ID
+ averagedmaterial = [x for x in grid.materials if x.ID == requiredID]
+ if averagedmaterial:
+ numID = averagedmaterial.numID
+ else:
+ numID = len(grid.materials)
+ m = Material(numID, requiredID)
+ m.type = 'dielectric-smoothed'
+ # Create dielectric-smoothed constituents for material
+ m.er = np.mean((materials[0].er, materials[1].er,
+ materials[2].er), axis=0)
+ m.se = np.mean((materials[0].se, materials[1].se,
+ materials[2].se), axis=0)
+ m.mr = np.mean((materials[0].mr, materials[1].mr,
+ materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
+
+ # Append the new material object to the materials list
+ grid.materials.append(m)
+
+ build_ellipsoid(xc, yc, zc, xr, yr, zr, grid.dx, grid.dy, grid.dz, numID,
+ numIDx, numIDy, numIDz, averaging, grid.solid,
+ grid.rigidE, grid.rigidH, grid.ID)
+
+ dielectricsmoothing = 'on' if averaging else 'off'
+ logger.info(self.grid_name(grid) + f"Ellipsoid with centre {p2[0]:g}m, " +
+ f"{p2[1]:g}m, {p2[2]:g}m, x-semiaxis {xr:g}m, y-semiaxis {yr:g}m and z-semiaxis {zr:g}m of material(s) " +
+ f"{', '.join(materialsrequested)} created, dielectric " +
+ f"smoothing is {dielectricsmoothing}.")
diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py
index 8eb827a3..55d18277 100644
--- a/gprMax/cmds_geometry/fractal_box.py
+++ b/gprMax/cmds_geometry/fractal_box.py
@@ -141,7 +141,7 @@ class FractalBox(UserObjectGeometry):
raise ValueError
# Create materials from mixing model as number of bins now known
# from fractal_box command.
- mixingmodel.calculate_debye_properties(nbins, grid)
+ mixingmodel.calculate_properties(nbins, grid)
elif not material:
logger.exception(self.__str__() + f' mixing model or material with ' +
'ID {mixing_model_id} does not exist')
diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py
index 887c4c70..5c804f71 100644
--- a/gprMax/cmds_geometry/fractal_box_builder.py
+++ b/gprMax/cmds_geometry/fractal_box_builder.py
@@ -83,7 +83,12 @@ class FractalBoxBuilder(UserObjectGeometry):
volume.fractalvolume *= materialnumID
else:
volume.generate_fractal_volume()
- volume.fractalvolume += volume.mixingmodel.startmaterialnum
+ #volume.fractalvolume += volume.mixingmodel.startmaterialnum
+ for i in range(0, volume.nx):
+ for j in range(0, volume.ny):
+ for k in range(0, volume.nz):
+ numberinbin = volume.fractalvolume[i,j,k]
+ volume.fractalvolume[i,j,k] = volume.mixingmodel.matID[int(numberinbin)]
volume.generate_volume_mask()
@@ -320,7 +325,13 @@ class FractalBoxBuilder(UserObjectGeometry):
raise ValueError
else:
volume.generate_fractal_volume()
- volume.fractalvolume += volume.mixingmodel.startmaterialnum
+ #volume.fractalvolume += volume.mixingmodel.startmaterialnum
+ for i in range(0, volume.nx):
+ for j in range(0, volume.ny):
+ for k in range(0, volume.nz):
+ numberinbin = volume.fractalvolume[i,j,k]
+ volume.fractalvolume[i,j,k] = volume.mixingmodel.matID[int(numberinbin)]
+
data = volume.fractalvolume.astype('int16', order='C')
build_voxels_from_array(volume.xs, volume.ys, volume.zs,
diff --git a/gprMax/cmds_geometry/sphere.py b/gprMax/cmds_geometry/sphere.py
index c7436abb..014ec3b0 100644
--- a/gprMax/cmds_geometry/sphere.py
+++ b/gprMax/cmds_geometry/sphere.py
@@ -72,8 +72,13 @@ class Sphere(UserObjectGeometry):
raise
# Centre of sphere
+
p2 = uip.round_to_grid_static_point(p1)
- xc, yc, zc = uip.round_to_grid(p1)
+ #xc, yc, zc = uip.round_to_grid(p1)
+ xc, yc, zc = uip.discretise_point(p1)
+
+
+
# Look up requested materials in existing list of material instances
materials = [y for x in materialsrequested for y in grid.materials if y.ID == x]
@@ -109,8 +114,8 @@ class Sphere(UserObjectGeometry):
materials[2].se), axis=0)
m.mr = np.mean((materials[0].mr, materials[1].mr,
materials[2].mr), axis=0)
- m.sm = np.mean((materials[0].mr, materials[1].mr,
- materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
# Append the new material object to the materials list
grid.materials.append(m)
diff --git a/gprMax/cmds_geometry/triangle.py b/gprMax/cmds_geometry/triangle.py
index 53d59b41..17ab6077 100644
--- a/gprMax/cmds_geometry/triangle.py
+++ b/gprMax/cmds_geometry/triangle.py
@@ -160,8 +160,8 @@ class Triangle(UserObjectGeometry):
materials[2].se), axis=0)
m.mr = np.mean((materials[0].mr, materials[1].mr,
materials[2].mr), axis=0)
- m.sm = np.mean((materials[0].mr, materials[1].mr,
- materials[2].mr), axis=0)
+ m.sm = np.mean((materials[0].sm, materials[1].sm,
+ materials[2].sm), axis=0)
# Append the new material object to the materials list
grid.materials.append(m)
diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py
index 7a5f639c..31be50ee 100644
--- a/gprMax/cmds_multiuse.py
+++ b/gprMax/cmds_multiuse.py
@@ -31,6 +31,8 @@ from .geometry_outputs import GeometryObjects as GeometryObjectsUser
from .materials import DispersiveMaterial as DispersiveMaterialUser
from .materials import Material as MaterialUser
from .materials import PeplinskiSoil as PeplinskiSoilUser
+from .materials import RangeMaterial as RangeMaterialUser
+from .materials import ListMaterial as ListMaterialUser
from .pml import CFS, CFSParameter
from .receivers import Rx as RxUser
from .snapshots import Snapshot as SnapshotUser
@@ -1438,6 +1440,137 @@ class SoilPeplinski(UserObjectMulti):
grid.mixingmodels.append(s)
+class MaterialRange(UserObjectMulti):
+ """Simple model for varying material properties to create stochastic models
+
+
+ Attributes:
+ er_lower_range: float required for lower relative permittivity value.
+ er_upper_range: float required for upper relative permittivity value.
+ sigma_lower_range: float required for lower conductivity value.
+ sigma_upper_range: float required for upper conductivity value.
+ mr_lower_range: float required for lower relative magnetic permeability value.
+ mr_upper_range: float required for upper relative magnetic permeability value.
+ ro_lower_range: float required for lower magnetic loss value.
+ ro_upper_range: float required for upper magnetic loss value.
+
+ id: string used as identifier for this variable material.
+ """
+
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
+ self.order = 17
+ self.hash = '#material_range'
+
+ def create(self, grid, uip):
+ try:
+ er_lower = self.kwargs['er_lower']
+ er_upper = self.kwargs['er_upper']
+ sigma_lower = self.kwargs['sigma_lower']
+ sigma_upper = self.kwargs['sigma_upper']
+ mr_lower = self.kwargs['mr_lower']
+ mr_upper = self.kwargs['mr_upper']
+ ro_lower = self.kwargs['ro_lower']
+ ro_upper = self.kwargs['ro_upper']
+ ID = self.kwargs['id']
+ except KeyError:
+ logger.exception(self.params_str() + ' requires at exactly nine '
+ 'parameters.')
+ raise
+
+ if er_lower < 1:
+ logger.exception(self.params_str() + ' requires a value greater or equal to 1 '
+ 'for the lower range of relative permittivity.')
+ raise ValueError
+ if mr_lower < 1:
+ logger.exception(self.params_str() + ' requires a value greater or equal to 1 '
+ 'for the lower range of relative magnetic permeability.')
+ raise ValueError
+ if sigma_lower < 0:
+ logger.exception(self.params_str() + ' requires a positive value '
+ 'for the lower limit of conductivity.')
+ raise ValueError
+ if ro_lower < 0:
+ logger.exception(self.params_str() + ' requires a positive value '
+ 'for the lower range magnetic loss.')
+ raise ValueError
+ if er_upper < 1:
+ logger.exception(self.params_str() + ' requires a value greater or equal to 1'
+ 'for the upper range of relative permittivity.')
+ raise ValueError
+
+ if mr_upper < 1:
+ logger.exception(self.params_str() + ' requires a value greater or equal to 1'
+ 'for the upper range of relative magnetic permeability')
+ raise ValueError
+
+ if sigma_upper < 0:
+ logger.exception(self.params_str() + ' requires a positive value '
+ 'for the upper range of conductivity.')
+ raise ValueError
+ if ro_upper < 0:
+ logger.exception(self.params_str() + ' requires a positive value '
+ 'for the upper range of magnetic loss.')
+
+
+ if any(x.ID == ID for x in grid.mixingmodels):
+ logger.exception(self.params_str() + f' with ID {ID} already exists')
+ raise ValueError
+
+ # Create a new instance of the Material class material
+ # (start index after pec & free_space)
+ s = RangeMaterialUser(ID, (er_lower, er_upper), (sigma_lower, sigma_upper), (mr_lower, mr_upper), (ro_lower, ro_upper))
+
+ logger.info(self.grid_name(grid) + 'Material properties used to '
+ f'create {s.ID} with range(s) {s.er[0]:g} to {s.er[1]:g}, relative permittivity '
+ f'{s.sig[0]:g} to {s.sig[1]:g}, S/m conductivity, {s.mu[0]:g} to {s.mu[1]:g} relative magnetic permeability '
+ f'{s.ro[0]:g} to {s.ro[1]:g} Ohm/m magnetic loss, created')
+
+
+ grid.mixingmodels.append(s)
+
+
+class MaterialList(UserObjectMulti):
+ """Simple model for varying material properties to create stochastic models
+
+
+ Attributes:
+ list_of_materials: list of material IDs
+
+ id: string used as identifier for this variable material.
+ """
+
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
+ self.order = 18
+ self.hash = '#material_list'
+
+ def create(self, grid, uip):
+ try:
+ list_of_materials = self.kwargs['list_of_materials']
+ ID = self.kwargs['id']
+ except KeyError:
+ logger.exception(self.params_str() + ' requires at at least 2 '
+ 'parameters.')
+ raise
+
+
+ if any(x.ID == ID for x in grid.mixingmodels):
+ logger.exception(self.params_str() + f' with ID {ID} already exists')
+ raise ValueError
+
+ # Create a new instance of the Material class material
+ # (start index after pec & free_space)
+ s = ListMaterialUser(ID, list_of_materials)
+
+ logger.info(self.grid_name(grid) + 'A list of materials used to '
+ f'create {s.ID} that includes {s.mat}, created')
+
+
+ grid.mixingmodels.append(s)
+
+
+
class GeometryView(UserObjectMulti):
"""Outputs to file(s) information about the geometry (mesh) of model.
diff --git a/gprMax/cython/geometry_primitives.pyx b/gprMax/cython/geometry_primitives.pyx
index 1a4003de..7c359bc7 100644
--- a/gprMax/cython/geometry_primitives.pyx
+++ b/gprMax/cython/geometry_primitives.pyx
@@ -917,6 +917,222 @@ cpdef void build_cylinder(
averaging, solid, rigidE, rigidH, ID)
+
+cpdef void build_cone(
+ float x1,
+ float y1,
+ float z1,
+ float x2,
+ float y2,
+ float z2,
+ float r1,
+ float r2,
+ float dx,
+ float dy,
+ float dz,
+ int numID,
+ int numIDx,
+ int numIDy,
+ int numIDz,
+ bint averaging,
+ np.uint32_t[:, :, ::1] solid,
+ np.int8_t[:, :, :, ::1] rigidE,
+ np.int8_t[:, :, :, ::1] rigidH,
+ np.uint32_t[:, :, :, ::1] ID
+):
+ """Builds cones which sets values in the solid, rigid and ID arrays for
+ a Yee voxel.
+
+ Args:
+ x1, y1, z1, x2, y2, z2: floats for coordinates of the centres of the cone
+ faces.
+ r1: float for radius of the first face of the cone.
+ r2: float for radius of the second face of the cone.
+ dx, dy, dz: floats for spatial discretisation.
+ numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+ averaging: bint for whether material property averaging will occur for
+ the object.
+ solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
+ """
+
+ cdef Py_ssize_t i, j, k
+ cdef int xs, xf, ys, yf, zs, zf, xc, yc, zc
+ cdef float f1f2mag, f2f1mag, f1ptmag, f2ptmag, dot1, dot2, factor1, factor2
+ cdef float theta1, theta2, distance1, distance2, R1, R2
+ cdef float height, distance_axis_1, distance_axis_2
+ cdef bint build, x_align, y_align, z_align
+ cdef np.ndarray f1f2, f2f1, f1pt, f2pt
+ cdef float Rmax
+
+ Rmax = np.amax([r1, r2])
+
+ # Check if cylinder is aligned with an axis
+ x_align = y_align = z_align = 0
+ # x-aligned
+ if (round_value(y1 / dy) == round_value(y2 / dy) and
+ round_value(z1 / dz) == round_value(z2 / dz)):
+ x_align = 1
+
+ # y-aligned
+ elif (round_value(x1 / dx) == round_value(x2 / dx) and
+ round_value(z1 / dz) == round_value(z2 / dz)):
+ y_align = 1
+
+ # z-aligned
+ elif (round_value(x1 / dx) == round_value(x2 / dx) and
+ round_value(y1 / dy) == round_value(y2 / dy)):
+ z_align = 1
+
+ # Calculate a bounding box for the cylinder
+ if x1 < x2:
+ if x_align:
+ xs = round_value(x1 / dx)
+ xf = round_value(x2 / dx)
+ else:
+ xs = round_value((x1 - Rmax) / dx) - 1
+ xf = round_value((x2 + Rmax) / dx) + 1
+ else:
+ if x_align:
+ xs = round_value(x2 / dx)
+ xf = round_value(x1 / dx)
+ else:
+ xs = round_value((x2 - Rmax) / dx) - 1
+ xf = round_value((x1 + Rmax) / dx) + 1
+ if y1 < y2:
+ if y_align:
+ ys = round_value(y1 / dy)
+ yf = round_value(y2 / dy)
+ else:
+ ys = round_value((y1 - Rmax) / dy) - 1
+ yf = round_value((y2 + Rmax) / dy) + 1
+ else:
+ if y_align:
+ ys = round_value(y2 / dy)
+ yf = round_value(y1 / dy)
+ else:
+ ys = round_value((y2 - Rmax) / dy) - 1
+ yf = round_value((y1 + Rmax) / dy) + 1
+ if z1 < z2:
+ if z_align:
+ zs = round_value(z1 / dz)
+ zf = round_value(z2 / dz)
+ else:
+ zs = round_value((z1 - Rmax) / dz) - 1
+ zf = round_value((z2 + Rmax) / dz) + 1
+ else:
+ if z_align:
+ zs = round_value(z2 / dz)
+ zf = round_value(z1 / dz)
+ else:
+ zs = round_value((z2 - Rmax) / dz) - 1
+ zf = round_value((z1 + Rmax) / dz) + 1
+
+ # Set bounds to domain if they outside
+ if xs < 0:
+ xs = 0
+ if xf > solid.shape[0]:
+ xf = solid.shape[0]
+ if ys < 0:
+ ys = 0
+ if yf > solid.shape[1]:
+ yf = solid.shape[1]
+ if zs < 0:
+ zs = 0
+ if zf > solid.shape[2]:
+ zf = solid.shape[2]
+
+ # x-aligned cylinder
+ if x_align:
+ for j in range(ys, yf):
+ for k in range(zs, zf):
+ for i in range(xs, xf):
+ if np.sqrt((j * dy + 0.5 * dy - y1)**2 + (k * dz + 0.5 * dz - z1)**2) <= ((i-xs)/(xf-xs))*(r2-r1) + r1:
+ build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
+ averaging, solid, rigidE, rigidH, ID)
+ # y-aligned cylinder
+ elif y_align:
+ for i in range(xs, xf):
+ for k in range(zs, zf):
+ for j in range(ys, yf):
+ if np.sqrt((i * dx + 0.5 * dx - x1)**2 + (k * dz + 0.5 * dz - z1)**2) <= ((j-ys)/(yf-ys))*(r2-r1) + r1:
+ build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
+ averaging, solid, rigidE, rigidH, ID)
+ # z-aligned cylinder
+ elif z_align:
+ for i in range(xs, xf):
+ for j in range(ys, yf):
+ for k in range(zs, zf):
+ if np.sqrt((i * dx + 0.5 * dx - x1)**2 + (j * dy + 0.5 * dy - y1)**2) <= ((k-zs)/(zf-zs))*(r2-r1) + r1:
+ build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
+ averaging, solid, rigidE, rigidH, ID)
+
+ # Not aligned with any axis
+ else:
+ # Vectors between centres of cylinder faces
+ f1f2 = np.array([x2 - x1, y2 - y1, z2 - z1], dtype=np.float32)
+ f2f1 = np.array([x1 - x2, y1 - y2, z1 - z2], dtype=np.float32)
+
+ # Magnitudes
+ f1f2mag = np.sqrt((f1f2*f1f2).sum(axis=0))
+ f2f1mag = np.sqrt((f2f1*f2f1).sum(axis=0))
+
+ height = f1f2mag
+
+
+ for i in range(xs, xf):
+ for j in range(ys, yf):
+ for k in range(zs, zf):
+ # Build flag - default false, set to True if point is in cylinder
+ build = 0
+ # Vector from centre of first cylinder face to test point
+ f1pt = np.array([i * dx + 0.5 * dx - x1,
+ j * dy + 0.5 * dy - y1,
+ k * dz + 0.5 * dz - z1], dtype=np.float32)
+ # Vector from centre of second cylinder face to test point
+ f2pt = np.array([i * dx + 0.5 * dx - x2,
+ j * dy + 0.5 * dy - y2,
+ k * dz + 0.5 * dz - z2], dtype=np.float32)
+ # Magnitudes
+ f1ptmag = np.sqrt((f1pt*f1pt).sum(axis=0))
+ f2ptmag = np.sqrt((f2pt*f2pt).sum(axis=0))
+ # Dot products
+ dot1 = np.dot(f1f2, f1pt)
+ dot2 = np.dot(f2f1, f2pt)
+
+ if f1ptmag == 0 or f2ptmag == 0:
+ build = 1
+ else:
+ factor1 = dot1 / (f1f2mag * f1ptmag)
+ factor2 = dot2 / (f2f1mag * f2ptmag)
+ # Catch cases where either factor1 or factor2 are 1
+ try:
+ theta1 = np.arccos(factor1)
+ except FloatingPointError:
+ theta1 = 0
+ try:
+ theta2 = np.arccos(factor2)
+ except FloatingPointError:
+ theta2 = 0
+ distance1 = f1ptmag * np.sin(theta1)
+ distance2 = f2ptmag * np.sin(theta2)
+ distance_axis_1 = f1ptmag * np.cos(theta1)
+ distance_axis_2 = f2ptmag * np.cos(theta2)
+ R1 = r1
+ R2 = r2
+
+ if ((distance1 <= (distance_axis_1/height)*(R2 - R1) + R1 or distance2 <= (distance_axis_2/height)*(R1 - R2) + R2) and
+ theta1 <= np.pi/2 and theta2 <= np.pi/2):
+ build = 1
+
+ if build:
+ build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
+ averaging, solid, rigidE, rigidH, ID)
+
+
+
+
+
+
cpdef void build_sphere(
int xc,
int yc,
@@ -983,6 +1199,77 @@ cpdef void build_sphere(
averaging, solid, rigidE, rigidH, ID)
+cpdef void build_ellipsoid(
+ int xc,
+ int yc,
+ int zc,
+ float xr,
+ float yr,
+ float zr,
+ float dx,
+ float dy,
+ float dz,
+ int numID,
+ int numIDx,
+ int numIDy,
+ int numIDz,
+ bint averaging,
+ np.uint32_t[:, :, ::1] solid,
+ np.int8_t[:, :, :, ::1] rigidE,
+ np.int8_t[:, :, :, ::1] rigidH,
+ np.uint32_t[:, :, :, ::1] ID
+):
+ """Builds ellipsoids which sets values in the solid, rigid and ID arrays for
+ a Yee voxel.
+
+ Args:
+ xc, yc, zc: ints for cell coordinates of the centre of the ellipsoid.
+ xr: float for x-semiaxis of the elliposid.
+ yr: float for y-semiaxis of the elliposid.
+ zr: float for z-semiaxis of the elliposid.
+ dx, dy, dz: floats for spatial discretisation.
+ numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
+ averaging: bint for whether material property averaging will occur for
+ the object.
+ solid, rigidE, rigidH, ID: memoryviews to access solid, rigid and ID arrays.
+ """
+
+ cdef Py_ssize_t i, j, k
+ cdef int xs, xf, ys, yf, zs, zf
+
+ # Calculate a bounding box for sphere
+ xs = round_value(((xc * dx) - xr) / dx) - 1
+ xf = round_value(((xc * dx) + xr) / dx) + 1
+ ys = round_value(((yc * dy) - yr) / dy) - 1
+ yf = round_value(((yc * dy) + yr) / dy) + 1
+ zs = round_value(((zc * dz) - zr) / dz) - 1
+ zf = round_value(((zc * dz) + zr) / dz) + 1
+
+ # Set bounds to domain if they outside
+ if xs < 0:
+ xs = 0
+ if xf > solid.shape[0]:
+ xf = solid.shape[0]
+ if ys < 0:
+ ys = 0
+ if yf > solid.shape[1]:
+ yf = solid.shape[1]
+ if zs < 0:
+ zs = 0
+ if zf > solid.shape[2]:
+ zf = solid.shape[2]
+
+ for i in range(xs, xf):
+ for j in range(ys, yf):
+ for k in range(zs, zf):
+ if (((i + 0.5 - xc)**2 * dx**2)/xr**2 +
+ ((j + 0.5 - yc)**2 * dy**2)/yr**2 +
+ ((k + 0.5 - zc)**2 * dz**2)/zr**2 <= 1):
+ build_voxel(i, j, k, numID, numIDx, numIDy, numIDz,
+ averaging, solid, rigidE, rigidH, ID)
+
+
+
cpdef void build_voxels_from_array(
int xs,
int ys,
diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py
index a50f26c7..805835e5 100644
--- a/gprMax/hash_cmds_file.py
+++ b/gprMax/hash_cmds_file.py
@@ -219,12 +219,12 @@ def check_cmd_names(processedlines, checkessential=True):
'#waveform', '#voltage_source',
'#hertzian_dipole', '#magnetic_dipole',
'#transmission_line', '#rx', '#rx_array',
- '#snapshot', '#include_file']}
+ '#snapshot', '#include_file', '#material_range', '#material_list']}
# Geometry object building commands that there can be multiple instances
# of in a model - these will be lists within the dictionary
geometrycmds = ['#geometry_objects_read', '#edge', '#plate', '#triangle',
- '#box', '#sphere', '#cylinder', '#cylindrical_sector',
+ '#box', '#sphere', '#ellipsoid', '#cone', '#cylinder', '#cylindrical_sector',
'#fractal_box', '#add_surface_roughness',
'#add_surface_water', '#add_grass']
# List to store all geometry object commands in order from input file
diff --git a/gprMax/hash_cmds_geometry.py b/gprMax/hash_cmds_geometry.py
index 3f6f0926..dbe8f9c8 100644
--- a/gprMax/hash_cmds_geometry.py
+++ b/gprMax/hash_cmds_geometry.py
@@ -24,12 +24,14 @@ from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness
from .cmds_geometry.add_surface_water import AddSurfaceWater
from .cmds_geometry.box import Box
from .cmds_geometry.cylinder import Cylinder
+from .cmds_geometry.cone import Cone
from .cmds_geometry.cylindrical_sector import CylindricalSector
from .cmds_geometry.edge import Edge
from .cmds_geometry.fractal_box import FractalBox
from .cmds_geometry.plate import Plate
from .cmds_geometry.sphere import Sphere
from .cmds_geometry.triangle import Triangle
+from .cmds_geometry.ellipsoid import Ellipsoid
from .utilities.utilities import round_value
logger = logging.getLogger(__name__)
@@ -217,6 +219,40 @@ def process_geometrycmds(geometry):
scene_objects.append(cylinder)
+
+ elif tmp[0] == '#cone:':
+ if len(tmp) < 10:
+ logger.exception("'" + ' '.join(tmp) + "'" +
+ ' requires at least nine parameters')
+ raise ValueError
+
+ p1 = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
+ p2 = (float(tmp[4]), float(tmp[5]), float(tmp[6]))
+ r1 = float(tmp[7])
+ r2 = float(tmp[8])
+
+ # Isotropic case with no user specified averaging
+ if len(tmp) == 10:
+ cone = Cone(p1=p1, p2=p2, r1=r1, r2=r2, material_id=tmp[9])
+
+ # Isotropic case with user specified averaging
+ elif len(tmp) == 11:
+ averaging = check_averaging(tmp[10].lower())
+ cone = Cone(p1=p1, p2=p2, r1=r1, r2=r2, material_id=tmp[9],
+ averaging=averaging)
+
+ # Uniaxial anisotropic case
+ elif len(tmp) == 12:
+ cone = Cone(p1=p1, p2=p2, r1=r1, r2=r2, material_ids=tmp[9:])
+
+ else:
+ logger.exception("'" + ' '.join(tmp) + "'" +
+ ' too many parameters have been given')
+ raise ValueError
+
+ scene_objects.append(cone)
+
+
elif tmp[0] == '#cylindrical_sector:':
if len(tmp) < 10:
logger.exception("'" + ' '.join(tmp) + "'" +
@@ -295,6 +331,41 @@ def process_geometrycmds(geometry):
scene_objects.append(sphere)
+
+ elif tmp[0] == '#ellipsoid:':
+ if len(tmp) < 8:
+ logger.exception("'" + ' '.join(tmp) + "'" +
+ ' requires at least seven parameters')
+ raise ValueError
+
+ p1 = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
+ xr = float(tmp[4])
+ yr = float(tmp[5])
+ zr = float(tmp[6])
+
+ # Isotropic case with no user specified averaging
+ if len(tmp) == 8:
+ ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, material_id=tmp[7])
+
+ # Isotropic case with user specified averaging
+ elif len(tmp) == 9:
+ averaging = check_averaging(tmp[8].lower())
+ ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, material_id=tmp[7],
+ averaging=averaging)
+
+ # Uniaxial anisotropic case
+ elif len(tmp) == 8:
+ ellipsoid = Ellipsoid(p1=p1, xr=xr, yr=yr, zr=zr, material_id=tmp[7:])
+
+ else:
+ logger.exception("'" + ' '.join(tmp) + "'" +
+ ' too many parameters have been given')
+ raise ValueError
+
+ scene_objects.append(ellipsoid)
+
+
+
elif tmp[0] == '#fractal_box:':
# Default is no dielectric smoothing for a fractal box
diff --git a/gprMax/hash_cmds_multiuse.py b/gprMax/hash_cmds_multiuse.py
index 160c9d06..a46431a4 100644
--- a/gprMax/hash_cmds_multiuse.py
+++ b/gprMax/hash_cmds_multiuse.py
@@ -22,7 +22,7 @@ from .cmds_multiuse import (AddDebyeDispersion, AddDrudeDispersion,
AddLorentzDispersion, GeometryObjectsWrite,
GeometryView, HertzianDipole, MagneticDipole,
Material, Rx, RxArray, Snapshot, SoilPeplinski,
- TransmissionLine, VoltageSource, Waveform)
+ TransmissionLine, VoltageSource, Waveform, MaterialRange, MaterialList)
logger = logging.getLogger(__name__)
@@ -342,4 +342,47 @@ def process_multicmds(multicmds):
gow = GeometryObjectsWrite(p1=p1, p2=p2, filename=tmp[6])
scene_objects.append(gow)
+ cmdname = '#material_range'
+ if multicmds[cmdname] is not None:
+ for cmdinstance in multicmds[cmdname]:
+ tmp = cmdinstance.split()
+
+ if len(tmp) != 9:
+ logger.exception("'" + cmdname + ': ' + ' '.join(tmp) + "'" +
+ ' requires at exactly nine parameters')
+ raise ValueError
+ material_range = MaterialRange(er_lower=float(tmp[0]),
+ er_upper=float(tmp[1]),
+ sigma_lower=float(tmp[2]),
+ sigma_upper=float(tmp[3]),
+ mr_lower=float(tmp[4]),
+ mr_upper=float(tmp[5]),
+ ro_lower=float(tmp[6]),
+ ro_upper=float(tmp[7]),
+ id=tmp[8])
+ scene_objects.append(material_range)
+
+
+ cmdname = '#material_list'
+ if multicmds[cmdname] is not None:
+ for cmdinstance in multicmds[cmdname]:
+ tmp = cmdinstance.split()
+
+ if len(tmp) < 2:
+ logger.exception("'" + cmdname + ': ' + ' '.join(tmp) + "'" +
+ ' requires at least 2 parameters')
+ raise ValueError
+
+ tokens = len(tmp)
+ lmats = []
+ for iter in range(0,tokens-1):
+ lmats.append(tmp[iter])
+
+
+
+ material_list = MaterialList(list_of_materials=lmats,
+ id=tmp[tokens-1])
+ scene_objects.append(material_list)
+
+
return scene_objects
diff --git a/gprMax/materials.py b/gprMax/materials.py
index 95e07ba1..6b29b30a 100644
--- a/gprMax/materials.py
+++ b/gprMax/materials.py
@@ -17,10 +17,14 @@
# along with gprMax. If not, see .
import numpy as np
+import logging
+
import gprMax.config as config
+logger = logging.getLogger(__name__)
+
class Material:
"""Super-class to describe generic, non-dispersive materials,
their properties and update coefficients.
@@ -309,9 +313,12 @@ class PeplinskiSoil:
self.rb = bulkdensity
self.rs = sandpartdensity
self.mu = watervolfraction
- self.startmaterialnum = 0
+ self.startmaterialnum = 0 #This is not used anymore and code that uses it can be removed
+ # store all of the material IDs in a list instead of storing only the first number of the material
+ # and assume that all must be sequentially numbered. This allows for more general mixing models
+ self.matID = []
- def calculate_debye_properties(self, nbins, G):
+ def calculate_properties(self, nbins, G):
"""Calculates the real and imaginery part of a Debye model for the soil
as well as a conductivity. It uses an approximation to a semi-empirical
model (http://dx.doi.org/10.1109/36.387598).
@@ -340,11 +347,16 @@ class PeplinskiSoil:
# sigf = -1.645 + 1.939 * self.rb - 2.25622 * self.S + 1.594 * self.C
# Generate a set of bins based on the given volumetric water fraction
- # values
- mubins = np.linspace(self.mu[0], self.mu[1], nbins)
+ # values. Changed to make sure mid points are contained completely within the ranges.
+ # The limiting values of the ranges are not included in this.
+
+ #mubins = np.linspace(self.mu[0], self.mu[1], nbins)
+ mubins = np.linspace(self.mu[0], self.mu[1], nbins+1)
# Generate a range of volumetric water fraction values the mid-point of
# each bin to make materials from
- mumaterials = mubins + (mubins[1] - mubins[0]) / 2
+ #mumaterials = mubins + (mubins[1] - mubins[0]) / 2
+ mumaterials = 0.5*(mubins[1:nbins+1] + mubins[0:nbins])
+
# Create an iterator
muiter = np.nditer(mumaterials, flags=['c_index'])
@@ -368,8 +380,9 @@ class PeplinskiSoil:
if muiter.index == 0:
if material:
self.startmaterialnum = material.numID
+ self.matID.append(material.numID)
else:
- self.startmaterialnum = len(G.materials)
+ self.startmaterialnum = len(G.materials)
if not material:
m = DispersiveMaterial(len(G.materials), requiredID)
m.type = 'debye'
@@ -382,8 +395,174 @@ class PeplinskiSoil:
m.deltaer.append(er - eri)
m.tau.append(DispersiveMaterial.watertau)
G.materials.append(m)
+ self.matID.append(m.numID)
muiter.iternext()
+
+
+
+class RangeMaterial:
+ """Material defined with a given range of parameters to be used for
+ factal spatial disttibutions
+ """
+
+ def __init__(self, ID, er_range, sigma_range, mu_range, ro_range):
+ """
+ Args:
+ ID: string for name of the material.
+ er_range: tuple of floats for relative permittivity range of the material.
+ sigma_range: tuple of floats for electric conductivity range of the material.
+ mu_range: tuple of floats for magnetic permeability of material.
+ ro_range: tuple of floats for magnetic loss range of material.
+ """
+
+ self.ID = ID
+ self.er = er_range
+ self.sig = sigma_range
+ self.mu = mu_range
+ self.ro = ro_range
+ self.startmaterialnum = 0 #This is not really needed anymore and code that uses it can be removed.
+ # store all of the material IDs in a list instead of storing only the first number of the material
+ # and assume that all must be sequentially numbered. This allows for more general mixing models
+ self.matID = []
+
+
+ def calculate_properties(self, nbins, G):
+ """Calculates the properties of the materials.
+
+ Args:
+ nbins: int for number of bins to use to create the different materials.
+ G: FDTDGrid class describing a grid in a model.
+ """
+
+ # Generate a set of relative permittivity bins based on the given range
+ erbins = np.linspace(self.er[0], self.er[1], nbins+1)
+ # Generate a range of relative permittivity values the mid-point of
+ # each bin to make materials from
+ #ermaterials = erbins + np.abs((erbins[1] - erbins[0])) / 2
+ ermaterials = 0.5*(erbins[1:nbins+1] + erbins[0:nbins])
+
+ # Generate a set of conductivity bins based on the given range
+ sigmabins = np.linspace(self.sig[0], self.sig[1], nbins+1)
+ # Generate a range of conductivity values the mid-point of
+ # each bin to make materials from
+ #sigmamaterials = sigmabins + (sigmabins[1] - sigmabins[0]) / 2
+ sigmamaterials = 0.5*(sigmabins[1:nbins+1] + sigmabins[0:nbins])
+
+ # Generate a set of magnetic permeability bins based on the given range
+ mubins = np.linspace(self.mu[0], self.mu[1], nbins+1)
+ # Generate a range of magnetic permeability values the mid-point of
+ # each bin to make materials from
+ #mumaterials = mubins + np.abs((mubins[1] - mubins[0])) / 2
+ mumaterials = 0.5*(mubins[1:nbins+1] + mubins[0:nbins])
+
+ # Generate a set of magnetic loss bins based on the given range
+ robins = np.linspace(self.ro[0], self.ro[1], nbins+1)
+ # Generate a range of magnetic loss values the mid-point of
+ # each bin to make materials from
+ #romaterials = robins + np.abs((robins[1] - robins[0])) / 2
+ romaterials = 0.5*(robins[1:nbins+1] + robins[0:nbins])
+
+
+ # Iterate over the bins
+ for iter in np.arange(0,nbins):
+
+ # Relative permittivity
+ er = ermaterials[iter]
+
+ # Effective conductivity
+ se = sigmamaterials[iter]
+
+ # magnetic permeability
+ mr = mumaterials[iter]
+
+ # magnetic loss
+ sm = romaterials[iter]
+
+ # Check to see if the material already exists before creating a new one
+ requiredID = '|{:.4f}+{:.4f}+{:.4f}+{:.4f}|'.format(float(er),float(se),float(mr),float(sm))
+ material = next((x for x in G.materials if x.ID == requiredID), None)
+ if iter == 0:
+ if material:
+ self.startmaterialnum = material.numID
+ self.matID.append(material.numID)
+ else:
+ self.startmaterialnum = len(G.materials)
+ if not material:
+ m = Material(len(G.materials), requiredID)
+ m.type = ''
+ m.averagable = True
+ m.er = er
+ m.se = se
+ m.mr = mr
+ m.sm = sm
+ G.materials.append(m)
+ self.matID.append(m.numID)
+
+
+
+class ListMaterial:
+ """A list of predefined materials to be used for
+ factal spatial disttibutions. This command does not create new materials but collects them to be used in a
+ stochastic distribution by a fractal box.
+ """
+
+ def __init__(self, ID, listofmaterials):
+ """
+ Args:
+ ID: string for name of the material.
+ listofmaterials: A list of material IDs.
+
+ """
+
+ self.ID = ID
+ self.mat = listofmaterials
+ self.startmaterialnum = 0 #This is not really needed anymore
+ # store all of the material IDs in a list instead of storing only the first number of the material
+ # and assume that all must be sequentially numbered. This allows for more general mixing models
+ # this is important here as this model assumes predefined materials.
+ self.matID = []
+
+
+ def calculate_properties(self, nbins, G):
+ """Calculates the properties of the materials. No Debye is used but name kept the same as used in other
+ class that needs Debye
+
+ Args:
+ nbins: int for number of bins to use to create the different materials.
+ G: FDTDGrid class describing a grid in a model.
+ """
+
+
+ # Iterate over the bins
+ for iter in np.arange(0,nbins):
+
+ # Check to see if the material already exists before creating a new one
+ #requiredID = '|{:}_in_{:}|'.format((self.mat[iter]),(self.ID))
+ requiredID = self.mat[iter]
+ material = next((x for x in G.materials if x.ID == requiredID), None)
+ self.matID.append(material.numID)
+
+ #if iter == 0:
+ # if material:
+ # self.startmaterialnum = material.numID
+ # else:
+ # self.startmaterialnum = len(G.materials)
+
+ #if not material:
+ # temp = next((x for x in G.materials if x.ID == self.mat[iter]), None)
+ # m = copy.deepcopy(temp) #This needs to import copy in order to work
+ # m.ID = requiredID
+ # m.numID = len(G.materials)
+ # G.materials.append(m)
+
+
+ if not material:
+ logger.exception(self.__str__() + f' material(s) {material} do not exist')
+ raise ValueError
+
+
+
def create_built_in_materials(G):