diff --git a/CHANGELOG.md b/CHANGELOG.md
index d35804e2..ea86361d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 This is a breaking release. It corrects the surface-EMG pipeline against Farina, Mesin, Martina & Merletti (2004). The problems were found through feedback from @joaohbbittar; thank you. Both the default route and `execution=` are affected.
 
+It also ports three fixes that never reached `main` (Gfluctdv noise, fibre conduction velocities, large fibre layouts) and corrects the units of the intramuscular SFAP.
+
 ### Breaking changes
 
 This release changes default outputs; regenerate anything calibrated against earlier versions.
@@ -16,7 +18,10 @@ This release changes default outputs; regenerate anything calibrated against ear
 - sEMG outputs change in amplitude, shape, duration and spectrum. Bipolar outputs from 0.5.0–0.13.0 are invalid. Durations, mean/median frequencies and conduction-velocity estimates from those versions are off by about 2× and must be regenerated.
 - Absolute sEMG amplitudes change by one to two orders of magnitude (amplitude no longer depends on conduction velocity). Anything calibrated in µV or mV against earlier outputs must be regenerated. To keep the earlier resistive (measured) skin layer, pass `skin_conductivity__S_per_m=4.88e-4 * pq.S / pq.m` to `Muscle`. To keep the earlier layers, pass `skin_thickness__mm=1.29 * pq.mm, fat_thickness__mm=0.3 * pq.mm`.
 - Pickled `SurfaceEMG` objects from earlier versions lack the new fibre properties; create them again.
-- The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed.
+- The `execution=` MUAP disk-cache key changes, so templates cached by earlier versions are recomputed (sEMG and iEMG).
+- Pools that call `insert_Gfluctdv()` now produce noisy discharge: spike trains change, and descending-drive parameters fitted while the noise was silent should be refitted.
+- iEMG outputs change through the fibre conduction velocities (mean ~2.0 → ~3.45 m/s at the default diameter). Muscles with more than ~313 000 fibres get a different, now circular, fibre layout.
+- iEMG amplitudes, shapes and spectra change through the SFAP unit fixes (typically 20–50× larger at a fine wire). Anything tuned against earlier iEMG amplitudes, including `add_noise` floors given in absolute units, must be redone. Values are in mV.
 - iEMG outputs change through the narrower end-plate spreads. Pass `IntramuscularEMG(..., nmj_coefficients="reference_code")` to reproduce earlier end-plates.
 
 ### Fixed
@@ -34,6 +39,15 @@ This release changes default outputs; regenerate anything calibrated against ear
 - **Default skin and fat thickness (FDI).** Skin 1.29 → 1.2 mm (hand dorsum median, Li et al. 2018) and fat 0.3 → 1.3 mm, so the total depth to muscle is 2.5 mm, the low end of the FDI depth measured by Nishikawa et al. (2024); the fat value is derived from that total, not measured. Surface amplitudes from the default `Muscle` fall by about 40 %.
 - **Failed motor units raise.** On the default `SurfaceEMG` route an exception while simulating one motor unit was logged and replaced by an all-zero MUAP, which passes for a valid, silent unit. The exception is now re-raised with the unit and array index; the `ExecutionPolicy` route already raised.
 - **Angular sampling overflow raises.** With `sampling_points_in_theta_domain` ≥ 256 the scaled modified Bessel functions leave the float64 range at the highest angular orders (I_n underflows to 0, K_n overflows). The volume conductor replaced the overflowing entries with zero, which gave a singular system on the CPU and NaN MUAPs on the GPU. It now raises a `ValueError` that names the setting.
+- **Gfluctdv membrane noise was silent.** `AlphaMN.insert_Gfluctdv()` inserted the Destexhe (2001) fluctuating-conductance mechanism but left its GLOBAL multipliers `multex_Gfluctdv` / `multin_Gfluctdv` at the `.mod` default of 0, so it added no noise. They are now set to 1 on insertion. Any pool that calls `insert_Gfluctdv()` now has conductance noise and more variable discharge; this also makes the `gfluctdv` search in `examples/02_finetune/01_optimize_dd_for_target_firing_rate.py` effective.
+- **Fibre conduction velocities.** `Muscle` mixed units in the velocity–diameter line (5.5 (m/s)/mm with diameters in mm), giving ~1.1 m/s, which the 2 m/s clip then raised: every fibre had exactly 2.0 m/s. Fibre velocities now follow 0.95 + 0.05·d[µm] m/s (Nandedkar & Stålberg 1983; ~3.45 m/s at 50 µm), the same line `SurfaceFiberDispersion` uses; noise (SD 0.2 m/s) and clip (2–6 m/s) are unchanged. Intramuscular EMG uses these velocities; surface EMG uses the mean velocity and is unaffected.
+- **Intramuscular SFAP units and polarity.** Four errors in the iEMG source and volume conductor (`bioelectric.py`, `MotorUnitSim.calc_sfaps`) are corrected on both the default and `execution=` routes (`runtime/density.py`, `runtime/iemg.py`); a single-fibre SFAP now matches an independent SI line-source solution within 3 % at 0.1–1 mm (`tests/test_iemg_sfap_units.py`).
+  - Muscle conductivities defaulted to 63 / 330 "S/m" (1000× the tissue values of 0.063 / 0.33 S/m), making potentials 1000× too small. Defaults are now 0.063 / 0.33 S/m.
+  - The Rosenfalck IAP was evaluated at twice the distance and differentiated with the chain rule applied twice: an IAP half the published length and 8× the current. It now uses λ = 1 mm, as the sEMG model does.
+  - The wave travelling towards the shorter fibre half was subtracted instead of added, so it had inverted polarity; the two halves cancelled near the end-plate and in the far field.
+  - The sum over the fibre omitted the spatial step, so amplitudes scaled with `1/spatial_resolution__mm`.
+  Together the bugs nearly cancelled at the default 0.01 mm step for a nearby electrode. On an FDI cylinder (bipolar fine wire, 2 mm IED, 15 units) MUAPs go from 0.03–1.8 mV to 0.9–37 mV peak-to-peak. These are point electrodes in an infinite homogeneous medium, which overestimates the near field (no electrode area, no near-field correction); see Breaking changes.
+- **Large circular fibre layouts.** Fibre centres came from a 399 499-point table covering a square, of which only ~313 000 points lie inside the circle. Beyond that the outer ring was filled only at the diagonals, and beyond 399 499 the fibre count was silently capped. Layouts the table covers are unchanged (FDI defaults ≈ 54 000 fibres); larger counts now use a scrambled Halton sequence in the disc, seeded from MyoGen's generator, with a warning.
 
 ### Added
 
