Centre a sphere of radius $r$ on the charge. With $\mathbf E=\frac{1}{4\pi\varepsilon_0}\frac{q}{r^2}\hat{\mathbf r}$ and $d\mathbf a=r^2\sin\theta\,d\theta\,d\varphi\,\hat{\mathbf r}$, the two factors of $r^2$ cancel:

The radius drops out, so the flux is identical through any enclosing surface (field lines are unbroken in charge-free space); a charge placed outside contributes nothing, since every line that enters the surface also leaves it. Hence $\Phi_E=q/\varepsilon_0$, independent of size and shape. (flux_through_sphere gives $\varepsilon_0\Phi_E=q=10^{-9}$ C at $R=0.1,0.5,2$ m, enclosed_charge reads $q$ back, and the charge at $(3,0,0)$ returns $\sim2\times10^{-15}\approx0$.)

Symmetry forces $\mathbf E=E(r)\hat{\mathbf r}$, so a concentric Gaussian sphere of radius $r$ gives $\oint\mathbf E\cdot d\mathbf a=E(r)\,4\pi r^2=Q_{\text{enc}}/\varepsilon_0$. Outside ($r\ge R$) the surface encloses all of $Q$:

the point-charge field. Inside ($r<R$) only the fraction $(r/R)^3$ of the uniform charge lies within, $Q_{\text{enc}}=Q(r/R)^3$:

At $r=R$ both branches give $kQ/R^2$, so $E$ is continuous. For $Q=2$ nC, $R=0.1$ m, field_magnitude returns $199.7$ V/m at $r=3R$ ($=kQ/9R^2$) and $898.8$ V/m at $r=R/2$ ($=kQ/2R^2$), matching the two formulas, and both branches meet at $kQ/R^2=1797.5$ V/m across $r=R$.

Cylindrical symmetry makes $\mathbf E=E(s)\hat{\mathbf s}$, radial and the same all along and around the axis. On a coaxial cylinder of radius $s$ and length $\ell$ the flat caps carry no flux ($\mathbf E\perp\hat{\mathbf n}$), and the curved wall (area $2\pi s\ell$) has $\mathbf E\parallel d\mathbf a$:

The length cancels, leaving

purely radial and falling as $1/s$. For $\lambda=2$ nC/m, field_magnitude of line_charge_field gives $359.5$ V/m at $s=0.1$ m and $179.8$ V/m at $s=0.2$ m — a clean factor-of-two $1/s$ drop — while the axial component $E_z=0$ off-axis, confirming the field is exactly radial.

By symmetry the field points straight away from the sheet with equal magnitude on each side, $\mathbf E=\pm E\,\hat{\mathbf n}$. A pillbox of face area $A$ straddling the sheet catches flux only through its two faces (the field grazes the side walls):

No $A$ and no distance survive, so the field is uniform. It reverses across the sheet, so the normal component jumps by

For $\sigma=1$ nC/m², plane_sheet_field returns the distance-independent $56.5$ V/m on either side and a jump of $112.9$ V/m $=\sigma/\varepsilon_0$ across $z=0$.

Write $Q_{\text{enc}}=\int_{\mathcal V}\rho\,d\tau$ and convert the flux with the divergence theorem (~MA-02):

Since this holds for every volume, the integrands match: $\nabla\cdot\mathbf E=\rho/\varepsilon_0$. On the P2 sphere, inside $\mathbf E=\tfrac{kQ}{R^3}\mathbf r$ (with $\mathbf r=(x,y,z)$), so

outside, the inverse-square $\mathbf E=kQ\hat{\mathbf r}/r^2$ has $\nabla\cdot\mathbf E=0$. The companion $\nabla\times\mathbf E=\mathbf 0$ completes the pair. Hence gauss_residual with $\rho=Q/\tfrac43\pi R^3$ returns $\sim1.4\times10^{-13}$ at the interior point — zero up to the finite-difference truncation, confirming $\nabla\cdot\mathbf E=\rho/\varepsilon_0$.

The Coulomb field is radial, $\mathbf E=\frac{1}{4\pi\varepsilon_0}\frac{q}{r^2}\hat{\mathbf r}$, so integrating in from infinity along a radial path ($d\mathbf l=dr\,\hat{\mathbf r}$):

Inverting with $\mathbf E=-\nabla V$ (here $V$ depends only on $r$):

recovering the EM-01 field — the integral and the gradient invert each other. Numerically potential_of_charge(q) gives $V(0.5)=35.9502\text{ V}=kq/r$, and field_from_potential returns $(102.944,\,68.629,\,-34.315)$ at $(0.3,0.2,-0.1)$, matching point_charge_field(q).

For any electrostatic field $\nabla\times\mathbf E=\mathbf 0$, so Stokes' theorem gives $\oint\mathbf E\cdot d\mathbf l=\int(\nabla\times\mathbf E)\cdot d\mathbf a=0$ around every closed loop. Split a loop into a forward path $\mathcal C_1$ and a return $\mathcal C_2$ between the same endpoints:

so $V(\mathbf b)-V(\mathbf a)=-\int_{\mathbf a}^{\mathbf b}\mathbf E\cdot d\mathbf l$ depends only on the endpoints. Shifting the reference point adds the same constant to $V(\mathbf p_1)$ and $V(\mathbf p_2)$, leaving the difference fixed. Two different references give the drop $12.0903$ V vs $12.0904$ V — equal up to the finite-step ($n=4000$) integration error — confirming path-independence.

The scalar sum of the two $q/\eta$ terms is

where $\eta_+$ is the distance to $+q$ at $(-a,0,0)$. Under $x\to-x$, $\eta_+\leftrightarrow\eta_-$, so $V\to-V$ — odd — and on the plane $x=0$ (where $\eta_+=\eta_-$) it vanishes identically. For $r\gg a$, expand $1/\eta_\pm\approx\frac1r\mp\frac{a x}{r^3}$:

a dipole of magnitude $2qa$ pointing from the $-q$ toward the $+q$ charge (the leading ~EM-05 term). On the $+x$ axis $V<0$ — you sit nearer $-q$ — matching potential_point_charges(...), which returns $0$ on the $x=0$ plane and $V(0.2)=-47.93\text{ V}=-V(-0.2)$.

Off the source $\rho=0$, so Poisson reduces to Laplace, $\nabla^2V=0$. For $V\propto1/r$ use the radial part of the spherical Laplacian, $\nabla^2 f=\frac{1}{r^2}\frac{d}{dr}\!\left(r^2\frac{df}{dr}\right)$:

So the point-charge potential is harmonic everywhere off the charge — it is $1/r$, not the field's $1/r^2$, that solves Laplace's equation. The code confirms it: poisson_residual(V, 0.0, point) returns $-6.4\times10^{-8}$, zero to finite-difference truncation error.

For the quadratic bowl $V=A(x^2+y^2+z^2)$ each second derivative is constant, $\partial_x^2V=\partial_y^2V=\partial_z^2V=2A$, so

Poisson's equation $\nabla^2V=-\rho/\varepsilon_0$ then forces a uniform source $\rho=-\varepsilon_0\nabla^2V=-6A\varepsilon_0$. Feeding this into poisson_residual(V, rho, point) gives $\approx-3\times10^{-9}$ (zero up to truncation), whereas the wrong rho=0.0 returns $6.0$: the residual $\nabla^2V+\rho/\varepsilon_0$ normalized by the curvature scale $|V|/L^2=A$ equals $6A/A=6$, exactly diagnosing the missing source.

On the plane $z=0$ both radicals reduce to the same value $\sqrt{x^2+y^2+d^2}$, so the bracket is $q/\sqrt{x^2+y^2+d^2}-q/\sqrt{x^2+y^2+d^2}=0$ and $V\equiv0$ there. Each term is a Coulomb potential $\propto1/r$, satisfying $\nabla^2(1/r)=0$ away from its source; the image sits at $(0,0,-d)$, outside the region $z>0$, so $\nabla^2V=0$ everywhere in the upper half-space except at the real charge. This candidate has the correct $1/r$ singularity at $q$, solves Laplace's equation, vanishes on $z=0$, and $\to0$ at infinity. By the first uniqueness theorem the boundary data fix the solution, so this is the potential for $z>0$. Hence image_potential_plane(q, d) returns $0$ at every $(x,y,0)$ and image_field_plane gives the matching $\mathbf E=-\nabla V$.

Differentiate the P1 potential in $z$; with $r_\mp^2=x^2+y^2+(z\mp d)^2$ the distances to the real ($-$) and image ($+$) charges,

Then $\sigma=-\varepsilon_0\,\partial V/\partial z|_{z=0}=-qd\big/\bigl[2\pi(x^2+y^2+d^2)^{3/2}\bigr]$, most negative directly under $q$ (at $s=0$) and dying as $s^{-3}$. Integrating over the plane in polar rings $da=2\pi s\,ds$ (with $s^2=x^2+y^2$),

so the entire image charge reappears as real induced charge. Numerically induced_surface_charge peaks at $\sigma(0,0)=-1.59\times10^{-8}$ C/m² and shrinks outward, while total_induced_charge $=-1.006\times10^{-9}\approx-q$.

