From b3426882c390f71427badbeaf54a5972c429531d Mon Sep 17 00:00:00 2001
From: Craig Warren <craig_warren@me.com>
Date: Tue, 11 Oct 2016 13:48:35 +0100
Subject: [PATCH] Corrected bug with adding dispersion analysis results to a
 dict.

---
 gprMax/grid.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/gprMax/grid.py b/gprMax/grid.py
index 18f64805..fa9b06ba 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -198,18 +198,24 @@ def dispersion_analysis(G):
                 print(Fore.RED + "\nWARNING: Duration of source waveform '{}' means it does not fit within specified time window and is therefore being truncated.".format(waveform.ID) + Style.RESET_ALL)
 
     if maxfreqs:
-        maxfreq = max(maxfreqs)
+        results['maxfreq'] = max(maxfreqs)
 
         # Find minimum wavelength (material with maximum permittivity)
-        ers = [x.er for x in G.materials if x.ID != 'pec']
+        ers = []
+        for x in G.materials:
+            if x.se != float('inf'):
+                er = x.er
+                if x.deltaer:
+                    er += max(x.deltaer)
+                ers.append(er)
         maxer = max(ers)
-        material = next(x for x in G.materials if x.er == maxer and x.ID != 'pec')
+        results['material'] = next(x for x in G.materials if x.er == maxer and x.ID != 'pec')
 
         # Minimum velocity
         minvelocity = c / np.sqrt(maxer)
 
         # Minimum wavelength
-        minwavelength = minvelocity / maxfreq
+        minwavelength = minvelocity / results['maxfreq']
         
         # Maximum spatial step
         delta = max(G.dx, G.dy, G.dz)
@@ -218,13 +224,13 @@ def dispersion_analysis(G):
         S = (c * G.dt) / delta
         
         # Grid sampling density
-        N = minwavelength / delta
+        results['N'] = minwavelength / delta
         
         # Numerical phase velocity
-        vp = np.pi / (N * np.arcsin((1 / S) * np.sin((np.pi * S) / N)))
+        vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])))
         
         # Physical phase velocity error (percentage)
-        deltavp = (((vp * c) - c) / c) * 100
+        results['deltavp'] = (((vp * c) - c) / c) * 100
 
     return results