diff --git a/myogen/simulator/core/emg/fiber_simulation.py b/myogen/simulator/core/emg/fiber_simulation.py
index b7b6e7e2..be665c10 100644
--- a/myogen/simulator/core/emg/fiber_simulation.py
+++ b/myogen/simulator/core/emg/fiber_simulation.py
@@ -55,8 +55,8 @@ def compute_intramuscular_kernel(
     z_grid: np.ndarray,
     electrode_z: np.ndarray,
     r: float,
-    sigma_r: float = 63.0,
-    sigma_z: float = 330.0,
+    sigma_r: float = 0.063,
+    sigma_z: float = 0.33,
 ) -> np.ndarray:
     """
     Compute the intramuscular volume conductor kernel using the Green's function.
@@ -72,9 +72,9 @@ def compute_intramuscular_kernel(
         Electrode z-positions in mm, shape (n_electrodes,).
     r : float
         Radial distance from fiber to electrode in mm.
-    sigma_r : float, default=63.0
+    sigma_r : float, default=0.063
         Radial conductivity in S/m.
-    sigma_z : float, default=330.0
+    sigma_z : float, default=0.33
         Longitudinal conductivity in S/m.
 
     Returns
diff --git a/myogen/simulator/core/emg/intramuscular/bioelectric.py b/myogen/simulator/core/emg/intramuscular/bioelectric.py
index 8a86e0e3..3325a8ce 100644
--- a/myogen/simulator/core/emg/intramuscular/bioelectric.py
+++ b/myogen/simulator/core/emg/intramuscular/bioelectric.py
@@ -4,14 +4,14 @@
 This module contains the core bioelectric modeling functions for simulating
 single fiber action potentials (SFAPs) and motor unit action potentials (MUAPs)
 in intramuscular EMG. The functions implement the volume conductor models
-from Farina et al. 2004 and the transmembrane current models from
-Rosenfalck 1969.
+from Farina et al. 2004 [farina2004] and the transmembrane current models from
+Rosenfalck 1969 [rosenfalck1969].
 
 References
 ----------
-[1] Farina, D., Merletti, R., 2001. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Transactions on Biomedical Engineering 48, 637–646. <br>
-[2] Rosenfalck, P., 1969. Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiologica Scandinavica Supplementum 321, 1–168. <br>
-[3] Nandedkar, S.D., Stålberg, E., 1983. Simulation of single muscle fibre action potentials. Medical & Biological Engineering & Computing 21, 158–165.
+[1] Farina, D., Merletti, R., 2001. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Transactions on Biomedical Engineering 48, 637–646. [farinamerletti2001] <br>
+[2] Rosenfalck, P., 1969. Intra- and extracellular potential fields of active nerve and muscle fibres. Acta Physiologica Scandinavica Supplementum 321, 1–168. [rosenfalck1969] <br>
+[3] Nandedkar, S.D., Stålberg, E., 1983. Simulation of single muscle fibre action potentials. Medical & Biological Engineering & Computing 21, 158–165. [nandedkar1983]
 """
 
 import numpy as np
@@ -22,7 +22,7 @@ def get_tm_current(z: np.ndarray, D1: float = 96.0, D2: float = -90.0) -> np.nda
     Calculate transmembrane current using Rosenfalck's model.
 
     This function implements the transmembrane current model from:
-    P. Rosenfalck "Intra and extracellular fields of active nerve and muscle fibers" (1969)
+    P. Rosenfalck "Intra and extracellular fields of active nerve and muscle fibers" (1969) [rosenfalck1969]
 
     Parameters
     ----------
@@ -45,7 +45,7 @@ def get_tm_current(z: np.ndarray, D1: float = 96.0, D2: float = -90.0) -> np.nda
 
 def get_tm_current_dz(z: np.ndarray, D1: float = 96.0, xp=np) -> np.ndarray:
     """
-    Calculate first derivative of transmembrane current (Rosenfalck model).
+    Calculate first derivative of transmembrane current (Rosenfalck model) [rosenfalck1969].
 
     This is the spatial derivative of the transmembrane current model used
     for action potential propagation simulation.
@@ -73,7 +73,7 @@ def get_tm_current_dz(z: np.ndarray, D1: float = 96.0, xp=np) -> np.ndarray:
 
 def get_tm_current_ddz(z: np.ndarray, D1: float = 96.0) -> np.ndarray:
     """
-    Calculate second derivative of transmembrane current (Rosenfalck model).
+    Calculate second derivative of transmembrane current (Rosenfalck model) [rosenfalck1969].
 
     Parameters
     ----------