In the region $z>0$ the real field is identical to that of $q$ at $+d$ together with the image $-q$ at $-d$, so the force on $q$ is simply its Coulomb attraction to that image, a distance $2d$ away:

directed toward the plane. The energy is not that of a real pair, however: the image field exists only in $z>0$, which holds just half of the $q$–image field energy. Integrating this force in from infinity (equivalently, halving the pair energy $-kq^2/2d$) gives

half of $-kq^2/2d$. Numerically image_force(q, d) $=-2.247\times10^{-7}$ N $=-k\,q^2/(2d)^2$ (with $k=$ K_E) — negative, i.e. attractive toward the conductor.

Separating $V=X(x)Y(y)$ turns $\nabla^2V=0$ into $X''/X=-Y''/Y=k^2$. The grounded plates $V(x,0)=V(x,a)=0$ select $Y=\sin(n\pi y/a)$ with $k=n\pi/a$, and decay as $x\to\infty$ picks $X=e^{-n\pi x/a}$. Imposing the mouth condition $V(0,y)=V_0$ is a Fourier sine problem; sine orthogonality gives $C_n=\tfrac2a\int_0^a V_0\sin\tfrac{n\pi y}{a}\,dy=4V_0/n\pi$ for odd $n$ and $0$ for even, so

Summing with $\sum_{n\,\text{odd}}\tfrac{w^n}{n}=\tfrac12\ln\tfrac{1+w}{1-w}$ at $w=e^{-\pi(x-iy)/a}$ and taking the imaginary part collapses this to

On $y=0,a$ the numerator vanishes so $V=0$; as $x\to\infty$, $\sinh\to\infty$ so $V\to0$; as $x\to0^+$ the ratio $\to\infty$ and $\arctan\to\tfrac\pi2$, giving $V\to V_0$. The two forms agree to $\sim10^{-4}$ (e.g. $V(0.3,0.5)=4.731$ both ways), and slot_potential_closed returns $0$ on $y=0$ and $\sim10^{-16}$ on $y=a$ — the grounded plates.

Discretizing $\nabla^2V=0$ with the 5-point Laplacian gives $V_{i+1,j}+V_{i-1,j}+V_{i,j+1}+V_{i,j-1}-4V_{i,j}=0$, i.e. the mean-value property

Each Gauss–Seidel sweep overwrites every interior node with its neighbour-average while the boundary stays pinned. Any error field is itself discrete-harmonic and zero on the boundary, so by the discrete maximum principle it can have no interior extremum and must vanish — the iteration contracts to the one interior consistent with the boundary, exactly the first uniqueness theorem. For boundary data $V=V_0x/L$, the linear $V=V_0x/L$ is harmonic ($\partial_x^2=\partial_y^2=0$) and already obeys the averaging, so it is that unique interior. Hence solve_laplace_2d(lambda x, y: 10.0*x) relaxes to $V(0.65,0.35)=6.500=10\times0.65$ at every node, independent of the zeroed initial guess.

The two charges sum to $Q=q+(-q)=0$, so the cloud is neutral. The dipole moment (Gr Eq. 3.98) is $\mathbf p=\sum_a q_a\mathbf r'_a=q(\tfrac a2\hat{\mathbf z})+(-q)(-\tfrac a2\hat{\mathbf z})=qa\,\hat{\mathbf z}$, pointing from $-q$ to $+q$. With $Q=0$ the monopole term drops and the $n=1$ term leads; using $P_1(\cos\alpha)=\cos\alpha$ and $\mathbf p\cdot\hat{\mathbf r}=qa\cos\theta$,

For $q=10^{-9},\,a=0.01$ this is $\mathbf p=(0,0,10^{-11})$ C·m. On the axis at $r=100a$ the dipole term gives $V=8.9876\times10^{-2}$ V against the exact $8.9878\times10^{-2}$ V — agreement to $2\times10^{-5}$, matching monopole_moment$=0$, dipole_moment$=(0,0,qa)$, and multipole_potential(dip,1)$\approx$potential_point_charges.

Since $Q=0$, the $n=0$ term is $V_{\text{mono}}=kQ/r=0$ — it captures none of the potential, so its relative error is exactly $100\%$. Each added term costs another factor $r'/r\lesssim a/r$, so truncating after the $\ell$-th term leaves a relative error $\sim(a/r)^{\ell+1}$. For the dipole ($\ell=1$) that is $(a/r)^{2}$; setting $(a/r)^{2}\lesssim0.01$ gives $r\gtrsim10a$. At $r=30a$ the measured monopole error is $1.0$ and the dipole-order error is $2.78\times10^{-4}$ ($0.028\%$) — strictly smaller and far under $1\%$. (The symmetric pair's even moments vanish, so the leading omitted term is the octupole and the error falls as $(a/2r)^{2}=2.78\times10^{-4}$ at $r=30a$, exactly the value multipole_potential prints.)

Under $\mathbf r'\to\mathbf r'-\mathbf a$,

The shift $-Q\mathbf a$ vanishes iff $Q=0$. For the neutral pair $Q=0$, and translating by $\mathbf a=(0.3,0.2,0.5)$ leaves $\mathbf p=(0,0,10^{-11})$ unchanged. For the like pair $(+q,+q)$, $Q=2q=2\times10^{-9}\neq0$, so $\mathbf p$ moves from $\mathbf 0$ to $-Q\mathbf a=(-6\times10^{-10},-4\times10^{-10},-10^{-9})$. This is exactly what the code shows: dipole_moment is invariant for the neutral pair but shifts by $-$monopole_moment$\cdot\mathbf a$ for the like pair.

Write $\mathbf p=p\,\hat{\mathbf z}$, so $\mathbf p\cdot\hat{\mathbf r}=p\cos\theta$. On the axis $\hat{\mathbf r}=\hat{\mathbf z}$ ($\theta=0$):

purely radial and along $\mathbf p$. On the perpendicular bisector $\hat{\mathbf r}\perp\mathbf p$ ($\theta=90^\circ$, $\mathbf p\cdot\hat{\mathbf r}=0$): $\mathbf E=\dfrac{k}{r^{3}}[\,0-p\,\hat{\mathbf z}\,]=-\dfrac{kp}{r^{3}}\hat{\mathbf z}$, antiparallel and half as strong — the $2:1$ ratio. With $p=10^{-11},\,r=0.1$, dipole_field(p) returns $E_z=179.75=2kp/r^{3}$ on the $z$-axis and $E_z=-89.876=-kp/r^{3}$ on the $x$-axis. $\mathbf E$ is purely radial on the axis ($\mathbf E\parallel\hat{\mathbf r}$) and purely transverse on the bisector ($\mathbf E\perp\hat{\mathbf r}$).

Monopole: $Q=q-2q+q=0$. Dipole: $\mathbf p=q(a\hat{\mathbf z})-2q(\mathbf 0)+q(-a\hat{\mathbf z})=\mathbf 0$. Both vanish, so the $n=2$ quadrupole leads and $V\sim1/r^{3}$. The charges sit on the $z$-axis, so $x_a=y_a=0$ and $r_a^{2}=z_a^{2}$. Then

while $Q_{xx}=Q_{yy}=\sum_a q_a(0-z_a^{2})=-2qa^{2}=-\tfrac12 Q_{zz}$ (traceless). With $q=10^{-9},\,a=0.01$: $Q_{zz}=4\times10^{-13}$ and $Q_{xx}=Q_{yy}=-2\times10^{-13}$, matching monopole_moment$=$dipole_moment$=0$ and quadrupole_moment $Q[2][2]=4qa^{2}$, $Q[0][0]=-2qa^{2}$.

Use the pairwise form $W=\tfrac12\sum_i q_iV_i=\sum_{i<j}\dfrac{kq_iq_j}{r_{ij}}$. The square has $\binom{4}{2}=6$ distinct pairs: the four edges at separation $a$ and the two diagonals at separation $a\sqrt2$. With all charges equal to $q$,

Because the sum runs only over $i<j$, no term ever has $i=j$, so the divergent $1/r_{ii}=1/0$ self-energy of each point charge never enters. With $q=a=1$ this is $k(4+\sqrt2)=8.988\times10^9\times5.4142=4.866\times10^{10}$ J, exactly what work_to_assemble returns for the four corner charges.

With $u=\tfrac{\varepsilon_0}{2}E^2$ and the spherical shell $d\tau=4\pi r^2\,dr$, and using $4\pi\varepsilon_0 k=1$ so that $\tfrac{\varepsilon_0}{2}\,4\pi k^2=\tfrac{k}{2}$:

Adding, $W=\bigl(\tfrac12+\tfrac1{10}\bigr)kQ^2/R=\tfrac35\,kQ^2/R$. Numerically (for $Q=10^{-9}$ C, $R=0.05$ m) self_energy_uniform_sphere gives $1.0785\times10^{-7}$ J, which field_energy_spherical reproduces as $1.0758\times10^{-7}$ J — agreeing to $0.25\%$ once the upper limit is pushed to $\sim2000R$ to capture the slow $1/r^4$ tail.

The gap field is uniform, $E=V/d$, so the energy density $u=\tfrac{\varepsilon_0}{2}E^2$ is constant and the stored energy is $u$ times the gap volume $Ad$:

