From 493db241d746b4d8e7a1c74f86e2f2ab91c9f2bb Mon Sep 17 00:00:00 2001
From: Nathan Mannall <n.mannall@epcc.ed.ac.uk>
Date: Tue, 15 Apr 2025 17:32:35 +0100
Subject: [PATCH] Remove shifting zero frequency from fractal surfaces

---
 gprMax/cython/fractals_generate.pyx              | 12 ++++++++++--
 gprMax/fractals.py                               |  9 +++++----
 gprMax/user_objects/cmds_geometry/fractal_box.py |  6 ------
 3 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/gprMax/cython/fractals_generate.pyx b/gprMax/cython/fractals_generate.pyx
index d7dd03c1..d00e7a37 100644
--- a/gprMax/cython/fractals_generate.pyx
+++ b/gprMax/cython/fractals_generate.pyx
@@ -26,6 +26,10 @@ from gprMax.config cimport float_or_double_complex
 cpdef void generate_fractal2D(
     int nx,
     int ny,
+    int ox,
+    int oy,
+    int gx,
+    int gy,
     int nthreads,
     float D,
     np.float64_t[:] weighting,
@@ -50,12 +54,16 @@ cpdef void generate_fractal2D(
 
     cdef Py_ssize_t i, j
     cdef double v2x, v2y, rr, B
+    cdef int sx, sy
+
+    sx = gx // 2
+    sy = gy // 2
 
     for i in prange(nx, nogil=True, schedule='static', num_threads=nthreads):
         for j in range(ny):
                 # Positional vector for current position
-                v2x = weighting[0] * i
-                v2y = weighting[1] * j
+                v2x = weighting[0] * ((i + ox + sx) % gx)
+                v2y = weighting[1] * ((j + oy + sy) % gy)
 
                 # Calulate norm of v2 - v1
                 rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2)**(1/2)
diff --git a/gprMax/fractals.py b/gprMax/fractals.py
index 5a79c7ce..8ca6492b 100644
--- a/gprMax/fractals.py
+++ b/gprMax/fractals.py
@@ -108,11 +108,13 @@ class FractalSurface:
 
         # 2D FFT
         A = fftpack.fftn(A)
-        # Shift the zero frequency component to the centre of the array
-        A = fftpack.fftshift(A)
 
         # Generate fractal
         generate_fractal2D(
+            surfacedims[0],
+            surfacedims[1],
+            0,
+            0,
             surfacedims[0],
             surfacedims[1],
             config.get_model_config().ompthreads,
@@ -122,8 +124,7 @@ class FractalSurface:
             A,
             self.fractalsurface,
         )
-        # Shift the zero frequency component to start of the array
-        self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
+
         # Set DC component of FFT to zero
         self.fractalsurface[0, 0] = 0
         # Take the real part (numerical errors can give rise to an imaginary part)
diff --git a/gprMax/user_objects/cmds_geometry/fractal_box.py b/gprMax/user_objects/cmds_geometry/fractal_box.py
index 1f020c91..4f718b1e 100644
--- a/gprMax/user_objects/cmds_geometry/fractal_box.py
+++ b/gprMax/user_objects/cmds_geometry/fractal_box.py
@@ -223,27 +223,21 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                 for surface in self.volume.fractalsurfaces:
                     if surface.surfaceID == "xminus":
                         if surface.fractalrange[0] < self.volume.xs:
-                            self.volume.nx += self.volume.xs - surface.fractalrange[0]
                             self.volume.xs = surface.fractalrange[0]
                     elif surface.surfaceID == "xplus":
                         if surface.fractalrange[1] > self.volume.xf:
-                            self.volume.nx += surface.fractalrange[1] - self.volume.xf
                             self.volume.xf = surface.fractalrange[1]
                     elif surface.surfaceID == "yminus":
                         if surface.fractalrange[0] < self.volume.ys:
-                            self.volume.ny += self.volume.ys - surface.fractalrange[0]
                             self.volume.ys = surface.fractalrange[0]
                     elif surface.surfaceID == "yplus":
                         if surface.fractalrange[1] > self.volume.yf:
-                            self.volume.ny += surface.fractalrange[1] - self.volume.yf
                             self.volume.yf = surface.fractalrange[1]
                     elif surface.surfaceID == "zminus":
                         if surface.fractalrange[0] < self.volume.zs:
-                            self.volume.nz += self.volume.zs - surface.fractalrange[0]
                             self.volume.zs = surface.fractalrange[0]
                     elif surface.surfaceID == "zplus":
                         if surface.fractalrange[1] > self.volume.zf:
-                            self.volume.nz += surface.fractalrange[1] - self.volume.zf
                             self.volume.zf = surface.fractalrange[1]
 
                 # If there is only 1 bin then a normal material is being used,