@@ -98,8 +98,8 @@ def get_elementary_current_response(
     z: np.ndarray,
     z_electrode: float,
     r: np.ndarray,
-    sigma_r: float = 63.0,  # S/m
-    sigma_z: float = 330.0,  # S/m
+    sigma_r: float = 0.063,  # S/m
+    sigma_z: float = 0.33,  # S/m
     xp=np,
 ) -> np.ndarray:
     """
@@ -107,7 +107,7 @@ def get_elementary_current_response(
 
     This function calculates the potential response at electrode location
     due to a unit current source at different positions along the muscle fiber.
-    Based on Nandedkar & Stålberg 1983.
+    Based on Nandedkar & Stålberg 1983 [nandedkar1983].
 
     Parameters
     ----------
@@ -117,22 +117,23 @@ def get_elementary_current_response(
         Electrode position along z-axis in mm
     r : np.ndarray
         Radial distance from fiber to electrode in mm
-    sigma_r : float, default=63.0
-        Radial conductivity in S/m (from Andreassen & Rosenfalck 1980)
-    sigma_z : float, default=330.0
-        Longitudinal conductivity in S/m (from Andreassen & Rosenfalck 1980)
+    sigma_r : float, default=0.063
+        Radial conductivity in S/m (from Andreassen & Rosenfalck 1980) [andreassenrosenfalck1980]
+    sigma_z : float, default=0.33
+        Longitudinal conductivity in S/m (from Andreassen & Rosenfalck 1980) [andreassenrosenfalck1980]
     xp : module, default=np
         Array backend (numpy or cupy)
 
     Returns
     -------
     np.ndarray
-        Elementary current response (transfer function)
+        Transfer function in Ω: potential (mV) per mA of point current. Multiply the current per unit length from ``get_current_density``
+        (mA/mm) by this and by the spatial step dz (mm), and sum over z, to get the SFAP in mV.
     """
-    # ---- FIXED UNIT CONVERSIONS ----
-    # Convert conductivities from S/m to S/mm to match spatial units (mm)
-    sigma_r_S_per_mm = sigma_r / 1000.0  # CORRECTED: convert S/m → S/mm
-    sigma_z_S_per_mm = sigma_z / 1000.0  # CORRECTED: convert S/m → S/mm
+    # Unit system of the iEMG kernel: mm, mV, mA. 1 S/m = 1e-3 S/mm. The defaults were
+    # 63 / 330 "S/m" until 0.13.0 (1000x the tissue values), which made SFAPs 1000x too small.
+    sigma_r_S_per_mm = sigma_r / 1000.0
+    sigma_z_S_per_mm = sigma_z / 1000.0
 
     # Normalize inputs to computation backend (prevents numpy/cupy mixing)
     z = xp.asarray(z)
@@ -228,7 +229,7 @@ def get_current_density(
 ):
     """
     Model the individual action potential (IAP) or single fiber action potential (SFAP) in space and time.
-    Translated from Farina & Merletti (2001) and Nandedkar & Stålberg (1998).
+    Translated from Farina & Merletti (2001) [farinamerletti2001] and Nandedkar & Stålberg (1998).
 
     Parameters
     ----------
@@ -276,14 +277,17 @@ def tendon_terminator(z_inline, L_inline):
 
     # Compute psi (transmembrane current derivative)
     if L1 >= L2:
-        psi = -4 * get_tm_current_dz(-2 * (Z - zi - v * T), xp=xp)
+        # psi = dV_m/dz for the wave travelling towards +z, zeta = v t - (z - zi) in mm (Rosenfalck
+        # length scale 1 mm, as in the sEMG model). Until 0.13.0 this read -4 * V'(-2 * arg): the chain
+        # rule counted twice and the IAP compressed to half its length (8x the current).
+        psi = -get_tm_current_dz(-(Z - zi - v * T), xp=xp)
         longest_wave = xp.diff(psi, axis=0) / dz
         longest_wave *= tendon_terminator(Z[:-1, :] - zi - L1 / 2, L1)
         # Explicit bool→float64 cast required: CuPy does not support
         # implicit multiplication of bool arrays with float arrays.
         longest_wave *= ((Z[:-1, :] - zi) / v > 0).astype(xp.float64)
     else:
-        psi = 4 * get_tm_current_dz(-2 * (-Z + zi - v * T), xp=xp)
+        psi = get_tm_current_dz(-(-Z + zi - v * T), xp=xp)
         longest_wave = xp.diff(psi, axis=0) / dz
         longest_wave *= tendon_terminator(Z[:-1, :] - zi + L2 / 2, L2)
         longest_wave *= ((-Z[:-1, :] + zi) / v > 0).astype(xp.float64)  # bool→float64
@@ -296,10 +300,13 @@ def tendon_terminator(z_inline, L_inline):
 
     if L1 >= L2:
         shortest_wave *= tendon_terminator(Z[:-1, :] - zi + L2 / 2, L2)
-        iap = longest_wave - shortest_wave
+        # The mirrored array is already d2V/dz2 of the opposite wave (a second derivative is
+        # unchanged by reflection), so the two waves add. Until 0.13.0 they were subtracted,
+        # which inverted the wave travelling away from the longer half.
+        iap = longest_wave + shortest_wave
     else:
         shortest_wave *= tendon_terminator(Z[:-1, :] - zi - L1 / 2, L1)
-        iap = shortest_wave - longest_wave
+        iap = shortest_wave + longest_wave
 
     # Suppress endplate density if required
     if suppress_endplate_density:
@@ -315,7 +322,8 @@ def tendon_terminator(z_inline, L_inline):
     # Compute cross-sectional area in mm²
     area_mm2 = xp.pi * (d / 2) ** 2  # CORRECTED: removed extra /4
 
-    # Scale current density by intracellular conductivity and fiber cross-section area
+    # Scale current density by intracellular conductivity and fiber cross-section area:
+    # transmembrane current per unit length, S/mm * mm^2 * mV/mm^2 = mA/mm
     iap *= sigma_i * area_mm2
 
     return iap
@@ -329,7 +337,7 @@ def get_current_density_fast(
     L1: float,
     L2: float,
     v: float,
-    d: float = 55e-6,
+    d: float = 55e-3,
     suppress_endplate_density: bool = True,
     endplate_width: float = 0.5,
 ) -> np.ndarray:
@@ -355,7 +363,7 @@ def get_current_density_fast(
         Length from endplate to negative tendon in mm
     v : float
         Conduction velocity in mm/s
-    d : float, default=55e-6
+    d : float, default=55e-3
         Fiber diameter in mm
     suppress_endplate_density : bool, default=True
         Whether to suppress endplate region
diff --git a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
index 468ae26a..938d0d46 100644
--- a/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
+++ b/myogen/simulator/core/emg/intramuscular/intramuscular_emg.py
@@ -314,7 +314,8 @@ def _calculate_muaps(self, n_jobs: int = -2, verbose: bool = True) -> INTRAMUSCU
         Returns
         -------
         INTRAMUSCULAR_MUAP__Block
-            Intramuscular MUAP shapes for all motor units stored in a neo.Block.
+            Intramuscular MUAP shapes for all motor units stored in a neo.Block, in mV
+            (tagged ``pq.dimensionless`` for backward compatibility).
         """
         if not self._motor_units:
             raise ValueError("Must call _initialize_motor_units() first")
