From f787ef4ba544de03858a621a1a03742de9af0aee Mon Sep 17 00:00:00 2001
From: Craig Warren <craig_warren@me.com>
Date: Mon, 4 Apr 2016 18:04:11 +0100
Subject: [PATCH] Updated dispersion check/warning to handle user-defined
 waveforms. Assumes max frequency in simulation is -60dB from frequency with
 max power.

---
 gprMax/gprMax.py |  2 +-
 gprMax/grid.py   | 78 ++++++++++++++++++++++++++++++++----------------
 2 files changed, 53 insertions(+), 27 deletions(-)

diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py
index a9ccd460..8752156b 100644
--- a/gprMax/gprMax.py
+++ b/gprMax/gprMax.py
@@ -330,7 +330,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
             print('{:3}\t{:12}\tepsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/m; '.format(material.numID, material.ID, material.er, material.se, material.mr, material.sm) + tmp + dielectricsmoothing)
 
     # Check to see if numerical dispersion might be a problem
-    if dispersion_check(G.waveforms, G.materials, G.dx, G.dy, G.dz):
+    if dispersion_check(G):
         print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present.')
     
     # Write files for any geometry views
diff --git a/gprMax/grid.py b/gprMax/grid.py
index aac79308..6a82f712 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -17,6 +17,7 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import numpy as np
+import matplotlib.pyplot as plt
 
 from gprMax.constants import c, floattype, complextype
 from gprMax.materials import Material
@@ -91,43 +92,68 @@ class FDTDGrid:
         self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 * Material.maxpoles), dtype=complextype)
 
 
-def dispersion_check(waveforms, materials, dx, dy, dz):
+def dispersion_check(G):
     """Check for potential numerical dispersion. Is the smallest wavelength present in the simulation discretised by at least a factor of 10
         
     Args:
-        waveforms (list): Waveforms present in the model.
-        materials (list): Materials present in the model.
-        dx, dy, dz (float): Spatial discretisation of the model.
+        G (class): Grid class instance - holds essential parameters describing the model.
     
     Returns:
         (boolean): Potential numerical dispersion
     """
     
     # Find maximum frequency
-    freqs = [waveform.freq for waveform in waveforms]
-    if not freqs:
-        return False
-    maxfreq = max(freqs)
-    
-    # Find minimum wavelength
-    ers = [material.er for material in materials]
-    miner = max(ers)
+    maxfreqs = []
+    for waveform in G.waveforms:
+        
+        # User-defined waveform
+        if waveform.uservalues is not None:
+            waveformvalues = waveform.uservalues
+        
+        # Built-in waveform
+        else:
+            time = np.linspace(0, 1, G.iterations)
+            time *= (G.iterations * G.dt)
+            waveformvalues = np.zeros(len(time))
+            timeiter = np.nditer(time, flags=['c_index'])
 
-    # Minimum velocity
-    minvelocity = c / np.sqrt(miner)
+            while not timeiter.finished:
+                waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt)
+                timeiter.iternext()
+
+        # Calculate magnitude of frequency spectra of waveform
+        power = 20 * np.log10(np.abs(np.fft.fft(waveformvalues))**2)
+        freqs = np.fft.fftfreq(power.size, d=G.dt)
+
+        # Shift powers so that frequency with maximum power is at zero decibels
+        power -= np.amax(power)
+        
+        # Set maximum frequency to -60dB from maximum power
+        freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
+        maxfreqs.append(freqs[freq])
+
+    if maxfreqs:
+        maxfreq = max(maxfreqs)
+        
+        # Find minimum wavelength
+        ers = [material.er for material in G.materials]
+        miner = max(ers)
+
+        # Minimum velocity
+        minvelocity = c / np.sqrt(miner)
+
+        # Minimum number of spatial steps to resolve smallest wavelength
+        resolution = 10
+        
+        # Minimum wavelength
+        minwavelength = minvelocity / maxfreq
+        
+        # Test for numerical dispersion
+        if max((G.dx, G.dy, G.dz)) > (minwavelength / resolution):
+            return True
+        else:
+            return False
 
-    # Minimum number of spatial steps to resolve smallest wavelength
-    resolution = 10
-    
-    # Multiplier for centre frequency (2-3 times the centre frequency of a waveform can be present in the simulation)
-    maxfreq *= 3
-    
-    # Minimum wavelength
-    minwavelength = minvelocity / maxfreq
-    
-    # Test for numerical dispersion
-    if max((dx, dy, dz)) > (minwavelength / resolution):
-        return True
     else:
         return False