Update to fractal behaviour. Fixes int/float bug & D now used directly.

gprMax/fractals.py

@@ -53,9 +53,7 @@ class FractalSurface(object):
         self.ny = yf - ys
         self.nz = zf - zs
         self.seed = seed
-        self.dimension = dimension
-        # Constant related to fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
-        self.b = -(2 * self.dimension - 7) / 2
+        self.dimension = dimension # Fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
         self.weighting = np.array([1, 1], dtype=np.float64)
         self.fractalrange = (0, 0)
         self.filldepth = 0
@@ -90,7 +88,7 @@ class FractalSurface(object):
         A = fftpack.fftshift(A)
 
         # Generate fractal
-        generate_fractal2D(surfacedims[0], surfacedims[1], G.nthreads, self.b, self.weighting, v1, A, self.fractalsurface)
+        generate_fractal2D(surfacedims[0], surfacedims[1], G.nthreads, self.dimension, self.weighting, v1, A, self.fractalsurface)
 
         # Shift the zero frequency component to start of the array
         self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
@@ -178,7 +176,7 @@ class FractalVolume(object):
         A = fftpack.fftshift(A)
 
         # Generate fractal
-        generate_fractal3D(self.nx, self.ny, self.nz, G.nthreads, self.b, self.weighting, v1, A, self.fractalvolume)
+        generate_fractal3D(self.nx, self.ny, self.nz, G.nthreads, self.dimension, self.weighting, v1, A, self.fractalvolume)
 
         # Shift the zero frequency component to the start of the array
         self.fractalvolume = fftpack.ifftshift(self.fractalvolume)

gprMax/fractals_generate_ext.pyx

@@ -21,13 +21,13 @@ cimport numpy as np
 from cython.parallel import prange
 
 
-cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, ::1] A, np.complex128_t[:, ::1] fractalsurface):
+cpdef void generate_fractal2D(int nx, int ny, int nthreads, float D, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, ::1] A, np.complex128_t[:, ::1] fractalsurface):
     """This function generates a fractal surface for a 2D array.
         
     Args:
         nx, ny (int): Fractal surface size in cells
         nthreads (int): Number of threads to use
-        b (int): Constant related to fractal dimension
+        D (float): Fractal dimension
         weighting (memoryview): Access to weighting vector
         v1 (memoryview): Access to positional vector at centre of array, scaled by weighting
         A (memoryview): Access to array containing random numbers (to be convolved with fractal function)
@@ -46,20 +46,20 @@ cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[
                 # Calulate norm of v2 - v1
                 rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2)**(1/2)
 
-                B = rr**b
+                B = rr**D
                 if B == 0:
                     B = 0.9
 
                 fractalsurface[i, j] = A[i, j] / B
 
 
-cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, :, ::1] A, np.complex128_t[:, :, ::1] fractalvolume):
+cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, float D, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, :, ::1] A, np.complex128_t[:, :, ::1] fractalvolume):
     """This function generates a fractal volume for a 3D array.
 
     Args:
         nx, ny, nz (int): Fractal volume size in cells
         nthreads (int): Number of threads to use
-        b (int): Constant related to fractal dimension
+        D (float): Fractal dimension
         weighting (memoryview): Access to weighting vector
         v1 (memoryview): Access to positional vector at centre of array, scaled by weighting
         A (memoryview): Access to array containing random numbers (to be convolved with fractal function)
@@ -79,7 +79,7 @@ cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.fl
 
                 # Calulate norm of v2 - v1
                 rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2 + (v2z - v1[2])**2)**(1/2)
-                B = rr**b
+                B = rr**D
                 if B == 0:
                     B = 0.9