@@ -516,8 +517,8 @@ def simulate_intramuscular_emg(
         Returns
         -------
         INTRAMUSCULAR_EMG__Block
-            Intramuscular EMG signals for the electrode array stored in a neo.Block.
-            Results are stored in the `intramuscular_emg__Block` property after execution.
+            Intramuscular EMG signals for the electrode array stored in a neo.Block, in mV
+            (tagged ``pq.dimensionless`` for backward compatibility). Results are stored in the `intramuscular_emg__Block` property after execution.
 
         Raises
         ------
diff --git a/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py b/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
index 86ab3af8..ded55f12 100644
--- a/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
+++ b/myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
@@ -321,7 +321,8 @@ def calc_sfaps(
                         ))
 
                     ecr_all = cp.column_stack(ecr_list)
-                    sfaps_dev[:, :, fiber_idx] = current_density.T @ ecr_all
+                    # mA/mm x Ω x dz (mm) summed over z: SFAP in mV
+                    sfaps_dev[:, :, fiber_idx] = (current_density.T @ ecr_all) * dz
 
             cp.cuda.Device(0).synchronize()
             self.sfaps = cp.asnumpy(sfaps_dev)
@@ -378,7 +379,8 @@ def calc_sfaps(
                     ))
 
                 ecr_all = np.column_stack(ecr_list)
-                self.sfaps[:, :, fiber_idx] = current_density.T @ ecr_all
+                # mA/mm x Ω x dz (mm) summed over z: SFAP in mV
+                self.sfaps[:, :, fiber_idx] = (current_density.T @ ecr_all) * dz
 
         self.shift_sfaps(dt)
 
