diff --git a/share/binde_simplified/params_binde_simplified.lua b/share/binde_simplified/params_binde_simplified.lua
index c138e16a099ed1212c62ae8e5a1f681ba31b95e5..8423a089baf81045ab9f4b6b4a3f2f3683c8552b 100644
--- a/share/binde_simplified/params_binde_simplified.lua
+++ b/share/binde_simplified/params_binde_simplified.lua
@@ -1,3 +1,10 @@
+--[[
+    This example is a simplified version of a more complicated model
+    of an electrostatic discharge applied to a metallic box.
+    The discharge is modeled as a surface current source, but this is
+    not physically correct. This model is more a test for interface
+    BCs and sources than anything else.
+--]]
 sim.name = "binde_simplified"           -- simulation name
 sim.dt = 1e-11                          -- timestep size
 sim.timesteps = 10000                   -- num of iterations
@@ -8,6 +15,8 @@ sim.approx_order = 1                    -- approximation order
 postpro.silo_output_rate = 10           -- rate at which to write silo files
 postpro.cycle_print_rate = 10           -- console print rate
 
+sim.time_integrator = "leapfrog"        -- use leapfrog with dissipative materials
+
 -- ** Materials **
 -- Aluminum
 local alu = {2}
@@ -16,11 +25,6 @@ for i,v in ipairs(alu) do
     materials[v].epsilon = 1
     materials[v].mu = 1
     materials[v].sigma = 1
-
-    local coeff = sim.dt*materials[v].sigma/(materials[v].epsilon*const.eps0)
-    if math.abs(1-coeff) > 1 then
-        print("WARNING: sigma too high on domain " .. v .. ", this simulation is going to explode")
-    end
 end
 
 local air = {3}
@@ -38,14 +42,17 @@ for i,v in ipairs(Z_bnds) do
     bndconds[v].kind = "impedance"
 end
 
+-- Linear interpolation
 function interp(t, t0, t1, y0, y1)
     return (t-t0)*(y1-y0)/(t1-t0) + y0
 end
 
+-- Shape of the electrostatic discharge given as time-value pairs
 local current_shape = { {0, 0}, {9.5e-9, 0.03}, {1.1e-8, 1}, 
                         {1.4e-8, 0.45}, {1.6e-8, 0.45}, {2.7e-8, 0.6},
                         {9e-8, 0.01}, {1e-7, 0} }
 
