New object building primitives and two new mixing models for fractal box, bug fixes

Introduced object primitives #ellipsoid and #cone. The #ellipsoid works with regular semi-axes
along x, y, and z. The #cone works in a similar way to the cylinder and can be truncated with circular faces of different radii. If one radius is zero then it is a normal cone.
Introduced two new mixing models for using with a fractal box, #material_range and #material_list
The #material_range allows materials to be constructed using ranges of permittivity, conductivity, magnetic permeability and magnetic loss. Ranges can be increasing or decreasing as needed.
The #material_list creates a list of predefined materials that can be used with a fractal box. The order of the listed materials reflects the order of the bins they are allocated to in a fractal box.
The fractal box and associated mixing models have been updated to use a list of material IDs and do not rely on the assumption that are generated sequentially. This gives greater flexibility and
adds little computation burden on the creation of the fractal volume.

Bug fixes for mixing anisotropic volumetric materials where "mr" was used instead of "ms".
Bug fixes on objects that use integer coordinates for centres but where passed floats instead.

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

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

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)

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)

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 <http://www.gnu.org/licenses/>.
+
+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}.")

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)

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)

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 <http://www.gnu.org/licenses/>.
+
+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}.")

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')

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, 

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)

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)

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.
 

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,

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

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
 

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

gprMax/materials.py

@@ -17,10 +17,14 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 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.
@@ -308,9 +312,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).
@@ -339,11 +346,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'])
@@ -367,6 +379,7 @@ class PeplinskiSoil:
             if muiter.index == 0:
                 if material:
                     self.startmaterialnum = material.numID
+                    self.matID.append(material.numID)                    
                 else:
                     self.startmaterialnum = len(G.materials)                
             if not material:
@@ -381,10 +394,176 @@ 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):
     """Creates pre-defined (built-in) materials.