diff --git a/myogen/simulator/core/emg/runtime/density.py b/myogen/simulator/core/emg/runtime/density.py
index c844f40f..cafac44e 100644
--- a/myogen/simulator/core/emg/runtime/density.py
+++ b/myogen/simulator/core/emg/runtime/density.py
@@ -22,13 +22,17 @@ def _density(t, z, active, mirror, short_mask, end_mask, zi, v, dz, positive, sc
         for col in range(columns):
             value = 0.0
             if needed:
+                # zeta = v t - |z - zi| in mm, Rosenfalck length scale 1 mm, as
+                # in bioelectric.get_current_density (until 0.13.0 both read
+                # -2 * (...) and scaled by 4: half the IAP length, 8x the current).
                 if positive:
-                    x = -2 * (z[row] - zi - v * t[col])
+                    x = -(z[row] - zi - v * t[col])
                 else:
-                    x = -2 * (-z[row] + zi - v * t[col])
+                    x = -(-z[row] + zi - v * t[col])
                 if x > 0:
                     value = 96.0 * (3 * x**2 - x**3) * math.exp(-x)
-                value *= -4 if positive else 4
+                if positive:
+                    value = -value
             wave[row, col] = value
     # Ascending rows retain the unmodified next endpoint until it is consumed.
     for row in range(rows):
@@ -42,7 +46,9 @@ def _density(t, z, active, mirror, short_mask, end_mask, zi, v, dz, positive, sc
         source = mirror[row]
         for col in range(columns):
             short = wave[source, col] if source >= 0 and short_mask[row] else 0.0
-            value = wave[row, col] - short if positive else short - wave[row, col]
+            # The mirrored wave is already d2V/dz2 of the opposite wave, so the
+            # two add (they were subtracted until 0.13.0).
+            value = wave[row, col] + short
             result[row, col] = (value if end_mask[row] else 0.0) * scale
     return result
 
diff --git a/myogen/simulator/core/emg/runtime/iemg.py b/myogen/simulator/core/emg/runtime/iemg.py
index 2b2fb92b..d3d6864f 100644
--- a/myogen/simulator/core/emg/runtime/iemg.py
+++ b/myogen/simulator/core/emg/runtime/iemg.py
@@ -10,7 +10,7 @@
 from .session import ExecutionSession
 from .cache import DiskCache, fingerprint
 
-IEMG_CACHE_NAMESPACE = "intramuscular-muap-density-v2"
+IEMG_CACHE_NAMESPACE = "intramuscular-muap-density-v3"
 
 
 def stream_muap(mu, dt, dz, electrodes, session, *, batch_size=8, label="fibres"):
@@ -56,7 +56,8 @@ def compute(start, stop):
                 if radius < minimum_radius:
                     radius = minimum_radius
                 responses.append(get_elementary_current_response(z, electrode[2], radius))
-            values[local] = current.T @ np.column_stack(responses)
+            # Line-source sum over the fibre: the spatial step is part of the integral.
+            values[local] = (current.T @ np.column_stack(responses)) * dz
             delay = float(delays[index] / dt)
             integer = int(np.floor(delay))
             fraction = delay - integer
diff --git a/myogen/simulator/core/emg/surface/fiber_dispersion.py b/myogen/simulator/core/emg/surface/fiber_dispersion.py
index ec6d2f7e..6caeb20a 100644
--- a/myogen/simulator/core/emg/surface/fiber_dispersion.py
+++ b/myogen/simulator/core/emg/surface/fiber_dispersion.py
@@ -17,6 +17,7 @@
 
 import numpy as np
 
+from ...muscle.muscle import CV_INTERCEPT__M_S, CV_PER_DIAMETER__M_S_PER_MM
 from ..nerve_branches import NMJ_SCHEDULE_COEFFICIENTS, nerve_delays__ms, nmj_branch_schedule, two_layer_nmj
 
 NMJ_MODELS = ("uniform_mu", "two_layer", "two_layer_intra")
@@ -97,8 +98,8 @@ class SurfaceFiberDispersion:
     fiber_cv_threshold_slope: float = 0.0
     fiber_cv_within_mu_cov: float = 0.0
     diameter_scales_source: bool = False
-    cv_per_diameter__m_s_per_mm: float = 50.0
-    cv_intercept__m_s: float = 0.95
+    cv_per_diameter__m_s_per_mm: float = CV_PER_DIAMETER__M_S_PER_MM
+    cv_intercept__m_s: float = CV_INTERCEPT__M_S
     tendon_end_band__mm: float = 0.0
     tendon_end_correlation__mm: float = 0.0
     tendon_end_jitter__mm: float = 0.0
diff --git a/myogen/simulator/core/muscle/muscle.py b/myogen/simulator/core/muscle/muscle.py
index 50e432b9..dc1092f6 100644
--- a/myogen/simulator/core/muscle/muscle.py
+++ b/myogen/simulator/core/muscle/muscle.py
@@ -1,4 +1,5 @@
 import inspect
+import warnings
 from pathlib import Path
 from typing import Optional
 
@@ -7,7 +8,7 @@
 import skfmm
 from joblib import Parallel, delayed
 from scipy.integrate import dblquad
-from scipy.stats import chi2, multivariate_normal
+from scipy.stats import chi2, multivariate_normal, qmc
 from sklearn.neighbors import NearestNeighbors
 from tqdm import tqdm
 import quantities as pq
@@ -24,6 +25,21 @@
 )
 from myogen.utils.decorators import beartowertype
 
+# Fibre conduction velocity against diameter: CV = 0.95 m/s + 0.05 m/s per um
+# (Nandedkar & Stalberg 1983; 3.7 m/s at 55 um). Shared with SurfaceFiberDispersion.
+CV_INTERCEPT__M_S = 0.95
+CV_PER_DIAMETER__M_S_PER_MM = 50.0
+
+
+def _halton_disc(n: int, rng: np.random.Generator) -> np.ndarray:
+    """``n`` quasi-uniform points in the unit disc: a scrambled Halton sequence in the square, kept inside the disc."""
+    sampler = qmc.Halton(d=2, scramble=True, rng=rng)
+    points = np.empty((0, 2))
+    while len(points) < n:
+        p = sampler.random(int((n - len(points)) * 4 / np.pi * 1.05) + 64) * 2 - 1
+        points = np.vstack([points, p[np.hypot(p[:, 0], p[:, 1]) <= 1]])
+    return points[:n]
+
 
 def _perform_fast_marching(speed_map: np.ndarray, seed_points: np.ndarray) -> np.ndarray:
     """
@@ -347,11 +363,9 @@ def _generate_fiber_properties(self) -> None:
             self._muscle_fiber_diameters__mm, 20e-3 * pq.mm, 80e-3 * pq.mm
         )
 
-        # Generate conduction velocities based on fiber diameter
-        # CV = k * diameter + c, where k ≈ 4.5-6.0 (m/s)/mm, c ≈ 0.5-1.0 m/s
-        # Based on Hakansson (1956) and later studies
-        k = 5.5 * (pq.m / pq.s) / pq.mm  # (m/s)/mm
-        c = 0.8 * pq.m / pq.s  # m/s
+        # Generate conduction velocities based on fiber diameter (diameters are in mm)
+        k = CV_PER_DIAMETER__M_S_PER_MM * (pq.m / pq.s) / pq.mm  # (m/s)/mm
+        c = CV_INTERCEPT__M_S * pq.m / pq.s  # m/s
 
         # Add some biological variability
         cv_base = k * self._muscle_fiber_diameters__mm + c
@@ -485,19 +499,22 @@ def generate_muscle_fiber_centers(self, verbose: bool = True) -> None:
         )
         sorted_indices = np.argsort(dists)
 
-        if len(sorted_indices) >= self._number_of_muscle_fibers + 1:
+        # The table covers the square [-1.25, 1.25]^2, so only points inside the inscribed
+        # disc (~313 k) form a full circle; beyond that the outer ring holds only the corners.
+        n = self._number_of_muscle_fibers
+        if n + 1 <= len(sorted_indices) and dists[sorted_indices[n]].magnitude <= 1.25:
             self._muscle_fiber_centers__mm = (
-                self._muscle_fiber_centers__mm[sorted_indices[: self._number_of_muscle_fibers], :]
-                / dists[sorted_indices[self._number_of_muscle_fibers]]
+                self._muscle_fiber_centers__mm[sorted_indices[:n], :]
+                / dists[sorted_indices[n]]
                 * self._radius__mm
             )
         else:
-            self._muscle_fiber_centers__mm = (
-                self._muscle_fiber_centers__mm[sorted_indices, :]
-                / dists[sorted_indices[-1]]
-                * self._radius__mm
+            warnings.warn(
+                f"{n} fibres exceed the circular part of voronoi_pi1e5.csv; placing them on a "
+                "scrambled Halton sequence in the disc instead.",
+                stacklevel=2,
             )
-            self._number_of_muscle_fibers = len(self._muscle_fiber_centers__mm)
+            self._muscle_fiber_centers__mm = _halton_disc(n, get_random_generator()) * self._radius__mm
 
         # Remove fibers inside the bone boundary
         # Fibers should only exist in the muscle tissue, not in the bone core
diff --git a/myogen/simulator/neuron/cells.py b/myogen/simulator/neuron/cells.py
index f9629522..8847d377 100644
--- a/myogen/simulator/neuron/cells.py
+++ b/myogen/simulator/neuron/cells.py
@@ -856,6 +856,9 @@ def insert_Gfluctdv(self):
         These values produce realistic subthreshold membrane fluctuations
         (~2-3 mV RMS) typical of in vivo motor neuron recordings.
         """
+        # Enable the global multipliers (default 0 in the .mod file)
+        h.multex_Gfluctdv = 1
+        h.multin_Gfluctdv = 1
         for d in self.dend:
             d.insert("Gfluctdv")
             # Excitatory fluctuations
diff --git a/tests/test_gfluctdv_noise.py b/tests/test_gfluctdv_noise.py
new file mode 100644
index 00000000..45cd2e27
--- /dev/null
+++ b/tests/test_gfluctdv_noise.py
@@ -0,0 +1,31 @@
+"""Gfluctdv (Destexhe 2001 fluctuating conductance) must produce noise once inserted.
+
+Gfluctdv.mod scales both conductances by GLOBAL multipliers ``multex``/``multin``
+that default to 0, so inserting the mechanism without enabling them gives zero
+conductance and no membrane noise.
+"""
+
+import numpy as np
+from neuron import h
+
+from myogen import load_nmodl_mechanisms
+from myogen.simulator.neuron.populations import AlphaMN__Pool
+
+
+def test_insert_gfluctdv_produces_conductance_noise():
+    load_nmodl_mechanisms()
+    h.multex_Gfluctdv = 0  # the .mod default; a previous test must not mask it
+    h.multin_Gfluctdv = 0
+    pool = AlphaMN__Pool(recruitment_thresholds__array=np.array([0.1, 1.0]))
+    cell = pool[0]
+    cell.insert_Gfluctdv()
+
+    d = cell.dend[0]
+    g_e = h.Vector().record(d(0.5)._ref_g_e_Gfluctdv)
+    h.dt = 0.025
+    h.finitialize(-70)
+    while h.t < 50:
+        h.fadvance()
+
+    g = np.asarray(g_e)
+    assert np.std(g) > 0.1 * d.std_e_Gfluctdv, f"g_e std {np.std(g):.3g} S/cm2: Gfluctdv inserted but silent"
diff --git a/tests/test_iemg_execution_route.py b/tests/test_iemg_execution_route.py
new file mode 100644
index 00000000..284d896c
--- /dev/null
+++ b/tests/test_iemg_execution_route.py
@@ -0,0 +1,55 @@
+"""The ``execution=`` iEMG route reproduces the default route's physics.
+
+``runtime/density.py`` is a numba port of ``bioelectric.get_current_density``
+(the reference); it must follow every change to the reference.
+"""
+import numpy as np
+import pytest
+import quantities as pq
+
+from myogen import set_random_seed
+from myogen.simulator import (
+    ExecutionPolicy,
+    IntramuscularElectrodeArray,
+    IntramuscularEMG,
+    Muscle,
+    RecruitmentThresholds,
+)
+from myogen.simulator.core.emg.intramuscular import bioelectric
+from myogen.simulator.core.emg.runtime import density
+
+
+@pytest.mark.parametrize("L1, L2", [(12.0, 10.0), (10.0, 12.0), (11.0, 11.0)])
+@pytest.mark.parametrize("suppress", [True, False])
+def test_kernel_matches_the_reference_current_density(L1, L2, suppress):
+    t = np.arange(0.0, 20.0 + 1e-9, 0.1)[:, None]
+    z = np.arange(-15.0, 15.0 + 1e-9, 0.05)[:, None]
+    args = (t, z, 0.0, L1, L2, 4.0, 0.055)
+    expected = bioelectric.get_current_density(*args, suppress_endplate_density=suppress)
+    actual = density.current_density(*args, suppress_endplate_density=suppress)
+    np.testing.assert_allclose(actual, expected, rtol=1e-10, atol=1e-12 * np.max(np.abs(expected)))
+
+
+@pytest.mark.usefixtures("ample_session_memory")
+def test_execution_route_muaps_match_the_default_route(monkeypatch):
+    from myogen.simulator.core.emg.intramuscular import intramuscular_emg
+
+    monkeypatch.setattr(intramuscular_emg, "HAS_CUPY", False)
+    set_random_seed(5)
+    thresholds, _ = RecruitmentThresholds(N=4, recruitment_range__ratio=30)
+    muscle = Muscle(recruitment_thresholds=thresholds, radius__mm=2.0 * pq.mm,
+                    radius_bone__mm=0.5 * pq.mm, fiber_density__fibers_per_mm2=50 * pq.mm**-2,
+                    grid_resolution=32, autorun=True)
+    electrode = IntramuscularElectrodeArray(num_electrodes=2, inter_electrode_distance__mm=0.5 * pq.mm,
+                                            position__mm=(0.0 * pq.mm, 0.0 * pq.mm, 15.0 * pq.mm))
+    muaps = []
+    for execution in (None, ExecutionPolicy()):
+        iemg = IntramuscularEMG(muscle_model=muscle, electrode_array=electrode, MUs_to_simulate=[0, 3],
+                                sampling_frequency__Hz=10240.0 * pq.Hz)
+        set_random_seed(11)
+        muaps.append(iemg.simulate_muaps(n_jobs=1, verbose=False, execution=execution))
+    for reference, optimized in zip(muaps[0].segments, muaps[1].segments):
+        x = np.asarray(reference.analogsignals[0].magnitude)
+        y = np.asarray(optimized.analogsignals[0].magnitude)
+        assert x.shape == y.shape
+        np.testing.assert_allclose(y, x, rtol=1e-7, atol=1e-9 * np.max(np.abs(x)))
diff --git a/tests/test_iemg_sfap_units.py b/tests/test_iemg_sfap_units.py
new file mode 100644
index 00000000..ee2718ae
--- /dev/null
+++ b/tests/test_iemg_sfap_units.py
@@ -0,0 +1,84 @@
+"""Physical units of the intramuscular SFAP (production path through ``MotorUnitSim.calc_sfaps``).
+
+Reference: an independent line source in an anisotropic infinite medium, in SI units.
+Rosenfalck V_m(zeta) = 96 zeta^3 exp(-zeta) - 90 mV with zeta in mm [rosenfalck1969];
+I_m = sigma_i pi (d/2)^2 d2V_m/dz2 with sigma_i = 1.01 S/m [rosenfalck1969];
+phi = 1/(4 pi sigma_r) sum I_m dz / sqrt(K r^2 + (z - z_e)^2), K = sigma_z / sigma_r,
+sigma_r = 0.063 S/m and sigma_z = 0.33 S/m [andreassenrosenfalck1980] (as cited in ``bioelectric.py``).
+Two waves leave the end-plate; the comparison window ends before they reach the fibre ends.
+"""
+import numpy as np
+import pytest
+
+from myogen.simulator.core.emg.intramuscular.motor_unit_sim import MotorUnitSim
+
+D_MM = 0.055
+V_MM_S = 4000.0
+L_MM = 100.0
+ZI_MM = 50.0
+ZE_MM = 60.0
+DT_S = 1 / 50_000
+T_END_S = 40.0 / V_MM_S  # waves still 10 mm from the fibre ends
+
+
+def _production_sfap(r_mm, dz_mm, use_gpu=False, ze_mm=ZE_MM, zi_mm=ZI_MM, t_end_s=T_END_S):
+    sim = MotorUnitSim(
+        muscle_fiber_centers__mm=np.array([[0.0, 0.0]]),
+        muscle_length__mm=L_MM,
+        muscle_fiber_diameters__mm=np.array([D_MM]),
+        muscle_fiber_conduction_velocity__mm_per_s=np.array([V_MM_S]),
+        nominal_center__mm=np.array([0.0, 0.0]),
+    )
+    sim._neuromuscular_z_coordinates__mm = np.array([[zi_mm]])
+    sim.nerve_paths = np.zeros((1, 2))
+    sim.calc_sfaps(0, DT_S, dz_mm, np.array([[r_mm, 0.0, ze_mm]]), verbose=False, use_gpu=use_gpu)
+    sfap = np.asarray(sim.sfaps)[:, 0, 0]
+    t = np.arange(sfap.size) * DT_S
+    return sfap[t < t_end_s]
+
+
+def _si_reference(r_mm, dz_mm=0.01):
+    z = np.arange(0.0, L_MM + dz_mm, dz_mm)
+    t = np.arange(0.0, T_END_S, DT_S)
+    zeta = V_MM_S * t[None, :] - np.abs(z[:, None] - ZI_MM)
+    zp = np.clip(zeta, 0.0, None)
+    vm = 96.0 * zp**3 * np.exp(-zp) - 90.0
+    d2v = np.gradient(np.gradient(vm, dz_mm, axis=0), dz_mm, axis=0) * 1e3  # mV/mm^2 -> V/m^2
+    d2v[np.abs(z - ZI_MM) < 0.5] = 0.0  # end-plate zone suppressed, as in get_current_density
+    i_m = 1.01 * np.pi * (D_MM * 1e-3 / 2) ** 2 * d2v  # A/m
+    g = 1.0 / (4 * np.pi * 0.063 * np.sqrt(0.33 / 0.063 * (r_mm * 1e-3) ** 2 + ((z - ZE_MM) * 1e-3) ** 2))
+    return (i_m * g[:, None]).sum(0) * dz_mm * 1e-3 * 1e3  # V -> mV
+
+
+@pytest.mark.parametrize("r_mm", [0.1, 0.5, 1.0])
+def test_sfap_matches_si_line_source(r_mm):
+    sim = np.ptp(_production_sfap(r_mm, 0.01))
+    ref = np.ptp(_si_reference(r_mm))
+    assert sim == pytest.approx(ref, rel=0.03), f"r {r_mm} mm: production {sim:.4g} mV, SI {ref:.4g} mV"
+
+
+def test_sfap_mirror_symmetric_about_endplate():
+    """Both waves leave the end-plate with the same polarity: V_m(v t - |z - zi|) is symmetric in z."""
+    plus = _production_sfap(0.3, 0.01, ze_mm=ZI_MM + 10.0)
+    minus = _production_sfap(0.3, 0.01, ze_mm=ZI_MM - 10.0)
+    np.testing.assert_allclose(minus, plus, rtol=0, atol=0.02 * np.ptp(plus))
+
+
+def test_sfap_reflection_with_unequal_halves():
+    """End-plate at 40 mm or at 60 mm (L1 < L2 and L1 > L2): reflected geometries record the same SFAP."""
+    a = _production_sfap(0.3, 0.01, ze_mm=50.0, zi_mm=40.0, t_end_s=np.inf)
+    b = _production_sfap(0.3, 0.01, ze_mm=50.0, zi_mm=60.0, t_end_s=np.inf)
+    np.testing.assert_allclose(a, b, rtol=0, atol=0.02 * np.ptp(a))
+
+
+def test_sfap_independent_of_spatial_step():
+    a = np.ptp(_production_sfap(0.3, 0.01))
+    b = np.ptp(_production_sfap(0.3, 0.02))
+    assert a == pytest.approx(b, rel=0.02)
+
+
+def test_gpu_sfap_matches_cpu():
+    pytest.importorskip("cupy")
+    cpu = _production_sfap(0.3, 0.01, use_gpu=False)
+    gpu = _production_sfap(0.3, 0.01, use_gpu=True)
+    np.testing.assert_allclose(gpu, cpu, rtol=1e-6, atol=1e-9 * np.ptp(cpu))
diff --git a/tests/test_muscle_fibre_properties.py b/tests/test_muscle_fibre_properties.py
new file mode 100644
index 00000000..7a34c2da
--- /dev/null
+++ b/tests/test_muscle_fibre_properties.py
@@ -0,0 +1,55 @@
+"""Per-fibre conduction velocities and the fibre layout of the circular ``Muscle``."""
+
+import numpy as np
+import pytest
+import quantities as pq
+
+from myogen import set_random_seed
+from myogen.simulator import Muscle, RecruitmentThresholds
+
+
+def _muscle(radius__mm, density=None):
+    thr, _ = RecruitmentThresholds(N=20, recruitment_range__ratio=50, mode="fuglevand")
+    kw = {} if density is None else {"fiber_density__fibers_per_mm2": density * pq.mm**-2}
+    return Muscle(recruitment_thresholds=thr, radius__mm=radius__mm * pq.mm, **kw)
+
+
+def test_fibre_cv_follows_the_diameter_line():
+    """CV = 0.95 m/s + 0.05 m/s per um of diameter (Nandedkar & Stalberg 1983): ~3.45 m/s at 50 um.
+
+    With diameters in mm, the old slope of 5.5 (m/s)/mm gave ~1.1 m/s and the 2 m/s clip
+    set every fibre to exactly 2.0 m/s.
+    """
+    set_random_seed(1)
+    m = _muscle(4)
+    m.generate_muscle_fiber_centers(verbose=False)
+    m._generate_fiber_properties()
+    d_um = np.asarray(m._muscle_fiber_diameters__mm.rescale(pq.mm)) * 1e3
+    cv = np.asarray(m._muscle_fiber_conduction_velocities__mm_per_s.rescale(pq.m / pq.s))
+
+    assert np.mean(np.isclose(cv, 2.0)) < 0.05
+    assert np.median(cv) == pytest.approx(0.95 + 0.05 * np.median(d_um), abs=0.05)
+    assert np.corrcoef(d_um, cv)[0, 1] > 0.9
+
+
+def _layout(radius__mm, density):
+    set_random_seed(1)
+    m = _muscle(radius__mm, density)
+    m.generate_muscle_fiber_centers(verbose=False)
+    return m, np.asarray(m._muscle_fiber_centers__mm.rescale(pq.mm))
+
+
+def test_large_fibre_counts_are_not_capped():
+    m, xy = _layout(20, 400)  # 502 655 fibres requested, more than the 399 499-point table
+    expected = int(np.rint(np.pi * 20**2 * 400))
+    assert len(xy) == m._number_of_muscle_fibers == expected
+
+
+def test_large_fibre_layouts_stay_circular():
+    """Fibres in the outer ring must be spread around the whole circle, not only the diagonals."""
+    _, xy = _layout(20, 300)  # 376 991 fibres: beyond the table's 313 k inscribed-disc points
+    r = np.hypot(*xy.T)
+    assert r.max() <= 20 + 1e-9
+    ring = xy[r > 0.97 * 20]
+    counts, _ = np.histogram(np.arctan2(ring[:, 1], ring[:, 0]), bins=16, range=(-np.pi, np.pi))
+    assert counts.min() > 0.5 * counts.mean()