Constants import from scipy

testing/analytical_solutions.py

@@ -17,8 +17,9 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import numpy as np
+from scipy.constants import c
+from scipy.constants import epsilon_0 as e0
 
-import gprMax.config as config
 from gprMax.waveforms import Waveform
 
 
@@ -78,35 +79,35 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
     Ex_y = y
     Ex_z = z - 0.5 * dz
     Er_x = np.sqrt((Ex_x**2 + Ex_y**2 + Ex_z**2))
-    tau_Ex = Er_x / config.sim_config.em_consts["c"]
+    tau_Ex = Er_x / c
 
     # Coordinates of Rx for Ey FDTD component
     Ey_x = x
     Ey_y = y + 0.5 * dy
     Ey_z = z - 0.5 * dz
     Er_y = np.sqrt((Ey_x**2 + Ey_y**2 + Ey_z**2))
-    tau_Ey = Er_y / config.sim_config.em_consts["c"]
+    tau_Ey = Er_y / c
 
     # Coordinates of Rx for Ez FDTD component
     Ez_x = x
     Ez_y = y
     Ez_z = z
     Er_z = np.sqrt((Ez_x**2 + Ez_y**2 + Ez_z**2))
-    tau_Ez = Er_z / config.sim_config.em_consts["c"]
+    tau_Ez = Er_z / c
 
     # Coordinates of Rx for Hx FDTD component
     Hx_x = x
     Hx_y = y + 0.5 * dy
     Hx_z = z
     Hr_x = np.sqrt((Hx_x**2 + Hx_y**2 + Hx_z**2))
-    tau_Hx = Hr_x / config.sim_config.em_consts["c"]
+    tau_Hx = Hr_x / c
 
     # Coordinates of Rx for Hy FDTD component
     Hy_x = x + 0.5 * dx
     Hy_y = y
     Hy_z = z
     Hr_y = np.sqrt((Hy_x**2 + Hy_y**2 + Hy_z**2))
-    tau_Hy = Hr_y / config.sim_config.em_consts["c"]
+    tau_Hy = Hr_y / c
 
     # Initialise fields
     fields = np.zeros((iterations, 6))
@@ -134,7 +135,7 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
         fdot_Hy = wdot.calculate_value((timestep * dt) - tau_Hy, dt) * dl
 
         # Ex
-        fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_x**5)) * (
+        fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * e0 * Er_x**5)) * (
             3 * (fint_Ex + (tau_Ex * f_Ex)) + (tau_Ex**2 * fdot_Ex)
         )
 
@@ -145,12 +146,12 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
             tmp = 0
         fields[timestep, 1] = (
             tmp
-            * ((Ey_x * Ey_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_y**5))
+            * ((Ey_x * Ey_z) / (4 * np.pi * e0 * Er_y**5))
             * (3 * (fint_Ey + (tau_Ey * f_Ey)) + (tau_Ey**2 * fdot_Ey))
         )
 
         # Ez
-        fields[timestep, 2] = (1 / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_z**5)) * (
+        fields[timestep, 2] = (1 / (4 * np.pi * e0 * Er_z**5)) * (
             (2 * Ez_z**2 - (Ez_x**2 + Ez_y**2)) * (fint_Ez + (tau_Ez * f_Ez))
             - (Ez_x**2 + Ez_y**2) * tau_Ez**2 * fdot_Ez
         )