Alipio soil model was wrong

Makefile

@@ -37,6 +37,9 @@ grid.o	:	src/grid.c
 linalg.o	:	src/linalg.c
 		$(CC) $(CFLAGS) -fPIC $(INCLUDE) -c src/linalg.c $(LINK)
 
+test	:	$(OBJECTS)
+		$(CC) $(CFLAGS) -fPIC $(INCLUDE) -Isrc/test -o test.exe src/test/testing.c $(OBJECTS) $(LINK)
+
 # Libraries
 dynamic_library	:	libhphem.so
 

examples/plot_results/plot_visacro57emc01.jl

@@ -69,19 +69,17 @@ end
 
 ## Step Voltage along a line in the +x direction
 begin
-    df_gpd = DataFrame(CSV.File("visacro57emc01_gpd.csv"))
     df_efield = DataFrame(CSV.File("visacro57emc01_efield.csv"))
-    t = sort(unique(df_gpd.t))
+    t = sort(unique(df_efield.t))
     nt = length(t)
-    ninj = (size(df_gpd)[2] - 3)  # number of injections
+    ninj = (size(df)[2] - 3) ÷ 2  # number of injections
     for i = 1:ninj
-        gpd = eval(Meta.parse("df_gpd.i$(i)"))
         ecx = eval(Meta.parse("df_efield.ecx$(i)"))
         encx = eval(Meta.parse("df_efield.ecx$(i)"))
         for yi in [0, 8]
             for k in time_steps
                 # integral of the E conservative field =======================
-                x = sort(unique(df.x))
+                x = sort(unique(df_efield.x))
                 dx = x[2] - x[1]
                 c = isapprox.(df_efield.t, t[k]) .& isapprox.(df_efield.y, yi)
                 ef = ecx[c] #+ encx[c]
@@ -97,7 +95,7 @@ begin
                          xlabel="x [m]", label="pot. dif.");
 
                 # integral of the E total field ==============================
-                x = sort(unique(df.x))
+                x = sort(unique(df_efield.x))
                 dx = x[2] - x[1]
                 c = isapprox.(df_efield.t, t[k]) .& isapprox.(df_efield.y, yi)
                 ef = ecx[c] + encx[c]

examples/visacro57emc01.c

@@ -55,9 +55,9 @@ run_case (double tmax, int nt, double Lmax, double* inj_t, const unsigned int NR
     double mur = 1.0;  // soil rel. magnetic permeability
     double sigma0 = 1.0 / 2000.0;  // soil conductivity in low frequency
     // parameters that I fitted:
-    double h_soil = 0.925 * pow(sigma0 * 1e3, -0.73);
-    double g_soil = 0.325;
-    double eps_ratio = 1.850;  // soil rel. permitivitty ratio
+    double h_soil = 2.1020 * pow(sigma0 * 1e3, -0.73);
+    double g_soil = 0.50890;
+    double eps_ratio = 3.6858;  // soil rel. permitivitty ratio
 
     // Laplace transform of the injection currents =============================
     // do it explicitly for performance reasons, see FFTW3 manual: http://fftw.org/fftw3_doc/
@@ -155,7 +155,6 @@ run_case (double tmax, int nt, double Lmax, double* inj_t, const unsigned int NR
             break;
         }
     }
-
     // define an array of points where to calculate ground scalar electric potential (GPD)
     // They form a grid of (1 x 1) [m^2] meshes
     double offset = 6.0;  // distance from groundig grid where to begin calculating GPD and fields
@@ -451,6 +450,7 @@ run_case (double tmax, int nt, double Lmax, double* inj_t, const unsigned int NR
                     }
                 }
             }
+
             if (i == 0) {
                 printf("Expected more time until completion of the frequency loop: %.2f min.\n",
                        (omp_get_wtime() - begin) * ns / 60.0 / omp_get_num_threads());

src/auxiliary.c

@@ -54,7 +54,7 @@ alipio_soil (_Complex double* sigma, _Complex double* epsr, double sigma0,
              _Complex double s, double h, double g, double eps_ratio)
 {
 
-    _Complex double f = s / TWO_PI;
+    _Complex double f = s / (TWO_PI * I);
     *sigma = (sigma0 + sigma0 * h * cpow(f/1e6, g));
     double t = tan(PI * g / 2) / (TWO_PI * EPS0 * pow(1e6, g));
     *epsr = eps_ratio + t * sigma0 * h * cpow(f, g - 1.0);