+-- Function that evaluates the source
 function interface_source_7(tag, x, y, z, t)
     for ii = 1,(#current_shape-1) do
         local t0 = current_shape[ii][1]
@@ -60,15 +67,10 @@ function interface_source_7(tag, x, y, z, t)
     return 0, 0, 0
 end
 
+-- Surface current interface condition
 local discharge_surf = 7
 ifaceconds[discharge_surf] = {}
 ifaceconds[discharge_surf].kind = "surface_current"
 ifaceconds[discharge_surf].source = interface_source_7
 
-function on_timestep(ts)
-    --local time = ts*sim.dt
-    --local curr = interface_source_7(0,0,0,0,time)
-    --print("Time: " .. time .. ", current: " .. curr)
-end
-
 
diff --git a/share/modeconv/params_modeconv.lua b/share/modeconv/params_modeconv.lua
index 67b816ce7d2421c1096bbb031090bc47d4d97c57..20d3c8e99f1610a9bcb8cb7add0a1013108f56f2 100644
--- a/share/modeconv/params_modeconv.lua
+++ b/share/modeconv/params_modeconv.lua
@@ -1,3 +1,6 @@
+--[[
+    This example simulates a modal converter at 15.5 GHz.
+]]--
 sim.name = "modeconv"               -- simulation name
 sim.dt = 1e-14                      -- timestep size
 sim.timesteps = 50000               -- num of iterations
diff --git a/share/yagi/params_yagi.lua b/share/yagi/params_yagi.lua
index 4d73f80e198595d559bdbd73c33462ead7cac9dc..e3d17fdad4abcd9c9a422b5ef7ce4a50415e67d6 100644
--- a/share/yagi/params_yagi.lua
+++ b/share/yagi/params_yagi.lua
@@ -1,3 +1,6 @@
+--[[
+    Model of a 3-element yagi antenna receiving a plane wave at 300 MHz
+]]--
 sim.name = "yagi"              -- simulation name
 sim.dt = 1e-12                      -- timestep size
 sim.timesteps = 2000               -- num of iterations
@@ -15,12 +18,7 @@ for i,v in ipairs(alu) do
     materials[v] = {}
     materials[v].epsilon = 1
     materials[v].mu = 1
-    materials[v].sigma = 17
-
-    local coeff = sim.dt*materials[v].sigma/(materials[v].epsilon*const.eps0)
-    if math.abs(1-coeff) > 1 then
-        print("WARNING: sigma too high on domain " .. v .. ", this simulation is going to explode")
-    end
+    materials[v].sigma = 3.69e7 
 end
 
 local air = {5}
diff --git a/share/yagi_rad/params_yagi_rad.lua b/share/yagi_rad/params_yagi_rad.lua
new file mode 100644
index 0000000000000000000000000000000000000000..de966b78518e76ccfac3d90223b98860353b2359
--- /dev/null
+++ b/share/yagi_rad/params_yagi_rad.lua
@@ -0,0 +1,58 @@
+--[[
+    Model of a 3-element Yagi antenna radiating at 300 MHz.
+    The source is modeled as a volumetric current between the
+    two halves of the radiator.
+--]]
+sim.name = "yagi_rad"              -- simulation name
+sim.dt = 1e-13                      -- timestep size
+sim.timesteps = 50000             -- num of iterations
+sim.gmsh_model = "yagi_rad.geo"    -- gmsh model filename
+sim.use_gpu = 1                     -- 0: cpu, 1: gpu
+sim.approx_order = 2                -- approximation order
+--sim_geom_order = 1                -- geometric order, only 1 for now
+sim.time_integrator = "leapfrog"
+postpro.silo_output_rate = 500       -- rate at which to write silo files
+postpro.cycle_print_rate = 10       -- console print rate
+
+postpro["H"].silo_mode = "none"
+postpro["J"].silo_mode = "zonal"
+
+-- ** Materials **
+-- Aluminum
+local alu = {2, 3, 6, 8}
+for i,v in ipairs(alu) do
+    materials[v] = {}
+    materials[v].epsilon = 1
+    materials[v].mu = 1
+    materials[v].sigma = 3.69e7 
+end
+
+local air = {4, 7}
+for i,v in ipairs(air) do
+    materials[v] = {}
+    materials[v].epsilon = 1
+    materials[v].mu = 1
+    materials[v].sigma = 0
+end
+
+-- ** Boundary conditions **
+local Z_bnds = { 17 }
+for i,v in ipairs(Z_bnds) do
+    bndconds[v] = {}
+    bndconds[v].kind = "impedance"
+end
+
+local freq = 3e8
+function source(tag, x, y, z, t)
+    local Ex = 0
+    local Ey = 0
+    local Ez = math.sin(2*math.pi*t*freq)
+    return Ex, Ey, Ez
+end
+
+-- Apply the source to subdomain with tag 4
+sources[4] = source
+
+enable_boundary_sources(false)
+enable_interface_sources(false)
+
diff --git a/share/yagi_rad/yagi_rad.geo b/share/yagi_rad/yagi_rad.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ed83f19a4aae802bab70cf9fffe8431ceb458b5b
--- /dev/null
+++ b/share/yagi_rad/yagi_rad.geo
@@ -0,0 +1,30 @@
+SetFactory("OpenCASCADE");
+
+Mesh.MeshSizeFromCurvature = 5;
+Mesh.MeshSizeMax=0.4;
+
+spradius = 3;
+MHz = 300;
+
+ER = 0.005;
+
+R = 0.5*(300/MHz);
+A = 0.44*(300/MHz);
+D = 0.43*(300/MHz);
+sp_RA = 0.25*(300/MHz);
+sp_AD = 0.13*(300/MHz);
+
+src = A/20;
+Ax = (A-src)/2;
+
+Sphere(1) = {0, 0, 0, spradius, -Pi/2, Pi/2, 2*Pi};
+Cylinder(2) = {-sp_RA, 0, -R/2, 0, 0, R, ER, 2*Pi};
+Cylinder(3) = {0, 0, src/2, 0, 0, Ax, ER, 2*Pi};
+Cylinder(4) = {0, 0, -src/2, 0, 0, src, ER, 2*Pi};
+Cylinder(5) = {0, 0, src/2-Ax, 0, 0, Ax, ER, 2*Pi};
+Cylinder(6) = {sp_AD, 0, -D/2, 0, 0, D, ER, 2*Pi};
+Coherence;
+
+//MeshSize { Volume{2}  } = 0.01;
+//MeshSize {15,16,17,18,19,20,21,22} = 0.1;
+//+