Identifying the parallel-plate capacitance $C=\varepsilon_0 A/d$,

so the field picture (Eq. 2.45) and the circuit picture (Eq. 2.55) land on one number. For $A=0.01$ m², $d=1$ mm, $V=12$ V this is $C=88.5$ pF and $\tfrac12CV^2=6.375\times10^{-9}$ J; field_energy over the $A\times d$ box returns the identical $6.375\times10^{-9}$ J — exact here because the constant integrand makes the midpoint quadrature exact.

Place $\pm Q$ on the conductors, integrate the ~EM-02 field across the gap for $V$, and form $C=Q/V$. Spheres ($E=kQ/r^2$ between $a$ and $b$):

As $b\to\infty$, $ab/(b-a)\to a$, recovering the isolated sphere $C=4\pi\varepsilon_0 a$. Cylinders ($E=\lambda/2\pi\varepsilon_0 s$, $\lambda=Q/L$):

which grows with $L$ and with a tighter gap (smaller $\ln(b/a)$). Numerically capacitance_spherical(0.05, b) runs $1.11\times10^{-11}\to5.62\times10^{-12}\to5.564\times10^{-12}$ F for $b=0.1,5,500$ m, closing on capacitance_isolated_sphere(0.05) $=5.563\times10^{-12}$ F; and capacitance_cylindrical rises from $5.06\times10^{-11}$ F at $L=1$ m to $1.01\times10^{-10}$ F at $L=2$ m (and to $8.03\times10^{-11}$ F as the gap tightens to $b=2$ mm).

The surface layer itself produces the jump from $0$ inside to $\sigma/\varepsilon_0$ outside; the field acting on the layer omits its own contribution and is the average of the two sides, $E_{\text{avg}}=\tfrac12(\sigma/\varepsilon_0+0)=\sigma/2\varepsilon_0$. The force per area is the charge density times this field:

Because it goes as $\sigma^2$ it is positive whatever the sign of $\sigma$, so the pull is always outward, into the field. For one plate of the P3 capacitor ($\sigma=\varepsilon_0 V/d$) this is the attraction toward the opposite plate. Numerically surface_pressure(1e-6) $=5.647\times10^{-2}$ Pa, equal to $\tfrac{\varepsilon_0}{2}(\sigma/\varepsilon_0)^2$ — the same result in field form.

With $\mathbf P=P\,\hat{\mathbf z}$ constant, the bound volume density vanishes everywhere,

so all bound charge sits on the surfaces through $\sigma_b=\mathbf P\cdot\hat{\mathbf n}$. On the top face $\hat{\mathbf n}=+\hat{\mathbf z}$ gives $\sigma_b=P(\hat{\mathbf z}\cdot\hat{\mathbf z})=+P$; on the bottom $\hat{\mathbf n}=-\hat{\mathbf z}$ gives $-P$; on any side $\hat{\mathbf n}\perp\hat{\mathbf z}$ gives $0$. The two sheets carry $\pm P A$ over equal-area faces, so the slab is neutral. With $P=2\times10^{-6}\ \mathrm{C/m^2}$ the code returns $\sigma_b=+2\times10^{-6}$ (top), $-2\times10^{-6}$ (bottom), $0$ (side), and bound_volume_charge $=0$ — the pair-of-sheets picture exactly.

Write $\mathbf P=k\,\mathbf r=k(x,y,z)$. Its divergence is

At the surface $\hat{\mathbf n}=\hat{\mathbf r}$, and $\mathbf P\cdot\hat{\mathbf r}=k\,\mathbf r\cdot\hat{\mathbf r}=kr=kR$, so $\sigma_b=kR$. The totals cancel:

as required for a dielectric with no free charge. The code (e.g. $k=3,\ R=2$) returns bound_volume_charge $=-9=-3k$ and bound_surface_charge $=6=kR$, with the weighted sum zero to rounding.

Spherical symmetry forces $\mathbf D=D(r)\,\hat{\mathbf r}$, so Gauss's law for $\mathbf D$ on a concentric sphere gives

No permittivity enters — $\mathbf D$ is sourced by free charge alone and is the same in every medium. Dividing by $\varepsilon=\varepsilon_r\varepsilon_0$ recovers the field,

the vacuum Coulomb field screened by $\varepsilon_r$. Because the enclosed free charge is the same at every radius, free_charge_enclosed returns $q_f$ ($\approx5.0003\times10^{-9}$ C for $q_f=5$ nC) at $R=0.1,0.2,0.5,1.0$ alike — the tiny excess is the 80-point flux quadrature, not physics.

Start from the definition $\mathbf D=\varepsilon_0\mathbf E+\mathbf P$ and substitute the linear law $\mathbf P=\varepsilon_0\chi_e\mathbf E$:

with $\varepsilon=\varepsilon_0(1+\chi_e)=\varepsilon_0\varepsilon_r$ and $\chi_e=\varepsilon_r-1$. The two expressions are not two fields but one vector regrouped: $\varepsilon\mathbf E$ bundles the vacuum term $\varepsilon_0\mathbf E$ together with the material response $\mathbf P=\varepsilon_0\chi_e\mathbf E$ into a single constant times $\mathbf E$. For $\varepsilon_r=4$ ($\chi_e=3$) both displacement_linear(eps_r, E) and displacement_field(E, polarization_linear(eps_r, E)) return the identical $D_x=1.5915\times10^{-7}$ at the test point — the equality the tests check.

With $\mathbf P=P\,\hat{\mathbf z}$ and $\hat{\mathbf n}=\hat{\mathbf r}$,

The net bound charge vanishes, $\oint\sigma_b\,da=PR^2\!\int_0^{2\pi}\!d\phi\int_0^\pi\cos\theta\sin\theta\,d\theta=0$, since $\int_0^\pi\cos\theta\sin\theta\,d\theta=0$. For the interior, model the sphere as overlapping balls of charge $\pm\rho$ displaced by $\mathbf s$ with $\mathbf P=\rho\mathbf s$. Inside a uniform ball Gauss gives $\mathbf E=\rho\mathbf r/3\varepsilon_0$; superposing the positive ball (centre $\mathbf s$) and the negative ball (centre $0$),

a uniform depolarizing field. The code confirms polarized_sphere_surface_charge(P)(θ) $=P\cos\theta$ (integrating to $0$) and polarized_sphere_inner_field((0,0,P)) $=(0,0,-P/3\varepsilon_0)$, e.g. $-3.7647\times10^{4}$ V/m for $P=10^{-6}$.

The cross product $\mathbf v\times\mathbf B$ is perpendicular to both factors, so $\mathbf F\cdot\mathbf v = Q(\mathbf v\times\mathbf B)\cdot\mathbf v = 0$ (a vector dotted with a cross product that contains it). Hence the work done is

so $|\mathbf v|$ is constant. For $\mathbf v = v\,\hat{\mathbf x}$ and $\mathbf B = B\,\hat{\mathbf z}$, $\hat{\mathbf x}\times\hat{\mathbf z} = -\hat{\mathbf y}$, giving $\mathbf F = -QvB\,\hat{\mathbf y}$. Numerically $QvB = (1.602\times10^{-19})(10^{6})(0.5) = 8.01\times10^{-14}\,$N, so $\mathbf F = (0,\,-8.01\times10^{-14},\,0)\,$N — exactly what lorentz_force returns, with dot(F, v) and dot(F, B) both zero. With the speed fixed, the magnetic force supplies the centripetal force, $Qv_\perp B = mv_\perp^2/r$, so $r = mv_\perp/QB$ and $\omega = v_\perp/r = QB/m$.

Each element feels $d\mathbf F=I\,d\boldsymbol\ell\times\mathbf B$, which for a straight segment $\mathbf L$ in uniform $\mathbf B$ integrates to $\mathbf F=I(\mathbf L\times\mathbf B)$. With $\mathbf L=L\hat{\mathbf y}$, $\mathbf B=B\hat{\mathbf z}$ and $\hat{\mathbf y}\times\hat{\mathbf z}=\hat{\mathbf x}$, this is $\mathbf F=ILB\,\hat{\mathbf x}$. Let $\mathbf B$ now be the field of a parallel wire $I_2$ a distance $d$ off, $B=\mu_0 I_2/2\pi d$; a length $\ell$ of wire $I_1$ then feels

and the right-hand rule points the force on $I_1$ toward $I_2$ — parallel currents attract. Numerically force_on_wire(10, (0,1,0), (0,0,0.5)) $=(5,0,0)$ N $=ILB$; feeding infinite_wire_field(10) at $d=0.05$ m (so $B=4\times10^{-5}$ T) into force_on_wire gives $F/\ell=4\times10^{-4}$ N/m pointing at the other wire — exactly $\mu_0 I_1 I_2/2\pi d$.

