Merge branch 'devel' of https://github.com/gprMax/gprMax into beta_1

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, 
+                m.sm = np.mean((materials[0].sm, materials[1].sm, 
-                                materials[2].mr), axis=0)
+                                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, 
+                m.sm = np.mean((materials[0].sm, materials[1].sm, 
-                                materials[2].mr), axis=0)
+                                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, 
+                    m.sm = np.mean((materials[0].sm, materials[1].sm, 
-                                    materials[2].mr), axis=0)
+                                    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, 
+                m.sm = np.mean((materials[0].sm, materials[1].sm, 
-                                materials[2].mr), axis=0)
+                                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, 
+                    m.sm = np.mean((materials[0].sm, materials[1].sm, 
-                                    materials[2].mr), axis=0)
+                                    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.
@@ -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
+        # values. Changed to make sure mid points are contained completely within the ranges. 
-        mubins = np.linspace(self.mu[0], self.mu[1], nbins)
+        # 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):