Put the wire on the $z$-axis and the field point at distance $s$ in the mid-plane. With $d\boldsymbol\ell'=dz'\,\hat{\mathbf z}$ and separation length $\eta=\sqrt{s^2+{z'}^2}$, the cross product $d\boldsymbol\ell'\times\hat{\boldsymbol\eta}$ is azimuthal with magnitude $\sin\theta\,dz'$ where $\sin\theta=s/\eta$:

since $\int_{-\infty}^{\infty}s\,dz'/(s^2+{z'}^2)^{3/2}=2/s$, with the field circling the wire azimuthally. For $I=10$ A, infinite_wire_field gives exactly $\mu_0 I/2\pi s$ ($1.000\times10^{-4}$ T at $s=0.02$ m); approximating the wire by biot_savart over $[-L,L]$ climbs toward it — $9.95\times10^{-5}$, $9.998\times10^{-5}$, $9.99992\times10^{-5}$ T for $L=0.2,1,5$ m — converging once $L\gg s$.

On the axis every element $I\,d\boldsymbol\ell'$ sits the same distance $\eta=\sqrt{R^2+z^2}$ from the field point with $d\boldsymbol\ell'\perp\boldsymbol\eta$, so $|d\mathbf B|=\tfrac{\mu_0 I}{4\pi}\,d\ell'/(R^2+z^2)$. Symmetry cancels the transverse parts around the loop, leaving only the $z$-projection, weighted by $\cos\alpha=R/\sqrt{R^2+z^2}$:

At the centre $z=0$ this is $\mu_0 I/2R$; for $z\gg R$ it falls as $\mu_0 I R^2/2z^3$ — a magnetic dipole. Numerically circular_loop_field (Biot–Savart) and loop_axis_field_closed agree on the axis — $3.142\times10^{-5}$ T at the centre ($=\mu_0 I/2R$ for $I=5$ A, $R=0.1$ m), $1.111\times10^{-5}$ at $z=0.1$, $9.935\times10^{-7}$ at $z=0.3$ — with the transverse components vanishing.

Symmetry makes $\mathbf B$ azimuthal with constant magnitude on a circle of radius $s$ about the wire, so $\oint\mathbf B\cdot d\boldsymbol\ell=B\oint d\ell=B(2\pi s)$. Ampère's law sets this equal to $\mu_0 I_{\text{enc}}=\mu_0 I$:

the Biot–Savart result in one line, the radius dropping out of the circulation. Its companion $\nabla\cdot\mathbf B=0$ says $\mathbf B$ has no sources or sinks — field lines never start or stop, so there are no magnetic monopoles (the contrast with $\nabla\cdot\mathbf E=\rho/\varepsilon_0$). Numerically ampere_circulation(infinite_wire_field(10), s) returns $\mu_0 I=1.257\times10^{-5}$ for every radius (e.g. $s=0.03$ m), and div_B_residual $\approx9.6\times10^{-8}\approx0$ — no monopoles.

For $\mathbf A=A_z(s)\hat{\mathbf z}$ the only surviving curl component is azimuthal, $\nabla\times\mathbf A=-\partial_s A_z\,\hat{\boldsymbol\varphi}$. With $A_z=-\frac{\mu_0 I}{2\pi}\ln(s/s_0)$,

the EM-08 wire field. Replacing $s_0\to s_0'$ shifts $A_z$ by the constant $\frac{\mu_0 I}{2\pi}\ln(s_0'/s_0)$, i.e. $\mathbf A\to\mathbf A+\nabla\lambda$ with $\lambda=\big[\tfrac{\mu_0 I}{2\pi}\ln(s_0'/s_0)\big]z$ linear in $z$; since $\nabla\times\nabla\lambda=0$, $\mathbf B$ is untouched. At $(0.05,0,0)$ with $I=10$ A, B_from_A(wire_vector_potential(I)) returns $(0,\,4.0\times10^{-5},\,0)$ T, matching infinite_wire_field for every s0.

Since $\mathbf m=I\int d\mathbf a=I\mathbf a$ is current times the oriented area, only the enclosed area enters — the boundary shape is irrelevant. (a) Circle: $\mathbf a=\pi R^2\hat{\mathbf z}$, so $\mathbf m=I\pi R^2\hat{\mathbf z}$. (b) Square: $\mathbf a=L^2\hat{\mathbf z}$, so $\mathbf m=IL^2\hat{\mathbf z}$. Imposing $\pi R^2=L^2$ makes the two areas equal, hence $\mathbf m$ identical. For $I=3$ A, $R=0.02$ m (so $L=\sqrt\pi\,R$), both magnetic_dipole_moment calls return $m_z=3.770\times10^{-3}$ A·m² — equal to machine precision, confirming the moment depends on area alone.

For $\mathbf m=m\hat{\mathbf z}$ use $\hat{\mathbf z}\times\hat{\mathbf r}=\sin\theta\,\hat{\boldsymbol\varphi}$:

purely azimuthal. Along a fixed direction $\theta=\text{const}$ the magnitude $|\mathbf A|=\frac{\mu_0}{4\pi}\frac{m\sin\theta}{r^2}\propto r^{-2}$, so doubling $r$ divides $|\mathbf A|$ by $2^2=4$. On the equatorial line $(r,0,0)$, $\theta=\pi/2$ and $\sin\theta=1$. With $m=1$, dipole_vector_potential gives $|\mathbf A|=1.0\times10^{-5}$ at $(0.1,0,0)$ and $2.5\times10^{-6}$ at $(0.2,0,0)$, ratio exactly $4.0$.

On axis ($\theta=0$) $\hat{\mathbf r}=\hat{\mathbf z}$ and $\mathbf m\cdot\hat{\mathbf r}=m$, so

On the equator ($\theta=\pi/2$) $\mathbf m\cdot\hat{\mathbf r}=0$, leaving $\mathbf B=-\frac{\mu_0 m}{4\pi r^3}\hat{\mathbf z}$ — half as large and antiparallel, the same $2{:}1$ ratio the electric dipole shows under $(1/4\pi\varepsilon_0,\mathbf p)\to(\mu_0/4\pi,\mathbf m)$. With $m=1$, $d=0.1$, dipole_B_field_closed gives $\mathbf B(0,0,d)=(0,0,2.0\times10^{-4})$ and $\mathbf B(d,0,0)=(0,0,-1.0\times10^{-4})$, reproduced pointwise by B_from_A(dipole_vector_potential(m)).

Write $\mathbf A_{\text{dip}}=\frac{\mu_0}{4\pi}\frac{\mathbf m\times\mathbf r}{r^3}$ and apply $\nabla\cdot(f\mathbf F)=f\,\nabla\cdot\mathbf F+\mathbf F\cdot\nabla f$ with $f=r^{-3}$, $\mathbf F=\mathbf m\times\mathbf r$. The first term vanishes because $\nabla\cdot(\mathbf m\times\mathbf r)=\mathbf r\cdot(\nabla\times\mathbf m)-\mathbf m\cdot(\nabla\times\mathbf r)=0$ (both curls are zero); the second vanishes because $(\mathbf m\times\mathbf r)\cdot\nabla r^{-3}=-3r^{-5}(\mathbf m\times\mathbf r)\cdot\mathbf r=0$ (the cross product is $\perp\mathbf r$). Hence $\nabla\cdot\mathbf A_{\text{dip}}=0$ identically — the Coulomb gauge already holds. A shift $\mathbf A\to\mathbf A+\nabla\lambda$ leaves $\mathbf B=\nabla\times\mathbf A$ unchanged for any $\lambda$, but preserves $\nabla\cdot\mathbf A=0$ only when $\nabla^2\lambda=0$ — harmonic $\lambda$, e.g. $\lambda\propto z$. Accordingly coulomb_gauge_residual returns $\sim10^{-8}$ at the off-axis points $(0.1,0.05,0.08)$ etc., numerically zero.

Use the bound currents (Eqs. 6.13–6.14) $\mathbf J_b=\nabla\times\mathbf M$ and $\mathbf K_b=\mathbf M\times\hat{\mathbf n}$. For $\mathbf M=M\,\hat{\mathbf z}$ uniform, every derivative of $\mathbf M$ is zero, so $\mathbf J_b=\nabla\times(M\,\hat{\mathbf z})=\mathbf 0$ — no interior current. On the side wall the outward normal is radial, $\hat{\mathbf n}=\hat{\mathbf s}$, and with $\hat{\mathbf z}\times\hat{\mathbf s}=\hat{\boldsymbol\phi}$,

This is precisely a solenoid's wall current with $nI=M$, so the long-solenoid result gives $\mathbf B=\mu_0 nI\,\hat{\mathbf z}=\mu_0 M\,\hat{\mathbf z}$ inside. Numerically the curl returns $(0,0,0)$ and bound_surface_current((0,0,M),(1,0,0)) $=(0,M,0)$ with norm $=M$.

For $\mathbf M=(-a\,y,\;a\,x,\;0)$ the bound volume current (Eq. 6.13) is the curl

The first two components vanish ($M_z=0$ and $M_x,M_y$ have no $z$-dependence); the third is $\partial_x(a x)-\partial_y(-a y)=a-(-a)=2a$, so

A non-uniform $\mathbf M$ thus carries a uniform axial volume current, unlike the uniform $\mathbf M$ of P1. For $a=5$, bound_volume_current(...) returns $(0,0,10)=2a\,\hat{\mathbf z}$.

Apply Ampère's law for $\mathbf H$ (Eq. 6.20), $\oint\mathbf H\cdot d\boldsymbol\ell=I_{f,\text{enc}}$, to a rectangular loop with one side of length $\ell$ running along the axis inside the solenoid. Only free current is enclosed — the $n\ell I_f$ turns it threads — and $\mathbf H$ is axial and uniform inside, zero outside, so $H\ell=n\ell I_f$ gives

Then $\mathbf B=\mu\mathbf H=\mu_0(1+\chi_m)nI_f\,\hat{\mathbf z}$ (Eq. 6.31) and $\mathbf M=\chi_m\mathbf H=\chi_m nI_f\,\hat{\mathbf z}$ (Eq. 6.29). Substituting back into the definition,

recovering the original field — exactly what feeding B_linear and magnetization_linear through auxiliary_field_H returns.

Inside the uniformly magnetized sphere the fields are uniform (Eq. 6.16), $\mathbf B_{\text{in}}=\tfrac{2}{3}\mu_0\mathbf M$. The auxiliary field follows from its definition (Eq. 6.18):

Check the constitutive relation: $\mu_0(\mathbf H_{\text{in}}+\mathbf M)=\mu_0\!\left(-\tfrac13\mathbf M+\mathbf M\right)=\tfrac23\mu_0\mathbf M=\mathbf B_{\text{in}}$. ✓ Because $\mathbf H_{\text{in}}=-\mathbf M/3$ points against $\mathbf M$ it is a demagnetizing field — sourced by the effective magnetic charge $-\nabla\cdot\mathbf M$ on the poles. For $\mathbf M=8\times10^5\,\hat{\mathbf z}$ A/m the code returns $\mathbf B_{\text{in}}=0.6702\,\hat{\mathbf z}$ T and $\mathbf H_{\text{in}}=-2.667\times10^5\,\hat{\mathbf z}$ A/m, matching magnetized_sphere_inner_B and magnetized_sphere_inner_H.

(a) The relative permeability is $\mu/\mu_0=1+\chi_m$ (Eq. 6.33). Reading the sign and size of $\chi_m$: water $\chi_m=-1.7\times10^{-5}<0$ → diamagnetic, $\mu/\mu_0=0.999983$; the metal $\chi_m=2.5\times10^{-4}$ (small, positive) → paramagnetic, $\mu/\mu_0=1.00025$; soft iron $\chi_m=5500$ (large) → ferromagnetic, $\mu/\mu_0=5501$. (b) The descending branch is $M_\downarrow(H)=M_s\tanh\!\big((H+H_c)/w\big)$. At $H=0$,

the remanent (permanent) magnetization, while $M_\downarrow(-H_c)=M_s\tanh 0=0$ fixes the coercivity $H_c=5\times10^3$ A/m. Since $M_\downarrow(0)>0>M_\uparrow(0)$ the loop is open — $\mathbf M$ at $H=0$ depends on whether the field was swept down or up, i.e. on history. These match remanence ($7.893\times10^5$), coercivity ($5\times10^3$), classify_material and permeability.

With $\mathbf B$ uniform and normal to the loop, $d\mathbf a=da\,\hat{\mathbf z}$, so

Faraday's law then gives $\mathcal{E}=-\dfrac{d\Phi}{dt}=-B_0A\omega\cos\omega t$, peaking at $|\mathcal{E}|=B_0A\omega$ when $\cos\omega t=\pm1$. While the flux rises ($\dot\Phi=B_0A\omega\cos\omega t>0$) Lenz's law forces the induced current to oppose it, so $\mathcal{E}<0$. Numerically, at $\sin\omega t=1$ the flux equals $B_0A=0.5\times0.01=5\times10^{-3}$ (matching magnetic_flux), and at $t=0$ faraday_emf returns $-0.5\,$V $=-B_0A\omega$; lenz_sign returns $-1$ for $\dot\Phi>0$ and $+1$ for $\dot\Phi<0$, flipping with $\cos\omega t$.

Each charge in the rod feels $\mathbf f=q\,\mathbf v\times\mathbf B$; with $\mathbf v$, $\mathbf B$ and the rod mutually perpendicular, $|\mathbf v\times\mathbf B|=vB$ points along the rod. The EMF is the work per unit charge around the loop,

Equivalently, in time $dt$ the rod sweeps area $dA=Lv\,dt$, so $d\Phi=B\,dA=BLv\,dt$ and $|{-}d\Phi/dt|=BvL$ — the flux rule gives the identical answer. By Lenz the induced current circulates so its own flux opposes the growing $\Phi$ (counter-clockwise if $\mathbf B$ points out of the page and the area is increasing). Numerically $\mathcal{E}=0.3\times2\times0.5=0.3\,$V, as motional_emf(0.3, 2.0, 0.5) returns.

Inside the solenoid the field is uniform, $B=\mu_0(N/l)I$, so the flux through one turn is $BA=\mu_0(N/l)IA$. The total flux linkage threads all $N$ turns,

Comparing with $\Phi=LI$ identifies $L=\mu_0 N^2 A/l$. Inserting the numbers,

exactly what solenoid_inductance(1000, 1e-4, 0.2) returns ($6.2832\times10^{-4}$ H).

The inductor picture gives $W=\tfrac12 LI^2=\tfrac12(6.28\times10^{-4})(3)^2=2.83\times10^{-3}\,$J. The field picture integrates $u=B^2/2\mu_0$ over the interior volume $V=Al$, where $B=\mu_0NI/l$ is uniform:

so the two are algebraically identical — the substitution turns the field integral straight back into $\tfrac12 LI^2$. Numerically both equal $2.8274\times10^{-3}\,$J, confirming energy_in_inductor(L, 3.0) $=$ magnetic_field_energy_solenoid(1000, 1e-4, 0.2, 3.0), the magnetic twin of the ~EM-06 charge-vs-field duality.

Current $I_1$ in coil 1 makes the uniform interior field $B_1=\mu_0(N_1/l)I_1$, which threads every one of coil 2's $N_2$ turns. The flux linkage of coil 2 is therefore

and $\mathcal{E}_2=-M\,dI_1/dt$ follows from Faraday's law. The expression is symmetric under $N_1\leftrightarrow N_2$, so driving coil 2 instead gives the identical $M$: this is reciprocity $M_{12}=M_{21}$ (a general theorem, here made explicit). Hence mutual_inductance_solenoids(N1, N2, A, l) equals the same call with $N_1,N_2$ swapped (e.g. both give $3.770\times10^{-3}$ H for $N_1{=}2000,N_2{=}3000$).

In series the same current flows through every element, so the impedances add: $Z=Z_R+Z_L+Z_C=R+i\omega L+\frac{1}{i\omega C}$. Using $\frac{1}{i\omega C}=-\frac{i}{\omega C}$,

The reactance $X$ vanishes when $\omega L=1/\omega C$, i.e. $\omega^2=1/LC$, giving $\omega_0=1/\sqrt{LC}$. There $Z=R$ is purely real, and since $\lvert Z\rvert=\sqrt{R^2+X^2}\ge R$ with equality only at $X=0$, the magnitude is minimal, $\lvert Z\rvert=R$. So series_rlc_impedance(R, L, C, resonant_frequency(L, C)) returns $R+0i$ (here $10+0i\,\Omega$), while off resonance the imaginary part is nonzero.

Phasor Ohm's law gives $\lvert I\rvert=V_0/\lvert Z\rvert$ with $\lvert Z\rvert=\sqrt{R^2+X^2}$, $X=\omega L-1/\omega C$. The denominator is smallest when the reactance vanishes ($X=0$ at $\omega_0$), so the current is largest there:

Any nonzero $X$ only adds to the denominator, so $\lvert I\rvert$ falls off on both sides — the Lorentzian resonance peak shared with the ~CM-15 oscillator, of height $V_0/R$ and width set by $R$. With $R=10,\,V_0=5$ the peak is $\lvert I\rvert=V_0/R=0.5$ A; at $0.7\,\omega_0$ and $1.3\,\omega_0$ the current has dropped to $0.199$ A and $0.256$ A — both below the peak, so current_amplitude is maximal exactly at resonant_frequency(L, C).

The phasor current $I_0=V_0/Z$ carries the phase $\phi=-\arg Z=-\arctan(X/R)$ with reactance $X=\omega L-1/\omega C$:

Below $\omega_0$ the capacitor dominates ($1/\omega C>\omega L$, so $X<0$): $\phi>0$, the current leads the drive — capacitive. Above $\omega_0$ the inductor dominates ($X>0$): $\phi<0$, the current lags — inductive. At $\omega_0$ the reactances cancel ($X=0$) and $\phi=0$, current in phase. So current_phase returns $+78.1^\circ$ at $0.5\,\omega_0$ (leading), $\approx0$ at $\omega_0$, and $-78.1^\circ$ at $2\,\omega_0$ (lagging) — equal magnitudes because $X(r\omega_0) =-X(\omega_0/r)$ is antisymmetric about $\omega_0$ on a log-frequency axis.

At resonance the inductive reactance is $\omega_0 L=L/\sqrt{LC}=\sqrt{L/C}$, so the quality factor — the reactance-to-resistance ratio at $\omega_0$ — has two equal forms:

The half-power width is $\Delta\omega=R/L$; dividing by $\omega_0=1/\sqrt{LC}$,

so $\Delta\omega=\omega_0/Q$ — small $R$ gives large $Q$ and a narrow, sharply selective peak. Under the ~CM-15 dictionary $L\to m,\;R\to b$ this is $Q=\omega_0 m/b$, the mechanical quality factor. With $R=10,\,L=10^{-3},\,C=10^{-6}$: $Q=3.162$ and $\Delta\omega=R/L=10^4$ rad/s, equal to $\omega_0/Q$ — exactly what quality_factor(R, L, C) and bandwidth(R, L) return.

Inductor and capacitor are lossless (their current runs $90^\circ$ out of phase with their voltage), so only $R$ dissipates: with $\lvert I\rvert=V_0/\lvert Z\rvert$,

since $\lvert Z\rvert^2=R^2+X^2$ is minimal ($=R^2$) at resonance. The power halves when $\lvert Z\rvert^2=2R^2$, i.e. $X^2=R^2$, i.e. $\lvert X\rvert=\lvert\omega L-1/\omega C\rvert=R$. The two roots of $\omega L-1/\omega C=\pm R$ are separated by exactly $\Delta\omega=R/L$, the bandwidth that defines $Q$, and lie at $\omega_0\pm\Delta\omega/2$ to leading order in $1/Q$. With $R=10,\,V_0=5$ the peak is $V_0^2/2R=1.25$ W, and at $\omega_0+\Delta\omega/2$ average_power returns $\approx0.67$ W against the half-peak $0.625$ W — the small excess being the finite-$Q$ asymmetry of the exact half-power roots about $\omega_0$.

Between the plates the conduction current is zero, so Ampère–Maxwell keeps only the displacement term: $\nabla\times\mathbf B=\mu_0\varepsilon_0\,\partial\mathbf E/\partial t=\mu_0\mathbf J_d$ with $\mathbf J_d=\varepsilon_0\,\partial\mathbf E/\partial t$ — the changing field is the sole source of the circulating B. Its magnitude follows directly:

So displacement_current_density((1e12,0,0)) returns $(8.854,0,0)$ — the $x$-component $\varepsilon_0\cdot10^{12}=8.854\ \mathrm{A/m^2}$, the displacement current that closes Ampère's loop.

Take the divergence of the uncorrected $\nabla\times\mathbf B=\mu_0\mathbf J$. The left side vanishes identically, $\nabla\cdot(\nabla\times\mathbf B)=0$, so the law forces $\nabla\cdot\mathbf J=0$. But continuity reads $\nabla\cdot\mathbf J=-\partial\rho/\partial t$, nonzero wherever charge piles up (a charging capacitor plate) — a flat contradiction. Restore consistency by adding $\mathbf J_d=\varepsilon_0\,\partial\mathbf E/\partial t$ and using Gauss $\nabla\cdot\mathbf E=\rho/\varepsilon_0$:

Now $\nabla\cdot(\nabla\times\mathbf B)=\mu_0\nabla\cdot(\mathbf J+\mathbf J_d)=0$ holds for any time-varying fields. The two laws whose consistency this proves are precisely Gauss (gauss_E_residual) and Ampère–Maxwell (ampere_maxwell_residual); together they enforce exactly $\nabla\cdot\mathbf J+\partial\rho/\partial t=0$.

$\mathbf E=E_0\cos(kz-\omega t)\,\hat{\mathbf x}$ depends only on $z,t$, so $\nabla\cdot\mathbf E=\partial E_x/\partial x=0$; likewise $\nabla\cdot\mathbf B=\partial B_y/\partial y=0$ — the two divergence laws hold with no condition. For the curls only the $z$-derivative survives:

so Faraday balances only if $k=\omega/c$. Ampère likewise gives $\nabla\times\mathbf B=(E_0k/c)\sin(kz-\omega t)\,\hat{\mathbf x}$ versus $\mu_0\varepsilon_0\,\partial\mathbf E/\partial t=(E_0\omega/c^2)\sin(kz-\omega t)\,\hat{\mathbf x}$, matching only if $\omega=ck$. All four therefore hold exactly at $\omega=ck$, which is what verify_vacuum_plane_wave(1.0, 1.0, 1.0, …) shows: gauss_E and gauss_B are $0$ identically, and faraday, ampere_maxwell come back $\sim10^{-16}$.

Keep $\mathbf E$ built with $c=1$ (so $\omega=k=1$) but give $\mathbf B$ the speed $c'=2$: amplitude $E_0/c'=\tfrac12$ and frequency $\omega'=c'k=2$. Faraday now pits $\nabla\times\mathbf E=-E_0k\sin(kz-\omega t)\,\hat{\mathbf y}=-\sin(z-t)\,\hat{\mathbf y}$ against

The two carry different space–time phases, $\sin(z-t)$ versus $\sin(z-2t)$, so they cannot cancel at every $(z,t)$ — Faraday is broken. Physically the curl of $\mathbf E$ is what drives $\partial\mathbf B/\partial t$ pointwise, so a self-consistent wave is locked to one frequency and one amplitude ratio $B_0=E_0/c$; mistuning the speed snaps that lock. Numerically faraday_residual(E, Bwrong, …)$=0.79$ and ampere_maxwell_residual$=0.58$, both far past the $0.1$ violation threshold.

In vacuum ($\rho=0,\,\mathbf J=0$) take the curl of Faraday and substitute Ampère–Maxwell:

The identity $\nabla\times(\nabla\times\mathbf E)=\nabla(\nabla\cdot\mathbf E)-\nabla^2\mathbf E$ loses its first term because $\nabla\cdot\mathbf E=0$, collapsing to $-\nabla^2\mathbf E$; the companion $\nabla\cdot\mathbf B=0$ does the same for $\mathbf B$, so both fields obey a genuine wave equation

With $\varepsilon_0=8.854\times10^{-12}$ (~EM-01) and $\mu_0=1.257\times10^{-6}$ (~EM-08), $c=2.998\times10^8$ m/s — the speed of light from two static constants. This is C_SI to machine precision, and gauss_B_residual$\approx0$ for the plane-wave $\mathbf B$ confirms the $\nabla\cdot\mathbf B=0$ the collapse leans on.

With $\mathbf E=E_0\cos(kz)\,\hat x$ and $\mathbf B=(E_0/c)\cos(kz)\,\hat y$, the cross product uses $\hat x\times\hat y=\hat z$:

For the energy density, $1/c^2=\mu_0\varepsilon_0$ makes the magnetic store equal the electric, $\frac1{2\mu_0}(E_0/c)^2=\tfrac{\varepsilon_0}{2}E_0^2$, so

Dividing, $S/u=(1/\mu_0 c)/\varepsilon_0=1/(\mu_0\varepsilon_0 c)=c$, i.e. $S=c\,u$. For $E_0=1000$ V/m at a crest the code returns $|\mathbf S|=2654$ W/m² and $C\cdot u=2654$ W/m² — equal, with $\mathbf S\parallel+\hat z$.

The plane-wave relation $B=E/c$ with $c^2=1/\mu_0\varepsilon_0$ turns the magnetic density into the electric one:

The two halves of $u$ are therefore equal, and

Numerically (with $E_0=1000$ V/m) the code reports electric $u_E=4.427\times10^{-6}$ J/m³ equal to magnetic $u_B=4.427\times10^{-6}$ J/m³, summing to $u=8.854\times10^{-6}$ J/m³ = energy_density(E,B).

The momentum density is $g=S/c^2=1361/(2.998\times10^8)^2=1.51\times10^{-14}\ \mathrm{kg/(m^2\,s)}$. A black panel absorbs the beam, gaining its momentum at rate $S/c$ per unit area, so

A mirror sends the beam back, reversing the momentum, so it must supply twice the impulse, $P_{\text{refl}}=2S/c=9.08\ \mu\text{Pa}$ — the factor 2 is exactly $\Delta p=2p$ for reflection versus $\Delta p=p$ for absorption. These reproduce radiation_pressure(1361) $=4.54\times10^{-6}$ Pa and radiation_pressure(1361, reflected=True) $=9.08\times10^{-6}$ Pa.

Only $E_x$ and $B_y$ are nonzero, so $E^2=E_x^2$, $B^2=B_y^2$, and every product $E_iE_j$ (and $B_iB_j$) vanishes unless $i=j=x$ (resp. $i=j=y$) — so all off-diagonal $T_{ij}=0$. Using equipartition $\tfrac12\varepsilon_0E_x^2=\tfrac1{2\mu_0}B_y^2$ (P2),

and likewise $T_{yy}=-\tfrac12\varepsilon_0E_x^2+\tfrac1{2\mu_0}B_y^2=0$, while

So $\mathbf T=\mathrm{diag}(0,0,-u)$; the code gives the symmetric matrix with T[2][2] $=-8.854\times10^{-6}=-$energy_density(E,B)(0,0,0).

At $(a,0,0)$ the azimuthal direction is $\hat\phi=\hat y$, so $\mathbf E=E_z\hat z$ and $\mathbf B=B\,\hat y$. With $\hat z\times\hat y=-\hat x$,

radially inward, magnitude $|\mathbf S|=E_zB/\mu_0$. Integrating the inward flux over the side area $2\pi a L$ and inserting $B=\mu_0I/2\pi a$,

since $V=E_zL$ is the potential drop ($V=IR$). The Joule heat enters through the sides, carried by the field. The code returns $\mathbf S(a,0,0)$ along $-\hat x$ (e.g. $(-79.6,0,0)$ for the sample $E_z,B$) and exactly $\mathbf 0$ when $\mathbf B=\mathbf 0$.

For $f=A\cos(kz-\omega t)$ the two derivatives pull down a factor each:

Substituting into $\partial_t^2 f=v^2\nabla^2 f$ gives $-\omega^2 f=-v^2k^2 f$, i.e. $\omega^2=v^2k^2\Rightarrow\boxed{\omega=vk}$ — true iff the speed matches. For a general $g(z-vt)$, set $\xi=z-vt$: then $\partial_t^2 g=v^2g''(\xi)$ and $\partial_z^2 g=g''(\xi)$, so $\partial_t^2 g=v^2\partial_z^2 g$ identically (d'Alembert). Physically $\omega=vk$ means the phase velocity $\omega/k=v$ is the same for every frequency: the medium is non-dispersive and a pulse keeps its shape. This is why wave_equation_residual $\approx2.2\times10^{-5}$ (numerically $0$) when omega = vk, but jumps to $\approx0.75=O(1)$ when fed 2v.

For $\mathbf E=\mathbf E_0e^{i(\mathbf k\cdot\mathbf r-\omega t)}$ the operators become $\nabla\to i\mathbf k$, $\partial_t\to-i\omega$. Gauss in vacuum, $\nabla\cdot\mathbf E=0$, gives $i\mathbf k\cdot\mathbf E_0=0$, so $\mathbf E_0\perp\hat{\mathbf k}$ (transverse). Faraday $\nabla\times\mathbf E=-\partial\mathbf B/\partial t$ becomes $i\mathbf k\times\mathbf E_0=i\omega\mathbf B_0$, hence

With $\hat{\mathbf k}=\hat{\mathbf z}$ and $\mathbf E_0=E_0\hat{\mathbf x}$, $\hat{\mathbf z}\times\hat{\mathbf x}=\hat{\mathbf y}$, so $\mathbf B_0$ points along $+\hat{\mathbf y}$: $(\hat{\mathbf E},\hat{\mathbf B},\hat{\mathbf k})=(\hat{\mathbf x},\hat{\mathbf y},\hat{\mathbf z})$ is right-handed. For $E_0=1000$ V/m, transverse_B returns $(0,\,3.336\times10^{-6},\,0)$ T, matching $E_0/c=3.336\times10^{-6}$, and is_transverse is True for both fields.

Taking real parts at $z=0$, the tip of $\mathbf E$ moves as

If $\delta=0$ or $\pi$ then $E_y=\pm(A_y/A_x)E_x$ — a fixed straight line: linear (also linear if one amplitude vanishes). For $\delta=\pi/2$, $E_y=A_y\sin\omega t$, so the tip traces $(A_x\cos\omega t,\,A_y\sin\omega t)$, the ellipse $(E_x/A_x)^2+(E_y/A_y)^2=1$; when $A_x=A_y$ this is a circle of radius $A$. Everything else is a tilted ellipse. So $A_x=A_y,\delta=\pi/2$ gives a circle while $A_x=2A_y,\delta=\pi/2$ gives a 2:1 ellipse — exactly classify_polarization returning "circular", "linear", "elliptical" for the three checks.

The index is $n=\sqrt{\varepsilon_r\mu_r}=\sqrt{2.25}=1.5$, so the phase velocity is $v=c/n=2.998\times10^8/1.5=1.999\times10^8$ m/s. Since $\omega$ (hence the frequency) is conserved across the interface, $\lambda=v/f=(c/n)/(c/\lambda_0)=\lambda_0/n =500/1.5=333.3$ nm. For the intensity, the Poynting vector is $\mathbf S=\frac1{\mu_0}\mathbf E\times\mathbf B$ with $B=E/c$, so $S=\frac{E^2}{\mu_0 c}=c\varepsilon_0E^2$ (using $1/(\mu_0c)=c\varepsilon_0$). Averaging $E^2=E_0^2\cos^2(kz-\omega t)$ over a period gives $\langle\cos^2\rangle=\tfrac12$, hence

This matches refractive_index(2.25)$=1.5$, phase_velocity(2.25)$=$C$/1.5$, and wavelength(...)$=3.333\times10^{-7}$ m $=333.3$ nm.

With $n_1=1$, $n_2=1.5$ the amplitude coefficients are

The power fractions are $R=r^2=0.04$ and $T=\dfrac{n_2}{n_1}t^2=1.5\times0.64=0.96$, so $R+T=0.04+0.96=1$ — energy is conserved with no absorption. The sign $r<0$ (because $n_2>n_1$) means the reflected $\mathbf E$ is inverted: a $\pi$ phase shift on reflecting off the denser medium. Reversing the trip, glass→air, $r=\dfrac{1.5-1}{2.5}=+0.2$, so $R=(0.2)^2=0.04$ — identical, since $R=|r|^2$ is insensitive to the sign flip. These are exactly fresnel_normal(1,1.5)$=(-0.2,0.8)$, reflectance/transmittance$=0.04/0.96$, and reflectance(1.5,1.0)$=0.04$.

Separation of variables on the $B_z$ Helmholtz equation gives a transverse profile $\cos(k_x x)\cos(k_y y)$; the wall conditions quantize $k_x=m\pi/a$, $k_y=n\pi/b$. Propagation needs real $k$ in $k=\sqrt{\omega^2/c^2-k_x^2-k_y^2}$, and $k=0$ marks the cutoff:

With $a\ge b$, the smallest nonzero cutoff is $(m,n)=(1,0)$, giving $\omega_{10}=c\pi/a$ and $f_{10}=\omega_{10}/2\pi=c/2a$. For WR-90, $f_{10}=\dfrac{2.998\times10^8}{2(0.0229)}=6.55\text{ GHz}$ — exactly what cutoff_frequency(1,0,a,b) and dominant_mode_cutoff(a,b)/2\pi return ($6.546$ GHz, $\omega_{10}=4.113\times10^{10}$ rad/s).

With $a=2b$ substitute into $\omega_{mn}=c\pi\sqrt{(m/a)^2+(n/b)^2}$:

So nothing else turns on between $\omega_{10}$ and $2\omega_{10}$: only $\mathrm{TE}_{10}$ carries power across that octave, and at $2\omega_{10}$ the degenerate $\mathrm{TE}_{20}$ and $\mathrm{TE}_{01}$ switch on together. Numerically (with $a=2b$) the code gives $\omega_{20}/\omega_{10}=\omega_{01}/\omega_{10}=2.0$, and at $\omega=1.5\,\omega_{10}$ is_propagating is True for $(1,0)$ but False for $(2,0)$ and $(0,1)$ — confirming the single-mode band.

The longitudinal wavenumber follows from $\omega^2/c^2=k^2+(\omega_{co}/c)^2$, i.e. $k=\frac1c\sqrt{\omega^2-\omega_{co}^2}$. For $\omega>\omega_{co}$ this is real and the mode propagates as $e^{i(kz-\omega t)}$. For $\omega<\omega_{co}$ the radicand goes negative, $k=i\kappa$ with

so $e^{ikz}=e^{-\kappa z}$ — pure exponential decay, no oscillation and (since $\kappa$ is real, not complex) no Joule loss. At $\omega=\tfrac12\omega_{co}$, $\kappa=\frac1c\sqrt{\omega_{co}^2-\tfrac14\omega_{co}^2}=\frac{\sqrt3}{2}\frac{\omega_{co}}{c}\approx0.866\,\omega_{co}/c$. The code confirms guide_wavenumber$=0$ below cutoff while evanescent_decay$=288.87$ m⁻¹ $=0.866\,\omega_{co}/c$ (for $\omega_{co}=10^{11}$ rad/s).

The phase velocity is $v_p=\omega/k=\omega\big/\frac1c\sqrt{\omega^2-\omega_{co}^2} =\dfrac{c}{\sqrt{1-(\omega_{co}/\omega)^2}}$, which exceeds $c$. Differentiating the dispersion $\omega^2=c^2k^2+\omega_{co}^2$ gives $2\omega\,d\omega=2c^2k\,dk$, so

Their product is $v_p v_g=c^2$ identically. The superluminal $v_p$ is the speed of a pure phase front, which carries no information; energy and signals travel at $v_g<c$, so relativity is safe. As $\omega\gg\omega_{co}$ the radical $\to1$ and both $\to c$. The code confirms $v_p v_g=8.988\times10^{16}=$C**2, and at $\omega=100\,\omega_{co}$ both velocities sit within $10^{-4}$ of C.

A TEM mode has $E_z=B_z=0$, so the transverse field obeys $\nabla\cdot\mathbf E_\perp=0$ and $\nabla\times\mathbf E_\perp=0$ in the cross-section — an electrostatic pattern, the gradient of a potential solving Laplace's equation there. Inside a single hollow conductor the wall is one connected equipotential, forcing the potential constant and $\mathbf E_\perp=0$: a hollow guide carries no TEM wave. A coax has two conductors at different potentials, so a nontrivial radial $\mathbf E$ (and azimuthal $\mathbf B$) lives between them. Setting $E_z=B_z=0$ also kills the Helmholtz constant, $\omega^2/c^2-k^2=0$, so

the free-space dispersion: every frequency, down to DC, travels at $c$. This is the opposite of P3's guided mode, where $\omega^2/c^2-k^2=k_c^2>0$ imposes a cutoff. Hence tem_line_speed()$=c= 2.998\times10^8$ m/s at all frequencies, whereas phase_velocity_guide/group_velocity_guide only approach $c$ in the $\omega\gg\omega_{co}$ limit.

The retarded time obeys the implicit condition $t_r=t-\eta(t_r)/c$. For a charge fixed at the origin the separation to the field point never changes, $\eta(t_r)=r$ (constant), so the condition is no longer implicit and solves at once:

the same delay for every $t$. Feeding this unchanging source into the retarded potential (Eq. 10.26), the integrand is just the static Coulomb density, so

The light-travel delay $r/c$ is physically real but invisible, because $\rho$ never changed. This is exactly what retarded_time (constant trajectory) gives, $t-r/c$, and retarded_potential_static returns $(K_E\,q/r,\;t-r/c)$.

The Liénard-Wiechert potentials (Eq. 10.46–10.47) are $V=\dfrac{1}{4\pi\varepsilon_0}\dfrac{q}{\eta\,(1-\hat{\boldsymbol\eta}\cdot\mathbf v/c)}$ and $\mathbf A=(\mathbf v/c^2)V$. Setting $\mathbf v=\mathbf 0$ collapses the beaming denominator,

so

A charge at rest carries no vector potential, and $V$ is the ordinary Coulomb value at the retarded separation $\eta$ — the static limit recovered. With a constant zero-velocity trajectory lienard_wiechert returns $V=K_E\,q/\eta$, $|\mathbf A|=0$, and $t_r=t-\eta/c$.

With $\mathbf v=0.6c\,\hat{\mathbf x}$ aimed straight at a field point on the $+\hat x$ axis, the retarded separation points the same way, $\hat{\boldsymbol\eta}=\hat{\mathbf x}$, so $\hat{\boldsymbol\eta}\cdot\mathbf v/c=0.6$. The Liénard-Wiechert $V$ (Eq. 10.46) is the Coulomb value times the beaming factor:

Because $\mathbf v\neq\mathbf 0$, the vector potential $\mathbf A=(\mathbf v/c^2)V\neq\mathbf 0$. Running lienard_wiechert for a head-on $0.6c$ charge returns a potential exceeding the instantaneous $K_E\,q/\eta$ by a ratio of exactly $2.5000$, with $|\mathbf A|=3.60\times10^{-9}>0$.

Larmor's formula $P=\mu_0 q^2 a^2/6\pi c$ is quadratic in $a$, so scaling $a\to 2a$ multiplies the power by $2^2=4$, while $a=0$ (uniform motion) gives $P=0$ — it is acceleration, not speed, that makes a charge radiate. For an electron at $a=10^{22}\ \mathrm{m/s^2}$,

of order $10^{-10}$ W. This matches larmor_power(q,2e22)/larmor_power(q,1e22)$=4$ and larmor_power(1.602e-19,1e22)$\approx5.71\times10^{-10}$ W.

(i) The pattern carries the factor $\sin^2\theta$: it vanishes on the axis ($\sin 0=0$, an exact null along the dipole) and peaks broadside ($\sin\tfrac{\pi}{2}=1$) — the doughnut. (ii) Integrate over the sphere, pulling out the constant prefactor:

(iii) The total scales as $\omega^4$ (Rayleigh's law — high frequencies dominate, so the sky is blue) and as $p_0^2$: doubling $\omega$ multiplies $P$ by $2^4=16$, doubling $p_0$ by $2^2=4$. Numerically dipole_angular_power(0,...)$=0$, the broadside value is largest, total_power_from_pattern matches dipole_radiated_power to a relative $\sim10^{-13}$ ($\ll0.1\%$), and the $\times16$ and $\times4$ scalings hold exactly.

Filling in $F^{0i}=E_i/c$ for the top row/column and $F^{ij}=-\varepsilon_{ijk}B_k$ for the spatial block,

Antisymmetry $F^{\mu\nu}=-F^{\nu\mu}$ forces each diagonal entry to satisfy $F^{\mu\mu}=-F^{\mu\mu}=0$, so the four diagonal slots vanish; the remaining $16-4=12$ entries are locked in $\pm$ pairs, leaving $6$ independent numbers — the three $E_i/c$ (top row) and the three $B_k$ (spatial block). The spatial $3\times3$ block holds $F^{xy}=B_z$, $F^{xz}=-B_y$, $F^{yz}=B_x$, so $B_x$ lives in the $(y,z)$ slot — row 2, col 3 in the $(ct,x,y,z)$ ordering, i.e. F[2][3]. The code confirms is_antisymmetric $=$ True, fields_from_tensor round-trips (E, B), and F[2][3] $=B_x$.

Insert $E_x=E_z=0$, $E_y=E_0$, $\mathbf B=0$ into the boost map (Eq. 12.109). The electric components give $E'_x=E_x=0$, $E'_y=\gamma(E_y-vB_z)=\gamma E_0$, $E'_z=\gamma(E_z+vB_y)=0$, so $\mathbf E'=(0,\gamma E_0,0)$. The magnetic components give $B'_x=B_x=0$, $B'_y=\gamma(B_y+\tfrac{v}{c^2}E_z)=0$, and

With $\mathbf v=v\hat{\mathbf x}$ and $\mathbf E'=\gamma E_0\hat{\mathbf y}$, $\;-\tfrac1{c^2}\mathbf v\times\mathbf E'=-\tfrac{v\gamma E_0}{c^2}(\hat{\mathbf x}\times\hat{\mathbf y})=-\gamma\tfrac v{c^2}E_0\,\hat{\mathbf z}$, matching $\mathbf B'$. As $v\to0$ the magnetic field vanishes while $E'_y\to E_0$. For $E_0=1000$ V/m, $\beta=0.6$ ($\gamma=1.25$): $E'_y=1250$ V/m and $B'_z=-2.502\times10^{-6}$ T, exactly what boost_fields((0, E0, 0), (0,0,0), beta) returns.

Take $\mathbf E'\cdot\mathbf B'=E'_xB'_x+E'_yB'_y+E'_zB'_z$. The $x$-terms are untouched ($E'_xB'_x=E_xB_x$); the transverse pieces are

Expanding, the genuinely mixed terms $\tfrac v{c^2}E_yE_z$ and $-\tfrac v{c^2}E_zE_y$ cancel, as do $-vB_zB_y$ and $+vB_yB_z$, leaving $\gamma^2\big(1-\tfrac{v^2}{c^2}\big)(E_yB_y+E_zB_z)=E_yB_y+E_zB_z$ since $\gamma^2(1-v^2/c^2)=1$. Adding back $E_xB_x$ recovers $\mathbf E\cdot\mathbf B$. The identical cancellation in $B'^2-E'^2/c^2$ leaves it invariant. Numerically field_invariants returns $\mathbf E\cdot\mathbf B=-2.9\times10^{-3}$ and $B^2-E^2/c^2=9.19\times10^{-14}$, unchanged through $\beta=0.2,0.5,0.9,0.99$ — matching field_invariants(*boost_fields(E, B, beta)).

For the wave, $\mathbf E\perp\mathbf B$ gives $\mathbf E\cdot\mathbf B=0$, and $E=cB$ gives

so both invariants are zero. By P3 they are zero in every frame: $\mathbf E'\cdot\mathbf B'=0$ keeps $\mathbf E'\perp\mathbf B'$, and $B'^2-E'^2/c^2=0$ keeps $|\mathbf E'|=c|\mathbf B'|$. A "purely electric" frame ($\mathbf B'=0$) would force $B'^2-E'^2/c^2=-E'^2/c^2<0$, and a "purely magnetic" frame ($\mathbf E'=0$) would force $+B'^2>0$ — both contradict the vanishing invariant. So no boost can strip a light wave to one field. With $\mathbf E=(0,E_0,0)$, $\mathbf B=(0,0,E_0/c)$, field_invariants returns $\approx(0,0)$ and stays $\approx(0,0)$ after boost_fields.

With $E_x=E_z=0,\ E_y=E_0$ and $B_x=B_z=0,\ B_y=B_0$, the boost map gives

since $E'_z=\gamma(E_z+vB_y)=\gamma vB_0$ and $B'_z=\gamma(B_z-\tfrac v{c^2}E_y)=-\gamma\tfrac v{c^2}E_0$. Their dot product is

Because $\mathbf E'\cdot\mathbf B'=E_0B_0\neq0$ in every frame, neither field can vanish (a zero field would make the product zero). For $E_0=1000$ V/m, $B_0=2\times10^{-6}$ T, field_invariants gives $\mathbf E\cdot\mathbf B=2.0\times10^{-3}=E_0B_0$, and scanning $\beta=0.3,0.6,0.9$ with boost_fields holds $\mathbf E'\cdot\mathbf B'=E_0B_0$ while $|\mathbf E'|,|\mathbf B'|$ stay strictly positive.