physicskit.fields #

Returns:

ndarray, shape (…, 3) – Magnetic field in tesla.

Examples

A long straight wire recovers \(\mu_0 I / (2\pi s)\):

>>> wire = [[0, 0, -1e4], [0, 0, 1e4]]
>>> B = biot_savart_field(wire, [[0.1, 0.0, 0.0]], current=1.0)
>>> round(float(B[0, 1] / infinite_wire_field(1.0, 0.1)), 6)
1.0

physicskit.fields.casimir_energy_1d(d, c=1.0, cutoff=None)[source]#

Regularized zero-point (Casimir) energy of the discrete 1D cavity spectrum.

Uses the standard exponential-cutoff regularization: the sum \(\sum_n n x^n = x/(1-x)^2\) (with \(x=e^{-a}\), \(a=\pi c\,\text{cutoff}/d\)) has an exact closed form, so no series truncation is needed. Its small-cutoff (small-a) expansion is \(\tfrac{1}{4\sinh^2(a/2)} = 1/a^2 - 1/12 + O(a^2)\); the \(1/a^2\) piece is the (unphysical, cutoff-scheme-dependent) divergence that a continuum reference calculation regularized the same way would also produce, so it is subtracted exactly, leaving the finite remainder that survives as cutoff -> 0:

\[E(d) \to -\frac{\pi c}{24 d}\]

the standard 1D massless-field Casimir energy (equivalently, the Casimir energy of a CFT strip with central charge 1).

Parameters:

d (float) – Plate separation.
c (float) – Wave speed.
cutoff (float | None) – Regulator scale. Defaults to 1e-3 * d / c. This subtraction is a difference of two large, nearly-equal floating-point terms (\(\propto 1/\text{cutoff}^2\)), so making cutoff too small loses precision to cancellation rather than gaining it; 1e-4*d/c to 1e-2*d/c is the well-behaved range in double precision.

Return type:

Returns:

float – The regularized (finite, cutoff-subtracted) zero-point energy.

See also

casimir_mode_frequencies: The (bare, un-regularized) mode spectrum being summed.

Examples

The cutoff-and-subtract result agrees with the known closed form \(-\pi c/(24d)\), and (within the well-behaved cutoff range noted above) agrees more closely as the cutoff shrinks:

>>> d = 3.0
>>> exact = -np.pi * 1.0 / (24 * d)
>>> bool(abs(casimir_energy_1d(d, cutoff=3e-3 * d) - exact) < 1e-6)
True
>>> bool(abs(casimir_energy_1d(d, cutoff=7e-4 * d) - exact) < 1e-8)
True

The energy grows less negative (weaker confinement of vacuum energy) as the plates separate, giving an attractive force \(-dE/dd < 0\):

>>> bool(casimir_energy_1d(1.0) < casimir_energy_1d(2.0) < 0)
True

physicskit.fields.casimir_mode_frequencies(d, c=1.0, n_max=200)[source]#

Discrete standing-wave mode frequencies of a 1D cavity of plate separation d.

Parameters:

d (float) – Plate separation.
c (float) – Wave speed (=1 in natural units; use the physical speed of light for SI-unit frequencies).
n_max (int) – Number of modes to return.

Return type:

Returns:

ndarray, shape (n_max,) – \(\omega_n = n\pi c/d\) for \(n=1,\dots,n_{max}\).

See also

casimir_energy_1d: The regularized zero-point energy of this mode spectrum.

Examples

>>> omega = casimir_mode_frequencies(d=2.0, c=1.0, n_max=3)
>>> [round(float(w), 4) for w in omega]
[1.5708, 3.1416, 4.7124]

physicskit.fields.circular_loop_path(radius, n_segments=256, center=(0.0, 0.0, 0.0))[source]#

Vertices of a closed circular loop in a plane of constant z.

The loop runs counter-clockwise seen from +z, so a positive current produces a field along +z at its center.

Parameters:

radius (float) – Loop radius in meters.
n_segments (int) – Number of straight segments; the polygon’s on-axis field differs from a true circle’s by \(O(1/n^2)\).
center (array_like, shape (3,), default=(0, 0, 0)) – Loop center.

Return type:

Returns:

ndarray, shape (n_segments + 1, 3) – Vertices, with the first repeated at the end.

Examples

>>> path = circular_loop_path(1.0, n_segments=4)
>>> path.shape, bool(np.allclose(path[0], path[-1]))
((5, 3), True)

physicskit.fields.coulomb_field(charges, positions, points, eps=8.8541878188e-12)[source]#

Electric field of a set of point charges (Coulomb’s inverse-square law).

\[\mathbf E(\mathbf r) = \frac{1}{4\pi\varepsilon}\sum_i q_i \frac{\mathbf r - \mathbf r_i}{|\mathbf r - \mathbf r_i|^3}\]

(Jackson eq. 1.5).

Parameters:

charges (array_like, shape (n,)) – Point charges in coulombs.
positions (array_like, shape (n, d)) – Charge positions in meters.
points (array_like, shape (..., d)) – Observation points.
eps (float) – Permittivity of the surrounding medium.

Return type:

Returns:

ndarray, shape (…, d) – Electric field in V/m; nan exactly at a charge.

Examples

>>> E = coulomb_field([1e-9], [[0.0, 0.0, 0.0]], [[2.0, 0.0, 0.0]])
>>> round(float(E[0, 0]), 4)  # a quarter of the field at 1 m
2.2469

physicskit.fields.coulomb_potential(charges, positions, points, eps=8.8541878188e-12)[source]#

Electrostatic potential of a set of point charges (Coulomb’s law plus superposition).

\[\phi(\mathbf r) = \frac{1}{4\pi\varepsilon}\sum_i \frac{q_i}{|\mathbf r - \mathbf r_i|}\]

(Jackson, Classical Electrodynamics, 3rd ed., eq. 1.17).

Parameters:

charges (array_like, shape (n,)) – Point charges in coulombs.
positions (array_like, shape (n, d)) – Charge positions in meters, d = 2 or 3 (2D positions are points in the plane z = 0 of ordinary 3D space).
points (array_like, shape (..., d)) – Observation points.
eps (float) – Permittivity of the surrounding medium.

Return type:

Returns:

ndarray, shape (…) – Potential in volts; inf exactly at a charge.

Examples

>>> phi = coulomb_potential([1e-9], [[0.0, 0.0, 0.0]], [[1.0, 0.0, 0.0]])
>>> round(float(phi[0]), 3)  # 1 nC at 1 m: about 8.988 V
8.988

physicskit.fields.count_vortices(psi, density_threshold=0.05)[source]#

Locate quantized vortices by summing the phase winding around each grid plaquette.

Plaquettes in very low density regions are excluded, since the phase of a near-zero-amplitude wavefunction is dominated by numerical noise and produces spurious windings there.

A vortex core that sits exactly on a grid vertex (rather than strictly inside a plaquette) can go undetected, since its circulation is then split ambiguously between the four plaquettes touching that vertex; in practice this is avoided by seeding vortices at positions not exactly on the grid.

Parameters:

psi (ndarray) – Wavefunction.
density_threshold (float) – Plaquettes where \(|\psi|^2\) (relative to its maximum) falls below this fraction are excluded from the search.

Return type:

Returns:

ndarray of int, shape (n, n) – Integer winding number around each plaquette’s corner (i, j); nonzero entries mark a vortex core (+1 or -1 for a singly-quantized vortex/antivortex) enclosed by that plaquette.

Examples

>>> import numpy as np
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(48, 10.0)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi_vortex = gpe_imprint_vortex(psi0, X, Y, [(0.0, 0.0)])
>>> winding = count_vortices(psi_vortex)
>>> int(np.sum(np.abs(winding)))
1

physicskit.fields.courant_limit_1d(dx)[source]#

Maximum stable time step for the 1D FDTD update.

Parameters:: dx (float) – Spatial grid spacing in meters.
Return type:: float
Returns:: float – The Courant stability limit \(\Delta t = \Delta x / c_0\).

Examples

>>> round(float(courant_limit_1d(1e-3) * C0 / 1e-3), 6)
1.0

physicskit.fields.courant_limit_2d(dx, dy)[source]#

Maximum stable time step for the 2D FDTD update.

Parameters:

dx (float) – Spatial grid spacing in meters along each axis.
dy (float) – Spatial grid spacing in meters along each axis.

Return type:

Returns:

float – The Courant stability limit \(\Delta t = 1 / (c_0 \sqrt{1/\Delta x^2 + 1/\Delta y^2})\).

Examples

>>> dt = courant_limit_2d(1e-3, 1e-3)
>>> round(float(dt * C0 / 1e-3), 6)
0.707107

physicskit.fields.dielectric_slab(shape, i_start, i_end, eps_r_slab=4.0)[source]#

A planar dielectric slab: eps_r_slab between grid columns i_start and i_end, vacuum elsewhere.

A minimal demo permittivity map for fdtd_2d_tmz() / fdtd_2d_tmz_evolve(): a wave launched from one side crossing i_start partially reflects and partially refracts/transmits into the slab, and again at i_end on the way out, illustrating propagation through a medium interface.

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
i_start (int) – Grid-index bounds of the slab along the x axis (rows i_start:i_end).
i_end (int) – Grid-index bounds of the slab along the x axis (rows i_start:i_end).
eps_r_slab (float) – Relative permittivity inside the slab.

Return type:

Returns:

ndarray, shape (Nx, Ny) – eps_r map: eps_r_slab inside the slab, 1.0 outside.

Examples

>>> eps_r = dielectric_slab((40, 20), i_start=15, i_end=25, eps_r_slab=4.0)
>>> float(eps_r[10, 5]), float(eps_r[20, 5])
(1.0, 4.0)

physicskit.fields.dipole_field(p, points, origin=(0.0, 0.0, 0.0), eps=8.8541878188e-12)[source]#

Electric field of an ideal point dipole.

\[\mathbf E(\mathbf x) = \frac{1}{4\pi\varepsilon} \frac{3\hat{\mathbf n}(\mathbf p\cdot\hat{\mathbf n}) - \mathbf p}{r^3}\]

(Jackson eq. 4.13, away from the origin). On the dipole axis this is \(2p/(4\pi\varepsilon r^3)\), twice the equatorial magnitude.

Parameters:

p (array_like, shape (3,)) – Dipole moment in C m.
points (array_like, shape (..., 3)) – Observation points.
origin (array_like, shape (3,), default=(0, 0, 0)) – Dipole location.
eps (float) – Permittivity of the medium.

Return type:

Returns:

ndarray, shape (…, 3) – Field in V/m.

Examples

>>> E = dipole_field([0, 0, 1.0], [[0, 0, 1.0], [1.0, 0, 0]], eps=1 / (4 * np.pi))
>>> E.round(12).tolist()
[[0.0, 0.0, 2.0], [0.0, 0.0, -1.0]]

physicskit.fields.electric_field_from_potential(phi, spacing)[source]#

Electric field \(\mathbf E = -\nabla\phi\) from a gridded potential.

Uses second-order central differences in the interior (one-sided at the faces), via numpy.gradient().

Parameters:

phi (array_like) – Potential on the grid, in volts.
spacing (float or sequence of float) – Grid spacing in meters, one value or one per axis.

Return type:

Returns:

ndarray, shape (phi.ndim, *phi.shape) – Field components E[0] (along axis 0), E[1], … in V/m.

Examples

>>> x = np.linspace(0.0, 1.0, 11)
>>> E = electric_field_from_potential(np.outer(-5.0 * x, np.ones(4)), 0.1)
>>> bool(np.allclose(E[0], 5.0)), bool(np.allclose(E[1], 0.0))
(True, True)

physicskit.fields.fdtd_1d(Ez0, Hy0, eps_r, mu_r, steps, dt, dx, sigma=None)[source]#

Evolve a 1D plane wave (Ez, Hy) on a Yee grid using leapfrog FDTD.

Discretizes Maxwell’s curl equations for fields varying only along x, \(\partial_t H_y = \partial_x E_z/\mu\) and \(\partial_t E_z = \partial_x H_y/\epsilon\) (the same signs as fdtd_2d_tmz()), so a wave moving toward +x has \(H_y = -E_z/\eta\) (Poynting vector \(-E_zH_y>0\)).

Parameters:

Ez0 (ndarray) – Initial electric field.
Hy0 (ndarray) – Initial magnetic field, staggered half a cell ahead of Ez0.
eps_r (ndarray) – Relative permittivity and permeability at each Ez grid point.
mu_r (ndarray) – Relative permittivity and permeability at each Ez grid point.
steps (int) – Number of time steps to advance.
dt (float) – Time step in seconds; must satisfy courant_limit_1d().
dx (float) – Spatial grid spacing in meters.
sigma (ndarray | None) – Electric conductivity (S/m) at each grid point, e.g. from pml_conductivity_profile(), for a lossy/absorbing boundary. Defaults to a lossless grid with hard (PEC) walls at both ends.

Return type:

Returns:

Ez, Hy (ndarray) – Fields at the final time step, same shapes as Ez0, Hy0.

See also

fdtd_2d_tmz: The 2D (TMz-mode) analog of this solver.

Examples

A smooth Gaussian pulse, launched with the impedance-matched initial condition of a purely right-moving wave (\(H_y=-E_z/\eta_0\)), propagates at exactly the vacuum speed of light:

>>> import numpy as np
>>> N = 800
>>> dx = 1e-3
>>> dt = 0.99 * courant_limit_1d(dx)
>>> eta0 = np.sqrt(MU0 / EPS0)
>>> x0, sigma_pulse = 100, 25
>>> Ez0 = np.exp(-((np.arange(N) - x0) ** 2) / (2 * sigma_pulse ** 2))
>>> xh = np.arange(N - 1) + 0.5
>>> Hy0 = -np.exp(-((xh - x0) ** 2) / (2 * sigma_pulse ** 2)) / eta0
>>> eps_r, mu_r = np.ones(N), np.ones(N)
>>> Ez, Hy = fdtd_1d(Ez0, Hy0, eps_r, mu_r, steps=360, dt=dt, dx=dx)
>>> peak = np.argmax(Ez)
>>> measured_speed = (peak - x0) * dx / (360 * dt)
>>> round(float(measured_speed / C0), 3)
1.002

physicskit.fields.fdtd_2d_tmz(Ez0, Hx0, Hy0, eps_r, mu_r, steps, dt, dx, dy)[source]#

Evolve a 2D TMz-mode field (Ez, Hx, Hy) on a Yee grid using leapfrog FDTD.

Parameters:

Ez0 (ndarray) – Initial (out-of-plane) electric field.
Hx0 (ndarray) – Initial in-plane magnetic field components.
Hy0 (ndarray) – Initial in-plane magnetic field components.
eps_r (ndarray) – Relative permittivity and permeability maps.
mu_r (ndarray) – Relative permittivity and permeability maps.
steps (int) – Number of time steps to advance.
dt (float) – Time step in seconds; must satisfy courant_limit_2d().
dx (float) – Spatial grid spacing in meters along each axis.
dy (float) – Spatial grid spacing in meters along each axis.

Return type:

Returns:

Ez, Hx, Hy (ndarray, shape (Nx, Ny)) – Fields at the final time step.

See also

fdtd_1d: The 1D analog of this solver.
poynting_vector_tmz: Energy flux from the evolved fields.

Examples

>>> import numpy as np
>>> Nx, Ny = 50, 50
>>> Ez0 = np.zeros((Nx, Ny))
>>> Hx0 = np.zeros((Nx, Ny))
>>> Hy0 = np.zeros((Nx, Ny))
>>> eps_r, mu_r = np.ones((Nx, Ny)), np.ones((Nx, Ny))
>>> dx = dy = 1e-3
>>> dt = 0.5 * courant_limit_2d(dx, dy)
>>> Ez0[Nx // 2, Ny // 2] = 1.0
>>> Ez, Hx, Hy = fdtd_2d_tmz(Ez0, Hx0, Hy0, eps_r, mu_r, steps=10, dt=dt, dx=dx, dy=dy)
>>> Ez.shape
(50, 50)

physicskit.fields.fdtd_2d_tmz_evolve(Ez0, Hx0, Hy0, eps_r, mu_r, steps, dt, dx, dy, sigma=None, source=None, snapshot_stride=1)[source]#

Time-step the 2D TMz FDTD update while recording Ez snapshots, for animation.

Unlike fdtd_2d_tmz(), which returns only the field at the final step, this driver records Ez every snapshot_stride steps and supports two features layered on top of the same base update: a lossy absorbing boundary (sigma, e.g. from pml_conductivity_profile_2d(), so an outgoing wave is damped instead of reflecting off the hard PEC walls) and a soft source callback injected into Ez after every field update (e.g. oscillating_dipole_source()).

Parameters:

Ez0 (ndarray) – As in fdtd_2d_tmz().
Hx0 (ndarray) – As in fdtd_2d_tmz().
Hy0 (ndarray) – As in fdtd_2d_tmz().
eps_r (ndarray) – As in fdtd_2d_tmz().
mu_r (ndarray) – As in fdtd_2d_tmz().
steps (int) – As in fdtd_2d_tmz().
dt (float) – As in fdtd_2d_tmz().
dx (float) – As in fdtd_2d_tmz().
dy (float) – As in fdtd_2d_tmz().
sigma (ndarray | None) – Electric conductivity map. None (default) reproduces fdtd_2d_tmz()’s lossless, hard-wall behavior exactly.
source (callable, optional) – source(Ez, t), called after each field update and expected to mutate Ez in place (a soft source: added to, not overwriting, the field already there).
snapshot_stride (int) – Record Ez every this many steps (plus the initial condition).

Return type:

Returns:

snapshots (ndarray, shape (n_recorded, Nx, Ny)) – Ez at t=0 and after every recorded step.
times (ndarray, shape (n_recorded,)) – Time of each recorded snapshot.

See also

fdtd_2d_tmz: The single-shot (final-state-only) solver this wraps.

Examples

>>> import numpy as np
>>> Nx, Ny = 40, 40
>>> Ez0 = np.zeros((Nx, Ny))
>>> Hx0 = Hy0 = np.zeros((Nx, Ny))
>>> eps_r = mu_r = np.ones((Nx, Ny))
>>> dx = dy = 1e-3
>>> dt = 0.5 * courant_limit_2d(dx, dy)
>>> source = oscillating_dipole_source(Nx // 2, Ny // 2, amplitude=1.0, freq=5e10)
>>> snaps, times = fdtd_2d_tmz_evolve(Ez0, Hx0, Hy0, eps_r, mu_r, steps=20, dt=dt, dx=dx, dy=dy, source=source, snapshot_stride=5)
>>> snaps.shape
(5, 40, 40)

physicskit.fields.flux_tube_energy_density_2d(shape, dx, dy, separation, flux_quantum=1.0, tube_width=None)[source]#

2D energy-density map of the confined flux tube, for visualizing it as a stretching heatmap.

Extrudes flux_tube_field_1d() (evaluated along the axis joining the two charges) into 2D by giving it a fixed-width Gaussian transverse profile – the “fixed cross-section” that produces confinement – and forms the field energy density \(\tfrac{1}{2}\,\text{field}(x)^2\times\text{transverse}(y)^2\).

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
dx (float) – Grid spacing; the charges are placed on the centered grid at \(x=\pm\text{separation}/2\), \(y=0\).
dy (float) – Grid spacing; the charges are placed on the centered grid at \(x=\pm\text{separation}/2\), \(y=0\).
separation (float) – Distance between the two charges.
flux_quantum (float) – Charge magnitude, as in flux_tube_field_1d().
tube_width (float | None) – Transverse (y) width of the flux tube. Defaults to 5 * dy.

Return type:

Returns:

ndarray, shape (Nx, Ny) – Energy density, largest along the tube connecting the charges and decaying away from its axis.

Examples

>>> e = flux_tube_energy_density_2d((200, 40), dx=0.2, dy=0.2, separation=10.0)
>>> bool(e.sum() > 0)
True

physicskit.fields.flux_tube_field_1d(x, separation, flux_quantum=1.0, wall_width=None)[source]#

Static 1D confined flux profile between two opposite point charges (toy confinement model).

Solves the 1D Gauss law \(d(\text{flux})/dx = \rho(x)\) for a field confined to a fixed-cross-section tube connecting two opposite charges at \(x=\pm\text{separation}/2\), each smoothed into a Gaussian of width wall_width (a numerical regularization of the point charge, not a physical effect). The field is computed by direct numerical integration (cumulative trapezoidal sum) of this charge density, so it is genuinely “solved for” at each separation rather than assumed.

Parameters:

x (ndarray) – 1D grid the charges live on.
separation (float) – Distance between the two charges.
flux_quantum (float) – Charge magnitude (equivalently, the plateau field strength inside the tube).
wall_width (float | None) – Smoothing width of each point charge. Defaults to 3 * (x[1] - x[0]), i.e. a few grid cells – small enough to look like a point charge but resolved on the grid.

Return type:

Returns:

ndarray, same shape as x – The confined field, approximately +flux_quantum between the charges and 0 outside, with smooth transitions of width wall_width.

See also

flux_tube_energy_density_2d: The corresponding 2D energy-density map, for animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-20, 20, 2000)
>>> field = flux_tube_field_1d(x, separation=10.0, flux_quantum=1.0)
>>> bool(abs(field[np.argmin(np.abs(x))] - 1.0) < 0.05)
True
>>> bool(abs(field[0]) < 0.05)
True

physicskit.fields.gpe_energy(psi, V, g, X, Y, K2)[source]#

Evaluate the energy and angular momentum of a GPE wavefunction.

Parameters:

psi (ndarray) – Wavefunction.
V (ndarray) – External trapping potential.
g (float) – Interaction strength.
X (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
Y (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
K2 (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).

Return type:

dict

Returns:

dict – {"kinetic", "potential", "interaction", "angular_momentum", "total"}, where total is the lab-frame energy (kinetic + potential + interaction) and angular_momentum is \(\langle L_z \rangle\). The rotating-frame energy at rotation rate Omega is total - Omega * angular_momentum.

Examples

>>> import numpy as np
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(48, 10.0)
>>> V = 0.5 * (X ** 2 + Y ** 2)
>>> psi = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi *= 1.0 / np.sqrt(np.sum(np.abs(psi) ** 2) * (X[1, 0] - X[0, 0]) ** 2)
>>> E = gpe_energy(psi, V, g=0.0, X=X, Y=Y, K2=K2)
>>> round(float(E["kinetic"] + E["potential"]), 4)
1.0

physicskit.fields.gpe_evolve(psi0, V, g, dt, steps, K2, snapshot_stride=1)[source]#

Real-time propagation of the 2D Gross-Pitaevskii / cubic NLS equation via split-step Fourier.

Solves \(i\partial_t\psi = [-\tfrac{1}{2}\nabla^2 + V + g|\psi|^2]\psi\) forward in real time (unlike gpe_relax()’s imaginary-time relaxation, which only finds stationary states): the kinetic term is exact in Fourier space, and the potential-plus-nonlinear term is exact as a pointwise phase rotation in real space, the same split-step structure as physicskit.fields.solitons.nls_evolve() generalized to 2D and to an external potential. Real-time evolution is unitary and conserves the norm on its own (no renormalization needed, unlike imaginary time).

This one stepper serves two purposes, distinguished only by V and the sign of g:

With a harmonic trap (V = 0.5*(X**2+Y**2)) and repulsive interactions (g > 0), an off-center vortex genuinely precesses around the trap under the density gradient – real vortex dynamics, as opposed to the static relaxed state gpe_relax() finds.
With V = 0 (or weak) and attractive interactions (g < 0 here, opposite sign convention from nls_evolve’s g since this module’s equation carries a +g|psi|^2 term), a sufficiently tall, narrow initial packet self-focuses, concentrating into a narrower, taller peak – a numerical stand-in for wave collapse. True collapse is a singularity in finite time; this integrator (like any finite-grid scheme) cannot resolve it and the calculation should be stopped once the peak density is still visibly growing, not carried through the blow-up itself.

Parameters:

psi0 (ndarray) – Initial wavefunction.
V (ndarray) – External potential (use np.zeros_like for the free/collapse case).
g (float) – Interaction strength and sign (see above).
dt (float) – Time step.
steps (int) – Number of steps to advance.
K2 (ndarray) – Wavenumber-squared grid from harmonic_trap_grid().
snapshot_stride (int) – Record a snapshot every this many steps (plus the initial condition).

Return type:

Returns:

snapshots (ndarray of complex, shape (n_recorded, n, n)) – Wavefunction at t=0 and after every recorded step.
times (ndarray, shape (n_recorded,)) – Time of each recorded snapshot.

See also

gpe_relax: Imaginary-time relaxation to a stationary state.
physicskit.fields.solitons.nls_evolve: The 1D analog (opposite sign convention for g).

Examples

Norm is conserved under real-time evolution, unlike the imaginary-time propagation in gpe_relax() (which requires explicit renormalization):

>>> import numpy as np
>>> n, length = 48, 12.0
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(n, length)
>>> dxg = X[1, 0] - X[0, 0]
>>> V = 0.5 * (X ** 2 + Y ** 2)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi0 *= 1.0 / np.sqrt(np.sum(np.abs(psi0) ** 2) * dxg * dxg)
>>> snaps, times = gpe_evolve(psi0, V, g=2.0, dt=1e-3, steps=200, K2=K2, snapshot_stride=50)
>>> norms = np.sum(np.abs(snaps) ** 2, axis=(1, 2)) * dxg * dxg
>>> bool(np.max(np.abs(norms - norms[0])) < 1e-6)
True

physicskit.fields.gpe_imprint_vortex(psi, X, Y, positions)[source]#

Imprint one singly-quantized vortex per position by multiplying in a phase winding.

Parameters:

psi (ndarray) – Wavefunction to imprint onto.
X (ndarray) – Real-space coordinate grids from harmonic_trap_grid().
Y (ndarray) – Real-space coordinate grids from harmonic_trap_grid().
positions (sequence of tuple(float, float)) – (x0, y0) center of each vortex to imprint.

Return type:

Returns:

ndarray of complex, shape (n, n) – psi multiplied by \(\prod_i [(x-x_{0,i}) + i(y-y_{0,i})]\), renormalized to preserve the input norm.

Examples

>>> import numpy as np
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(32, 12.0)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi = gpe_imprint_vortex(psi0, X, Y, [(0.0, 0.0)])
>>> bool(abs(psi[16, 16]) < 1e-10)
True

physicskit.fields.gpe_relax(psi0, V, g, dtau, steps, X, Y, K2, Omega=0.0, n_particles=1.0)[source]#

Relax a GPE initial condition toward a stationary state via imaginary-time propagation.

Parameters:

psi0 (ndarray) – Initial wavefunction.
V (ndarray) – External trapping potential.
g (float) – Interaction (nonlinearity) strength.
dtau (float) – Imaginary-time step.
steps (int) – Number of steps to advance.
X (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
Y (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
K2 (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
Omega (float) – Rotation frequency of the trap.
n_particles (float) – Total particle number; the wavefunction is renormalized to this value after every step (imaginary-time propagation does not conserve norm on its own).

Return type:

Returns:

ndarray of complex, shape (n, n) – The relaxed (approximately stationary) wavefunction.

See also

gpe_energy: Evaluate the energy of a relaxed state.
gpe_imprint_vortex: Seed an initial condition with quantized vortices.

Notes

Energy decreases monotonically under this propagation whenever it is implemented correctly (a basic sanity check worth verifying on any new potential or parameter regime).

Examples

Comparing the rotating-frame energy of the vortex-free ground state against a state seeded with one centered vortex reproduces the standard vortex-nucleation criterion: \(\Omega_c = \Delta E/\Delta L_z\), above which the vortex state has lower energy in the rotating frame:

>>> import numpy as np
>>> n, length, g = 64, 12.0, 4.0
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(n, length)
>>> V = 0.5 * (X ** 2 + Y ** 2)
>>> psi_vf = gpe_relax(np.exp(-0.5 * (X**2 + Y**2)).astype(complex), V, g, dtau=5e-4, steps=4000, X=X, Y=Y, K2=K2)
>>> psi0_v = gpe_imprint_vortex(np.exp(-0.5 * (X**2 + Y**2)).astype(complex), X, Y, [(0.0, 0.0)])
>>> psi_v = gpe_relax(psi0_v, V, g, dtau=5e-4, steps=4000, X=X, Y=Y, K2=K2)
>>> E0, E1 = gpe_energy(psi_vf, V, g, X, Y, K2), gpe_energy(psi_v, V, g, X, Y, K2)
>>> Omega_c = (E1["total"] - E0["total"]) / (E1["angular_momentum"] - E0["angular_momentum"])
>>> bool(0.0 < Omega_c < 1.0)
True

physicskit.fields.harmonic_trap_grid(n, length)[source]#

Build a centered real-space and wavenumber grid for a 2D harmonic trap.

Parameters:

n (int) – Number of grid points along each axis.
length (float) – Physical domain size (the domain is [-length/2, length/2)).

Return type:

Returns:

X, Y (ndarray, shape (n, n)) – Real-space coordinate grids, centered at the trap origin.
KX, KY (ndarray, shape (n, n)) – Wavenumber grids.
K2 (ndarray, shape (n, n)) – \(K_X^2 + K_Y^2\).

Examples

>>> X, Y, KX, KY, K2 = harmonic_trap_grid(32, 12.0)
>>> float(X[0, 0])
-6.0

physicskit.fields.image_charges_plane(q, position)[source]#

Image charge making the grounded conducting plane z = 0 an equipotential.

A charge \(q\) at \((x, y, d)\) above a grounded plane sees, in the region \(z > 0\), exactly the potential of itself plus an image \(-q\) at \((x, y, -d)\) (Jackson §2.1). The pair is Green’s function for the half-space with Dirichlet boundary.

Parameters:

q (float) – The real charge, in coulombs.
position (array_like, shape (3,)) – Its position (x, y, d) with d != 0.

Return type:

Returns:

charges (ndarray, shape (1,)) – The image charge [-q].
positions (ndarray, shape (1, 3)) – The image position [[x, y, -d]].

Examples

>>> qi, ri = image_charges_plane(1e-9, [0.0, 0.0, 0.5])
>>> phi = coulomb_potential(np.r_[1e-9, qi], np.vstack([[0, 0, 0.5], ri]), [[0.3, -0.2, 0.0]])
>>> abs(float(phi[0])) < 1e-12
True

physicskit.fields.image_charges_sphere(q, position, radius, center=(0.0, 0.0, 0.0), sphere_charge=None)[source]#

Image charge(s) for a point charge and a conducting sphere.

For a sphere of radius \(R\) and a charge \(q\) at distance \(a\) from its center, the image

\[q' = -q\frac{R}{a} \quad\text{at distance}\quad \frac{R^2}{a}\]

along the same ray makes the sphere’s surface an equipotential at zero potential (Jackson §2.2, eq. 2.4). This holds with the charge outside (\(a > R\), image inside) or inside a spherical cavity (\(a < R\), image outside). For an isolated sphere carrying total charge \(Q\) (charge outside only), a second image \(Q - q'\) at the center restores the net charge while keeping the surface an equipotential (Jackson §2.3).

Parameters:

q (float) – The real charge, in coulombs.
position (array_like, shape (3,)) – Its position, in meters.
radius (float) – Sphere radius R.
center (array_like, shape (3,), default=(0, 0, 0)) – Sphere center.
sphere_charge (float | None) – Total charge of an isolated sphere. None (default) means the sphere is grounded.

Return type:

Returns:

charges (ndarray, shape (1,) or (2,)) – Image charge(s).
positions (ndarray, shape (1, 3) or (2, 3)) – Image position(s).

Examples

>>> qi, ri = image_charges_sphere(1.0, [2.0, 0.0, 0.0], radius=1.0)
>>> float(qi[0]), float(ri[0, 0])
(-0.5, 0.5)

physicskit.fields.induced_charge_density_plane(q, d, rho)[source]#

Surface charge induced on the grounded plane z = 0 by a charge at height d.

\[\sigma(\varrho) = -\frac{q\, d}{2\pi(\varrho^2 + d^2)^{3/2}}\]

(Griffiths, Introduction to Electrodynamics, 4th ed., eq. 3.10), which integrates to exactly \(-q\) over the plane.

Parameters:

q (float) – The real charge, in coulombs.
d (float) – Its height above the plane, in meters.
rho (array_like) – Distance along the plane from the foot of the perpendicular.

Return type:

Returns:

ndarray – Surface charge density in C/m^2.

Examples

>>> round(float(induced_charge_density_plane(1.0, 1.0, 0.0) * 2 * np.pi), 12)
-1.0

physicskit.fields.induced_charge_density_sphere(q, a, radius, theta)[source]#

Surface charge induced on a grounded sphere by an outside point charge.

\[\sigma(\gamma) = -\frac{q\,(a^2 - R^2)}{4\pi R\,(R^2 + a^2 - 2aR\cos\gamma)^{3/2}}\]

(Jackson eq. 2.5, rewritten in the sphere radius \(R\) and charge distance \(a\)), where \(\gamma\) is the angle from the ray through the charge. It integrates to the image charge \(-qR/a\).

Parameters:

q (float) – The real charge, in coulombs.
a (float) – Distance of the charge from the center, a > radius.
radius (float) – Sphere radius R.
theta (array_like) – Polar angle(s) gamma on the sphere, measured from the charge’s direction.

Return type:

Returns:

ndarray – Surface charge density in C/m^2.

Examples

>>> from scipy.integrate import trapezoid
>>> th = np.linspace(0.0, np.pi, 20001)
>>> sigma = induced_charge_density_sphere(1.0, 2.0, 1.0, th)
>>> total = trapezoid(sigma * 2 * np.pi * np.sin(th), th)
>>> round(float(total), 6)  # the image charge -q R / a
-0.5

physicskit.fields.infinite_wire_field(current, s, mu=1.25663706127e-06)[source]#

Field magnitude at distance s from an infinite straight wire, \(\mu I / (2\pi s)\).

(Griffiths eq. 5.38.) The field circles the wire in the right-hand sense about the current.

Parameters:

current (float) – Current in amperes.
s (array_like) – Perpendicular distance from the wire, in meters.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray – \(|\mathbf B|\) in tesla.

Examples

>>> round(float(infinite_wire_field(1.0, 1.0)), 10)  # 2e-7 T per ampere at 1 m
2e-07

physicskit.fields.kdv_evolve(u0, x, dt, steps)[source]#

Evolve a KdV initial condition forward in time on a periodic domain.

Parameters:

u0 (ndarray) – Initial field, sampled on the periodic grid x.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps (int) – Number of steps to advance.

Return type:

Returns:

ndarray – Field after steps * dt time units.

Notes

Because KdV solitons are true solitons, two colliding solitons of different speed pass through each other elastically: each emerges from the collision with its original amplitude and speed, shifted only in phase – unlike generic nonlinear waves, which would break up or change shape on collision.

Examples

A fast, tall soliton launched behind a slower, shorter one overtakes it; both retain their original amplitudes after the collision:

>>> import numpy as np
>>> from scipy.signal import find_peaks
>>> N, L = 512, 60.0
>>> x = np.linspace(-L / 2, L / 2, N, endpoint=False)
>>> u0 = kdv_soliton(x, c=9.0, x0=-20) + kdv_soliton(x, c=4.0, x0=-8)
>>> u = kdv_evolve(u0, x, dt=0.0005, steps=6000)
>>> peaks, _ = find_peaks(u, height=1.0)
>>> heights = sorted(u[peaks])
>>> [round(float(h), 1) for h in heights]
[2.0, 4.5]

physicskit.fields.kdv_evolve_frames(u0, x, dt, steps_per_frame, n_frames)[source]#

Evolve a KdV field and record a snapshot every steps_per_frame steps, for animation.

A thin looping wrapper around kdv_evolve() – no new physics, just repeated short calls with each chunk’s output fed back in as the next chunk’s initial condition (valid since kdv_evolve() is a deterministic one-step-at-a-time integrator), so the trajectory can be animated frame by frame.

Parameters:

u0 (ndarray) – Initial field.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps_per_frame (int) – Number of integration steps between recorded frames.
n_frames (int) – Number of frames to record after the initial condition.

Return type:

Returns:

frames (ndarray, shape (n_frames + 1, len(x))) – Field at t=0 and after each recorded chunk.
times (ndarray, shape (n_frames + 1,))

See also

physicskit.fields.visualizers.animate_field_1d: Render these frames as a line-plot animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-30, 30, 256, endpoint=False)
>>> u0 = kdv_soliton(x, c=4.0, x0=-15)
>>> frames, times = kdv_evolve_frames(u0, x, dt=0.001, steps_per_frame=200, n_frames=5)
>>> frames.shape
(6, 256)

physicskit.fields.kdv_soliton(x, c, x0=0.0)[source]#

Exact single-soliton solution of \(u_t + 6uu_x + u_{xxx} = 0\) at \(t=0\).

Parameters:

x (ndarray) – Spatial grid.
c (float) – Soliton speed (equal to twice its amplitude: amplitude \(=c/2\)).
x0 (float) – Initial center position.

Return type:

Returns:

ndarray – \(u(x, 0) = \tfrac{c}{2}\,\mathrm{sech}^2\!\big(\tfrac{\sqrt{c}}{2}(x-x_0)\big)\).

See also

kdv_evolve: Propagate this (or any) initial condition forward in time.

Examples

>>> import numpy as np
>>> round(float(kdv_soliton(0.0, c=4.0)), 6)
2.0

physicskit.fields.kdv_step(u_hat, k, dt)[source]#

Advance the KdV equation in Fourier space by one Strang-split step.

Parameters:

u_hat (ndarray) – Fourier coefficients of \(u\) (as from numpy.fft.fft).
k (ndarray) – Angular wavenumbers matching u_hat.
dt (float) – Time step.

Return type:

Returns:

ndarray – Fourier coefficients after one step of size dt.

physicskit.fields.loop_axial_field(current, radius, z, mu=1.25663706127e-06)[source]#

On-axis field of a circular current loop.

\[B_z(z) = \frac{\mu I R^2}{2 (R^2 + z^2)^{3/2}}\]

(Griffiths Example 5.6; Jackson eq. 5.40 on the axis).

Parameters:

current (float) – Current in amperes.
radius (float) – Loop radius R in meters.
z (array_like) – Axial distance from the loop’s center.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray – \(B_z\) in tesla.

Examples

>>> round(float(loop_axial_field(1.0, 1.0, 0.0) / MU0), 12)  # mu0 I / (2R)
0.5

physicskit.fields.magnetic_dipole_field(m, points, origin=(0.0, 0.0, 0.0), mu=1.25663706127e-06)[source]#

Magnetic field of an ideal point dipole.

\[\mathbf B(\mathbf x) = \frac{\mu}{4\pi} \frac{3\hat{\mathbf n}(\mathbf m\cdot\hat{\mathbf n}) - \mathbf m}{r^3}\]

(Jackson eq. 5.56, away from the origin). A small loop of area \(A\) carrying current \(I\) has \(m = IA\) along its normal, and its Biot-Savart field approaches this one far away.

Parameters:

m (array_like, shape (3,)) – Magnetic moment in A m^2.
points (array_like, shape (..., 3)) – Observation points.
origin (array_like, shape (3,), default=(0, 0, 0)) – Dipole location.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray, shape (…, 3) – Field in tesla.

Examples

>>> B = magnetic_dipole_field([0, 0, 1.0], [[0, 0, 1.0]], mu=4 * np.pi)
>>> B.round(12).tolist()
[[0.0, 0.0, 2.0]]

physicskit.fields.multipole_moments(charges, positions, origin=(0.0, 0.0, 0.0))[source]#

Monopole, dipole and (traceless) quadrupole moments of a charge distribution.

Uses Jackson’s conventions (§4.1, eq. 4.9):

\[q = \sum_i q_i, \quad \mathbf p = \sum_i q_i \mathbf x_i, \quad Q_{jk} = \sum_i q_i\,(3 x_{ij} x_{ik} - r_i^2 \delta_{jk}),\]

with positions measured from origin. A continuous density on a grid is handled the same way: pass charges = rho * dV and the cell centers as positions.

Parameters:

charges (array_like, shape (n,)) – Point charges, in coulombs.
positions (array_like, shape (n, 3) or (n, 2)) – Positions in meters.
origin (array_like, shape (3,), default=(0, 0, 0)) – Expansion center.

Return type:

MultipoleMoments

Returns:

MultipoleMoments

Examples

A physical dipole +-q separated by d along z:

>>> m = multipole_moments([1.0, -1.0], [[0, 0, 0.5], [0, 0, -0.5]])
>>> float(m.monopole), m.dipole.tolist()
(0.0, [0.0, 0.0, 1.0])

physicskit.fields.multipole_potential(moments, points, order=2, eps=8.8541878188e-12)[source]#

Far-field potential from a truncated multipole expansion.

\[\phi(\mathbf x) = \frac{1}{4\pi\varepsilon}\left[\frac{q}{r} + \frac{\mathbf p\cdot\mathbf x}{r^3} + \frac{1}{2}\sum_{jk} Q_{jk}\frac{x_j x_k}{r^5} + \cdots\right]\]

(Jackson eq. 4.10), with \(\mathbf x\) measured from moments.origin. The error of the order-\(\ell\) truncation falls off as \(r^{-(\ell+2)}\).

Parameters:

moments (MultipoleMoments) – From multipole_moments().
points (array_like, shape (..., 3) or (..., 2)) – Observation points.
order (int) – Highest multipole kept: monopole, dipole, or quadrupole.
eps (float) – Permittivity of the medium.

Return type:

Returns:

ndarray, shape (…) – Potential in volts.

Examples

>>> m = multipole_moments([1.0, -1.0], [[0, 0, 0.5], [0, 0, -0.5]])
>>> far = [[0.0, 0.0, 100.0]]
>>> approx = multipole_potential(m, far, order=1, eps=1.0)
>>> exact = coulomb_potential([1.0, -1.0], [[0, 0, 0.5], [0, 0, -0.5]], far, eps=1.0)
>>> bool(abs(approx[0] / exact[0] - 1) < 1e-4)
True

physicskit.fields.nls_bright_soliton(x, t, A=1.0, v=0.0, x0=0.0)[source]#

Exact bright-soliton solution of the focusing NLS equation.

Solves \(i\psi_t + \tfrac{1}{2}\psi_{xx} + |\psi|^2\psi = 0\).

Parameters:

x (ndarray) – Spatial grid.
t (float) – Time.
A (float) – Soliton amplitude.
v (float) – Soliton velocity.
x0 (float) – Initial center position.

Return type:

Returns:

ndarray of complex – \(\psi(x,t) = A\,\mathrm{sech}(A(x-x_0-vt))\, e^{i[v(x-x_0) + (A^2-v^2)t/2]}\).

See also

nls_dark_soliton: The defocusing-equation counterpart.
nls_evolve: Propagate this (or any) initial condition numerically.

Examples

>>> import numpy as np
>>> x = np.array([0.0])
>>> abs(nls_bright_soliton(x, t=0.0, A=1.0))
array([1.])

physicskit.fields.nls_dark_soliton(x, t, rho0=1.0, x0=0.0)[source]#

Exact stationary dark (black) soliton solution of the defocusing NLS equation.

Solves \(i\psi_t + \tfrac{1}{2}\psi_{xx} - |\psi|^2\psi = 0\): a density notch that vanishes at its center, sitting on a uniform background of density rho0.

Parameters:

x (ndarray) – Spatial grid.
t (float) – Time.
rho0 (float) – Background density far from the soliton.
x0 (float) – Center position.

Return type:

Returns:

ndarray of complex – \(\psi(x,t) = \sqrt{\rho_0}\,\tanh\!\big(\sqrt{\rho_0}\,(x-x_0)\big)\, e^{-i\rho_0 t}\), with healing length \(\xi = 1/\sqrt{\rho_0}\).

See also

nls_bright_soliton: The focusing-equation counterpart.

Examples

>>> import numpy as np
>>> x = np.array([0.0])
>>> abs(nls_dark_soliton(x, t=0.0, rho0=1.0))
array([0.])

physicskit.fields.nls_evolve(psi0, x, dt, steps, g=1.0)[source]#

Evolve a nonlinear Schrodinger initial condition via split-step Fourier.

Solves \(i\psi_t + \tfrac{1}{2}\psi_{xx} + g|\psi|^2\psi = 0\) on a periodic domain: g > 0 is focusing (bright solitons), g < 0 is defocusing (dark solitons).

Parameters:

psi0 (ndarray) – Initial wavefunction, sampled on the periodic grid x.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps (int) – Number of steps to advance.
g (float) – Nonlinearity strength and sign.

Return type:

Returns:

ndarray of complex – Field after steps * dt time units.

Examples

A bright soliton propagates without dispersing – its envelope \(|\psi|\) is unchanged after 2 time units of evolution:

>>> import numpy as np
>>> N, L = 1024, 80.0
>>> x = np.linspace(-L / 2, L / 2, N, endpoint=False)
>>> psi0 = nls_bright_soliton(x, t=0.0, A=1.0)
>>> psi = nls_evolve(psi0, x, dt=0.001, steps=2000, g=1.0)
>>> shape_error = np.max(np.abs(np.abs(psi) - np.abs(psi0)))
>>> bool(shape_error < 1e-3)
True

physicskit.fields.nls_evolve_frames(psi0, x, dt, steps_per_frame, n_frames, g=1.0)[source]#

Evolve an NLS field and record a snapshot every steps_per_frame steps, for animation.

A thin looping wrapper around nls_evolve(), analogous to kdv_evolve_frames(): no new physics, just repeated chunked calls with the state carried forward, to produce an animatable sequence of frames.

Parameters:

psi0 (ndarray) – Initial wavefunction.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps_per_frame (int) – Number of integration steps between recorded frames.
n_frames (int) – Number of frames to record after the initial condition.
g (float) – Nonlinearity strength and sign.

Return type:

Returns:

frames (ndarray of complex, shape (n_frames + 1, len(x))) – Field at t=0 and after each recorded chunk.
times (ndarray, shape (n_frames + 1,))

See also

physicskit.fields.visualizers.animate_field_1d: Render these frames as a line-plot animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-40, 40, 512, endpoint=False)
>>> psi0 = nls_bright_soliton(x, t=0.0, A=1.0)
>>> frames, times = nls_evolve_frames(psi0, x, dt=0.001, steps_per_frame=200, n_frames=5, g=1.0)
>>> frames.shape
(6, 512)

physicskit.fields.oscillating_dipole_source(i0, j0, amplitude, freq)[source]#

Build a soft, sinusoidally-oscillating point source for fdtd_2d_tmz_evolve().

Models a simple oscillating dipole antenna: a source current concentrated at a single Yee cell, injected as a soft source (added to the existing field rather than overwriting it, so outgoing waves already at that cell are not clobbered).

Parameters:

i0 (int) – Grid indices of the source cell.
j0 (int) – Grid indices of the source cell.
amplitude (float) – Peak field amplitude injected each step.
freq (float) – Oscillation frequency in Hz.

Returns:

callable – source(Ez, t), suitable for fdtd_2d_tmz_evolve()’s source argument: adds amplitude * sin(2*pi*freq*t) to Ez[i0, j0] in place.

Examples

>>> import numpy as np
>>> Ez = np.zeros((10, 10))
>>> source = oscillating_dipole_source(5, 5, amplitude=2.0, freq=1.0)
>>> source(Ez, t=0.25)
>>> round(float(Ez[5, 5]), 6)
2.0

physicskit.fields.plot_bec_density(X, Y, psi, ax=None)[source]#

Plot the condensate density \(|\psi(\mathbf{r})|^2\) as a heatmap.

Parameters:

X (ndarray) – Real-space coordinate grids.
Y (ndarray) – Real-space coordinate grids.
psi (ndarray) – Condensate wavefunction.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

See also

plot_bec_phase: Plot the phase, which reveals vortex cores as singularities.

Examples

>>> import numpy as np
>>> from physicskit.fields.quantum_fields import harmonic_trap_grid
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(40, 10.0)
>>> psi = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> fig, ax = plot_bec_density(X, Y, psi)
>>> isinstance(fig, plt.Figure)
True

physicskit.fields.plot_bec_phase(X, Y, psi, ax=None)[source]#

Plot the condensate phase \(\arg\psi(\mathbf{r})\); vortex cores appear as \(2\pi\) singularities.

Parameters:

X (ndarray) – Real-space coordinate grids.
Y (ndarray) – Real-space coordinate grids.
psi (ndarray) – Condensate wavefunction.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

See also

plot_bec_density: Plot the density, which shows vortex cores as zeros.
physicskit.fields.quantum_fields.count_vortices: Quantitative vortex detection.

Examples

>>> import numpy as np
>>> from physicskit.fields.quantum_fields import harmonic_trap_grid, gpe_imprint_vortex
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(40, 10.0)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi = gpe_imprint_vortex(psi0, X, Y, [(0.0, 0.0)])
>>> fig, ax = plot_bec_phase(X, Y, psi)
>>> isinstance(fig, plt.Figure)
True

physicskit.fields.plot_field_1d(x, u, ax=None, label=None)[source]#

Plot a 1D field snapshot (a KdV, NLS envelope, or Sine-Gordon profile).

Parameters:

x (ndarray) – Spatial grid.
u (ndarray) – Field values (real; pass numpy.abs(psi) for a complex envelope).
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.
label (str | None) – Legend label for this trace.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.fields.solitons import kdv_soliton
>>> x = np.linspace(-20, 20, 200)
>>> fig, ax = plot_field_1d(x, kdv_soliton(x, c=4.0))
>>> isinstance(fig, plt.Figure)
True

physicskit.fields.plot_poynting_field(X, Y, Sx, Sy, ax=None, stride=4)[source]#

Quiver-plot the Poynting energy-flux field over a 2D grid.

Parameters:

X (ndarray) – Real-space coordinate grids.
Y (ndarray) – Real-space coordinate grids.
Sx (ndarray) – Poynting vector components, as from physicskit.fields.electrodynamics.poynting_vector_tmz().
Sy (ndarray) – Poynting vector components, as from physicskit.fields.electrodynamics.poynting_vector_tmz().
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.
stride (int) – Subsample every stride-th grid point, to keep the quiver plot legible.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> X, Y = np.meshgrid(np.linspace(-1, 1, 20), np.linspace(-1, 1, 20))
>>> Sx, Sy = -Y, X
>>> fig, ax = plot_poynting_field(X, Y, Sx, Sy)
>>> isinstance(fig, plt.Figure)
True

physicskit.fields.pml_conductivity_profile(n_cells, pml_width, dx, order=3, sigma_max=None)[source]#

Build a polynomial-graded electric conductivity profile absorbing at both ends of a 1D grid.

A simplified (non-split-field) approximation to Berenger’s Perfectly Matched Layer: rather than solving auxiliary split-field equations, this grades the ordinary medium’s electric conductivity smoothly up from zero over pml_width cells at each boundary, damping outgoing waves via the standard lossy-dielectric FDTD update (see sigma in fdtd_1d()) instead of reflecting them off a hard wall.

Parameters:

n_cells (int) – Total number of grid points.
pml_width (int) – Number of cells over which the conductivity ramps from 0 to sigma_max.
dx (float) – Spatial grid spacing in meters, used to scale the default sigma_max.
order (int) – Polynomial grading order (Sadiku’s rule of thumb uses a cubic ramp).
sigma_max (float | None) – Peak conductivity at the outermost cell. Defaults to the standard empirical optimum (order + 1) / (150 * pi * dx).

Return type:

Returns:

ndarray, shape (n_cells,) – Conductivity profile, zero in the interior and ramping up at both edges.

Examples

>>> sigma = pml_conductivity_profile(200, pml_width=20, dx=1e-3)
>>> bool(np.all(sigma[20:-20] == 0))
True
>>> bool(sigma[0] > sigma[10] > 0)
True

physicskit.fields.pml_conductivity_profile_2d(shape, pml_width, dx, dy, order=3, sigma_max=None)[source]#

Build a 2D graded-conductivity absorbing boundary by combining two 1D PML profiles.

Reuses pml_conductivity_profile() along each axis and takes the pointwise maximum, so a cell absorbs at whichever rate is larger – correct in the corners, where both the x and y ramps are active.

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
pml_width (int) – Number of cells over which each axis’s conductivity ramps up.
dx (float) – Spatial grid spacing along each axis.
dy (float) – Spatial grid spacing along each axis.
order (int) – Polynomial grading order, passed to pml_conductivity_profile().
sigma_max (float | None) – Peak conductivity; defaults (per axis) to the standard empirical optimum.

Return type:

Returns:

ndarray, shape (Nx, Ny) – Conductivity map, zero in the interior and ramping up at all four edges.

Examples

>>> sigma = pml_conductivity_profile_2d((100, 80), pml_width=10, dx=1e-3, dy=1e-3)
>>> bool(np.all(sigma[10:-10, 10:-10] == 0))
True
>>> bool(sigma[0, 0] > 0)
True

physicskit.fields.poynting_vector_tmz(Ez, Hx, Hy)[source]#

Poynting vector \(\mathbf{S} = \mathbf{E} \times \mathbf{H}\) for a TMz-mode field.

Parameters:

Ez (ndarray) – Out-of-plane electric field.
Hx (ndarray) – In-plane magnetic field components, same shape as Ez.
Hy (ndarray) – In-plane magnetic field components, same shape as Ez.

Return type:

Returns:

Sx, Sy (ndarray) – In-plane energy flux density components, \(S_x = -E_z H_y\), \(S_y = E_z H_x\).

Examples

>>> import numpy as np
>>> Ez = np.array([[2.0]])
>>> Hx = np.array([[3.0]])
>>> Hy = np.array([[5.0]])
>>> poynting_vector_tmz(Ez, Hx, Hy)
(array([[-10.]]), array([[6.]]))

physicskit.fields.sine_gordon_evolve(u0, u0_prev, x, dt, steps)[source]#

Evolve the Sine-Gordon equation with explicit leapfrog finite differences.

Parameters:

u0 (ndarray) – Field at the starting time.
u0_prev (ndarray) – Field one step before the starting time (needed to seed the two-level leapfrog scheme); for a traveling-wave initial condition, use the exact solution evaluated at t=-dt.
x (ndarray) – Uniformly spaced spatial grid, with fixed (Dirichlet) boundaries.
dt (float) – Time step; must satisfy the CFL condition dt <= dx (wave speed 1).
steps (int) – Number of steps to advance.

Return type:

Returns:

u, u_prev (ndarray) – The field at the final step and the step before it (the pair needed to continue the integration further).

Examples

A kink launched at speed \(v=0.5\) arrives at the expected Lorentz-contracted position after propagating:

>>> import numpy as np
>>> N, L, v, x0 = 4000, 200.0, 0.5, -50.0
>>> x = np.linspace(-L / 2, L / 2, N)
>>> dx = x[1] - x[0]
>>> dt = 0.4 * dx
>>> u_prev = sine_gordon_kink(x, -dt, v, x0)
>>> u0 = sine_gordon_kink(x, 0.0, v, x0)
>>> steps = 1000
>>> u, _ = sine_gordon_evolve(u0, u_prev, x, dt, steps)
>>> expected = sine_gordon_kink(x, steps * dt, v, x0)
>>> bool(np.max(np.abs(u - expected)) < 0.01)
True

physicskit.fields.sine_gordon_evolve_frames(u0, u0_prev, x, dt, steps_per_frame, n_frames)[source]#

Evolve a Sine-Gordon field and record a snapshot every steps_per_frame steps, for animation.

A thin looping wrapper around sine_gordon_evolve(), analogous to kdv_evolve_frames(): the (u, u_prev) leapfrog pair is carried forward chunk to chunk, with no change to the underlying physics.

Parameters:

u0 (ndarray) – Field at the starting time.
u0_prev (ndarray) – Field one step before the starting time (see sine_gordon_evolve()).
x (ndarray) – Uniformly spaced spatial grid.
dt (float) – Time step.
steps_per_frame (int) – Number of integration steps between recorded frames.
n_frames (int) – Number of frames to record after the initial condition.

Return type:

Returns:

frames (ndarray, shape (n_frames + 1, len(x))) – Field at t=0 and after each recorded chunk.
times (ndarray, shape (n_frames + 1,))

See also

physicskit.fields.visualizers.animate_field_1d: Render these frames as a line-plot animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-50, 50, 800)
>>> dt = 0.4 * (x[1] - x[0])
>>> u_prev = sine_gordon_kink(x, -dt, v=0.5, x0=-20)
>>> u0 = sine_gordon_kink(x, 0.0, v=0.5, x0=-20)
>>> frames, times = sine_gordon_evolve_frames(u0, u_prev, x, dt, steps_per_frame=100, n_frames=4)
>>> frames.shape
(5, 800)

physicskit.fields.sine_gordon_kink(x, t, v=0.0, x0=0.0, polarity=1)[source]#

Exact kink (or antikink) solution of the Sine-Gordon equation \(u_{tt} - u_{xx} + \sin u = 0\).

Parameters:

x (ndarray) – Spatial grid.
t (float) – Time.
v (float) – Kink velocity (\(|v| < 1\), the wave speed of the linearized equation).
x0 (float) – Initial center position.
polarity (int) – 1 for a kink (field jumps by \(2\pi\)), -1 for an antikink.

Return type:

Returns:

ndarray – \(u(x,t) = 4\arctan\!\big(\exp[\text{polarity}\cdot\gamma(x-x_0-vt)]\big)\), with \(\gamma = 1/\sqrt{1-v^2}\).

See also

sine_gordon_evolve: Propagate this (or any) initial condition numerically.

Examples

>>> import numpy as np
>>> round(float(sine_gordon_kink(0.0, t=0.0)), 6)
3.141593

physicskit.fields.solenoid_axial_field(current, n_turns, radius, length, z, mu=1.25663706127e-06)[source]#

On-axis field of a finite solenoid (uniform current sheet approximation).

\[B_z(z) = \frac{\mu n I}{2}\left[ \frac{z + L/2}{\sqrt{R^2 + (z + L/2)^2}} - \frac{z - L/2}{\sqrt{R^2 + (z - L/2)^2}}\right], \qquad n = N/L,\]

obtained by integrating loop_axial_field() over the winding (Griffiths Problem 5.11). Tends to \(\mu n I\) deep inside a long solenoid and to half that at either end.

Parameters:

current (float) – Current in amperes.
n_turns (int) – Total number of turns N.
radius (float) – Winding radius R in meters.
length (float) – Solenoid length L in meters, centered on z = 0.
z (array_like) – Axial position.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray – \(B_z\) in tesla.

Examples

>>> ratio = solenoid_axial_field(1.0, 1000, 0.01, 1.0, [0.0, 0.5]) / (MU0 * 1000)
>>> ratio.round(4).tolist()  # center, and one end of a long solenoid
[0.9998, 0.5]

physicskit.fields.solenoid_path(radius, length, n_turns, points_per_turn=64, center=(0.0, 0.0, 0.0))[source]#

Vertices of a finite helical solenoid wound along z.

Parameters:

radius (float) – Winding radius in meters.
length (float) – Axial length in meters, from z = -length/2 to +length/2 about center.
n_turns (int) – Number of turns.
points_per_turn (int) – Polyline vertices per turn.
center (array_like, shape (3,), default=(0, 0, 0)) – Solenoid center.

Return type:

Returns:

ndarray, shape (n_turns * points_per_turn + 1, 3) – Helix vertices, winding counter-clockwise seen from +z (field inside along +z for positive current).

Examples

>>> path = solenoid_path(0.01, 0.1, n_turns=10)
>>> round(float(path[-1, 2] - path[0, 2]), 12)
0.1

physicskit.fields.solve_poisson(rho, spacing, *, phi0=None, fixed=None, method='sor', omega=None, tol=1e-08, max_iter=100000, eps=8.8541878188e-12)[source]#

Solve Poisson’s equation on a 2D or 3D grid by Jacobi or SOR relaxation.

Iterates the standard second-order finite-difference form of \(\nabla^2\phi = -\rho/\varepsilon\),

\[\phi_{ij}^{\text{new}} = \frac{(\phi_{i+1,j}+\phi_{i-1,j})/h_x^2 + (\phi_{i,j+1}+\phi_{i,j-1})/h_y^2 + \rho_{ij}/\varepsilon} {2/h_x^2 + 2/h_y^2},\]

either all at once from the previous sweep (Jacobi) or in place with over-relaxation \(\phi \leftarrow \phi + \omega(\phi^{\text{GS}} - \phi)\) (SOR; Press et al., Numerical Recipes, 3rd ed., §20.5). For an N x N grid SOR needs \(O(N)\) sweeps at the optimal \(\omega = 2/(1+\sin(\pi/N))\), against Jacobi’s \(O(N^2)\).

The outermost faces of the grid, and every cell where fixed is true, are Dirichlet cells held at their phi0 values – a grounded box by default, or electrodes/conductors at set voltages. In 2D the geometry is invariant along z, so rho is still a volume density and point sources are really line charges.

Parameters:

rho (array_like, shape (nx, ny) or (nx, ny, nz)) – Charge density in C/m^3.
spacing (float or sequence of float) – Grid spacing in meters, one value or one per axis.
phi0 (array_like, optional) – Initial guess and Dirichlet values (same shape as rho). Defaults to zero everywhere.
fixed (array_like of bool, optional) – Interior cells to hold at phi0 (conductors). The box faces are always held.
method (str) – Relaxation scheme.
omega (float | None) – SOR over-relaxation factor in (0, 2); defaults to 2 / (1 + sin(pi / max(shape))). Ignored for Jacobi.
tol (float) – Stop once the largest single-sweep update is below tol times the largest |phi|.
max_iter (int) – Maximum number of sweeps.
eps (float) – Permittivity of the medium.

Return type:

PoissonSolution

Returns:

PoissonSolution – (phi, n_iter, converged).

Examples

Two plates at +-1 V, grounded box: the potential is odd about the midplane.

>>> phi0 = np.zeros((41, 41))
>>> fixed = np.zeros((41, 41), dtype=bool)
>>> phi0[10:31, 15], phi0[10:31, 25] = 1.0, -1.0
>>> fixed[10:31, 15] = fixed[10:31, 25] = True
>>> sol = solve_poisson(np.zeros((41, 41)), 0.01, phi0=phi0, fixed=fixed)
>>> sol.converged, round(float(sol.phi[20, 20]), 6)
(True, 0.0)

physicskit.fields.solve_poisson_fft(rho, spacing, *, bc='periodic', eps=8.8541878188e-12)[source]#

Solve Poisson’s equation directly with a fast Fourier or sine transform.

bc="periodic" treats the grid as one period of an infinite lattice and solves \(\hat\phi(\mathbf k) = \hat\rho(\mathbf k) / (\varepsilon k^2)\) with the exact continuum \(k^2\) (spectral accuracy for smooth sources). The \(\mathbf k = 0\) mode is set to zero, which is the same as adding a uniform neutralizing background: only a net-neutral periodic system has a periodic potential.

bc="dirichlet" holds every face of the box at \(\phi = 0\) (the grid’s first and last nodes along each axis are the grounded walls) and diagonalizes the second-order finite-difference Laplacian with a type-I discrete sine transform, whose eigenvalues are \(\sum_a (2 - 2\cos(\pi m_a/(n_a-1)))/h_a^2\). It therefore returns the exact solution of the same linear system solve_poisson() iterates towards with a zero phi0 and no fixed cells, in one \(O(N\log N)\) step (Press et al., Numerical Recipes, 3rd ed., §20.4).

Parameters:

rho (array_like) – Charge density in C/m^3, any number of dimensions.
spacing (float or sequence of float) – Grid spacing in meters, one value or one per axis.
bc (str) – Boundary condition, as described above.
eps (float) – Permittivity of the medium.

Return type:

Returns:

ndarray – Potential in volts, same shape as rho.

Examples

A sinusoidal charge wave has potential \(\rho/(\varepsilon k^2)\):

>>> x = np.arange(64) * (2 * np.pi / 64)
>>> phi = solve_poisson_fft(np.sin(x)[:, None] * np.ones((64, 8)), 2 * np.pi / 64, eps=1.0)
>>> bool(np.allclose(phi[:, 0], np.sin(x)))
True

physicskit.fields.tmz_cavity_mode(shape, dx, dy, m, n)[source]#

The analytic \(TM_{mn}\) standing-wave mode of an ideal rectangular PEC cavity.

fdtd_2d_tmz() (and fdtd_2d_tmz_evolve() without sigma) already enforce Ez=0 on all four edges every step – exactly a perfectly-conducting (PEC) cavity wall, no absorbing boundary needed. Seeding the grid with this analytic mode shape (which already vanishes on the boundary) launches a standing wave that should oscillate in place at omega_mn rather than propagate.

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
dx (float) – Grid spacing; the cavity spans Lx=(Nx-1)*dx by Ly=(Ny-1)*dy.
dy (float) – Grid spacing; the cavity spans Lx=(Nx-1)*dx by Ly=(Ny-1)*dy.
m (int) – Mode indices (number of half-wavelengths along each axis).
n (int) – Mode indices (number of half-wavelengths along each axis).

Return type:

Returns:

Ez0 (ndarray, shape (Nx, Ny)) – \(\sin(m\pi x/L_x)\sin(n\pi y/L_y)\).
omega_mn (float) – The mode’s angular frequency, \(c_0\pi\sqrt{(m/L_x)^2+(n/L_y)^2}\).

Examples

>>> Ez0, omega = tmz_cavity_mode((41, 41), dx=1e-3, dy=1e-3, m=1, n=1)
>>> bool(np.all(Ez0[0, :] == 0) and np.all(Ez0[:, 0] == 0))
True
>>> bool(omega > 0)
True

Electrostatics: Coulomb’s law, Poisson’s equation, the method of images, and multipoles.

Everything here solves (or sidesteps solving) Poisson’s equation for the electrostatic potential,

\[\nabla^2 \phi = -\rho / \varepsilon_0, \qquad \mathbf E = -\nabla\phi,\]

in SI units, like the rest of physicskit.fields (EPS0 is the vacuum permittivity):

coulomb_potential() / coulomb_field() – the free-space Green’s function \(1/(4\pi\varepsilon_0 |\mathbf r - \mathbf r'|)\) summed over point charges (Coulomb 1785, plus superposition).
solve_poisson() – Jacobi or successive-over-relaxation (SOR) iteration of the second-order finite-difference Laplacian on a 2D or 3D grid, with Dirichlet values on the box faces and on any interior “conductor” cells (Press et al., Numerical Recipes, 3rd ed., §20.5).
solve_poisson_fft() – a direct spectral solve, periodic (FFT) or grounded-box (type-I discrete sine transform).
image_charges_plane() / image_charges_sphere() – Green’s functions for a grounded plane and a sphere by the method of images (Jackson, Classical Electrodynamics, 3rd ed., §2.1-2.4).
multipole_moments() / multipole_potential() – the Cartesian multipole expansion to quadrupole order (Jackson §4.1, eqs. 4.9-4.10).

class physicskit.fields.electrostatics.MultipoleMoments(monopole, dipole, quadrupole, origin)[source]#

Bases: NamedTuple

Cartesian multipole moments of a charge distribution about origin.

Variables:

monopole (float) – Total charge \(q = \sum_i q_i\).
dipole (ndarray, shape (3,)) – Dipole moment \(\mathbf p = \sum_i q_i \mathbf x_i\).
quadrupole (ndarray, shape (3, 3)) – Traceless quadrupole tensor \(Q_{jk} = \sum_i q_i (3x_{ij}x_{ik} - r_i^2\delta_{jk})\).
origin (ndarray, shape (3,)) – Expansion center the moments are taken about.

Parameters:

monopole (float)
dipole (ndarray)
quadrupole (ndarray)
origin (ndarray)

dipole: ndarray#: Alias for field number 1

monopole: float#: Alias for field number 0

origin: ndarray#: Alias for field number 3

quadrupole: ndarray#: Alias for field number 2

class physicskit.fields.electrostatics.PoissonSolution(phi, n_iter, converged)[source]#

Bases: NamedTuple

Result of solve_poisson().

Variables:

phi (ndarray) – Converged (or last) potential on the grid, in volts.
n_iter (int) – Number of sweeps performed.
converged (bool) – Whether the relative update fell below tol within max_iter.

Parameters:

phi (ndarray)
n_iter (int)
converged (bool)

converged: bool#: Alias for field number 2

n_iter: int#: Alias for field number 1

phi: ndarray#: Alias for field number 0

physicskit.fields.electrostatics.coulomb_field(charges, positions, points, eps=8.8541878188e-12)[source]#

Electric field of a set of point charges (Coulomb’s inverse-square law).

\[\mathbf E(\mathbf r) = \frac{1}{4\pi\varepsilon}\sum_i q_i \frac{\mathbf r - \mathbf r_i}{|\mathbf r - \mathbf r_i|^3}\]

(Jackson eq. 1.5).

Parameters:

charges (array_like, shape (n,)) – Point charges in coulombs.
positions (array_like, shape (n, d)) – Charge positions in meters.
points (array_like, shape (..., d)) – Observation points.
eps (float) – Permittivity of the surrounding medium.

Return type:

Returns:

ndarray, shape (…, d) – Electric field in V/m; nan exactly at a charge.

Examples

>>> E = coulomb_field([1e-9], [[0.0, 0.0, 0.0]], [[2.0, 0.0, 0.0]])
>>> round(float(E[0, 0]), 4)  # a quarter of the field at 1 m
2.2469

physicskit.fields.electrostatics.coulomb_potential(charges, positions, points, eps=8.8541878188e-12)[source]#

Electrostatic potential of a set of point charges (Coulomb’s law plus superposition).

\[\phi(\mathbf r) = \frac{1}{4\pi\varepsilon}\sum_i \frac{q_i}{|\mathbf r - \mathbf r_i|}\]

(Jackson, Classical Electrodynamics, 3rd ed., eq. 1.17).

Parameters:

charges (array_like, shape (n,)) – Point charges in coulombs.
positions (array_like, shape (n, d)) – Charge positions in meters, d = 2 or 3 (2D positions are points in the plane z = 0 of ordinary 3D space).
points (array_like, shape (..., d)) – Observation points.
eps (float) – Permittivity of the surrounding medium.

Return type:

Returns:

ndarray, shape (…) – Potential in volts; inf exactly at a charge.

Examples

>>> phi = coulomb_potential([1e-9], [[0.0, 0.0, 0.0]], [[1.0, 0.0, 0.0]])
>>> round(float(phi[0]), 3)  # 1 nC at 1 m: about 8.988 V
8.988

physicskit.fields.electrostatics.dipole_field(p, points, origin=(0.0, 0.0, 0.0), eps=8.8541878188e-12)[source]#

Electric field of an ideal point dipole.

\[\mathbf E(\mathbf x) = \frac{1}{4\pi\varepsilon} \frac{3\hat{\mathbf n}(\mathbf p\cdot\hat{\mathbf n}) - \mathbf p}{r^3}\]

(Jackson eq. 4.13, away from the origin). On the dipole axis this is \(2p/(4\pi\varepsilon r^3)\), twice the equatorial magnitude.

Parameters:

p (array_like, shape (3,)) – Dipole moment in C m.
points (array_like, shape (..., 3)) – Observation points.
origin (array_like, shape (3,), default=(0, 0, 0)) – Dipole location.
eps (float) – Permittivity of the medium.

Return type:

Returns:

ndarray, shape (…, 3) – Field in V/m.

Examples

>>> E = dipole_field([0, 0, 1.0], [[0, 0, 1.0], [1.0, 0, 0]], eps=1 / (4 * np.pi))
>>> E.round(12).tolist()
[[0.0, 0.0, 2.0], [0.0, 0.0, -1.0]]

physicskit.fields.electrostatics.electric_field_from_potential(phi, spacing)[source]#

Electric field \(\mathbf E = -\nabla\phi\) from a gridded potential.

Uses second-order central differences in the interior (one-sided at the faces), via numpy.gradient().

Parameters:

phi (array_like) – Potential on the grid, in volts.
spacing (float or sequence of float) – Grid spacing in meters, one value or one per axis.

Return type:

Returns:

ndarray, shape (phi.ndim, *phi.shape) – Field components E[0] (along axis 0), E[1], … in V/m.

Examples

>>> x = np.linspace(0.0, 1.0, 11)
>>> E = electric_field_from_potential(np.outer(-5.0 * x, np.ones(4)), 0.1)
>>> bool(np.allclose(E[0], 5.0)), bool(np.allclose(E[1], 0.0))
(True, True)

physicskit.fields.electrostatics.image_charges_plane(q, position)[source]#

Image charge making the grounded conducting plane z = 0 an equipotential.

A charge \(q\) at \((x, y, d)\) above a grounded plane sees, in the region \(z > 0\), exactly the potential of itself plus an image \(-q\) at \((x, y, -d)\) (Jackson §2.1). The pair is Green’s function for the half-space with Dirichlet boundary.

Parameters:

q (float) – The real charge, in coulombs.
position (array_like, shape (3,)) – Its position (x, y, d) with d != 0.

Return type:

Returns:

charges (ndarray, shape (1,)) – The image charge [-q].
positions (ndarray, shape (1, 3)) – The image position [[x, y, -d]].

Examples

>>> qi, ri = image_charges_plane(1e-9, [0.0, 0.0, 0.5])
>>> phi = coulomb_potential(np.r_[1e-9, qi], np.vstack([[0, 0, 0.5], ri]), [[0.3, -0.2, 0.0]])
>>> abs(float(phi[0])) < 1e-12
True

physicskit.fields.electrostatics.image_charges_sphere(q, position, radius, center=(0.0, 0.0, 0.0), sphere_charge=None)[source]#

Image charge(s) for a point charge and a conducting sphere.

For a sphere of radius \(R\) and a charge \(q\) at distance \(a\) from its center, the image

\[q' = -q\frac{R}{a} \quad\text{at distance}\quad \frac{R^2}{a}\]

along the same ray makes the sphere’s surface an equipotential at zero potential (Jackson §2.2, eq. 2.4). This holds with the charge outside (\(a > R\), image inside) or inside a spherical cavity (\(a < R\), image outside). For an isolated sphere carrying total charge \(Q\) (charge outside only), a second image \(Q - q'\) at the center restores the net charge while keeping the surface an equipotential (Jackson §2.3).

Parameters:

q (float) – The real charge, in coulombs.
position (array_like, shape (3,)) – Its position, in meters.
radius (float) – Sphere radius R.
center (array_like, shape (3,), default=(0, 0, 0)) – Sphere center.
sphere_charge (float | None) – Total charge of an isolated sphere. None (default) means the sphere is grounded.

Return type:

Returns:

charges (ndarray, shape (1,) or (2,)) – Image charge(s).
positions (ndarray, shape (1, 3) or (2, 3)) – Image position(s).

Examples

>>> qi, ri = image_charges_sphere(1.0, [2.0, 0.0, 0.0], radius=1.0)
>>> float(qi[0]), float(ri[0, 0])
(-0.5, 0.5)

physicskit.fields.electrostatics.induced_charge_density_plane(q, d, rho)[source]#

Surface charge induced on the grounded plane z = 0 by a charge at height d.

\[\sigma(\varrho) = -\frac{q\, d}{2\pi(\varrho^2 + d^2)^{3/2}}\]

(Griffiths, Introduction to Electrodynamics, 4th ed., eq. 3.10), which integrates to exactly \(-q\) over the plane.

Parameters:

q (float) – The real charge, in coulombs.
d (float) – Its height above the plane, in meters.
rho (array_like) – Distance along the plane from the foot of the perpendicular.

Return type:

Returns:

ndarray – Surface charge density in C/m^2.

Examples

>>> round(float(induced_charge_density_plane(1.0, 1.0, 0.0) * 2 * np.pi), 12)
-1.0

physicskit.fields.electrostatics.induced_charge_density_sphere(q, a, radius, theta)[source]#

Surface charge induced on a grounded sphere by an outside point charge.

\[\sigma(\gamma) = -\frac{q\,(a^2 - R^2)}{4\pi R\,(R^2 + a^2 - 2aR\cos\gamma)^{3/2}}\]

(Jackson eq. 2.5, rewritten in the sphere radius \(R\) and charge distance \(a\)), where \(\gamma\) is the angle from the ray through the charge. It integrates to the image charge \(-qR/a\).

Parameters:

q (float) – The real charge, in coulombs.
a (float) – Distance of the charge from the center, a > radius.
radius (float) – Sphere radius R.
theta (array_like) – Polar angle(s) gamma on the sphere, measured from the charge’s direction.

Return type:

Returns:

ndarray – Surface charge density in C/m^2.

Examples

>>> from scipy.integrate import trapezoid
>>> th = np.linspace(0.0, np.pi, 20001)
>>> sigma = induced_charge_density_sphere(1.0, 2.0, 1.0, th)
>>> total = trapezoid(sigma * 2 * np.pi * np.sin(th), th)
>>> round(float(total), 6)  # the image charge -q R / a
-0.5

physicskit.fields.electrostatics.multipole_moments(charges, positions, origin=(0.0, 0.0, 0.0))[source]#

Monopole, dipole and (traceless) quadrupole moments of a charge distribution.

Uses Jackson’s conventions (§4.1, eq. 4.9):

\[q = \sum_i q_i, \quad \mathbf p = \sum_i q_i \mathbf x_i, \quad Q_{jk} = \sum_i q_i\,(3 x_{ij} x_{ik} - r_i^2 \delta_{jk}),\]

with positions measured from origin. A continuous density on a grid is handled the same way: pass charges = rho * dV and the cell centers as positions.

Parameters:

charges (array_like, shape (n,)) – Point charges, in coulombs.
positions (array_like, shape (n, 3) or (n, 2)) – Positions in meters.
origin (array_like, shape (3,), default=(0, 0, 0)) – Expansion center.

Return type:

MultipoleMoments

Returns:

MultipoleMoments

Examples

A physical dipole +-q separated by d along z:

>>> m = multipole_moments([1.0, -1.0], [[0, 0, 0.5], [0, 0, -0.5]])
>>> float(m.monopole), m.dipole.tolist()
(0.0, [0.0, 0.0, 1.0])

physicskit.fields.electrostatics.multipole_potential(moments, points, order=2, eps=8.8541878188e-12)[source]#

Far-field potential from a truncated multipole expansion.

\[\phi(\mathbf x) = \frac{1}{4\pi\varepsilon}\left[\frac{q}{r} + \frac{\mathbf p\cdot\mathbf x}{r^3} + \frac{1}{2}\sum_{jk} Q_{jk}\frac{x_j x_k}{r^5} + \cdots\right]\]

(Jackson eq. 4.10), with \(\mathbf x\) measured from moments.origin. The error of the order-\(\ell\) truncation falls off as \(r^{-(\ell+2)}\).

Parameters:

moments (MultipoleMoments) – From multipole_moments().
points (array_like, shape (..., 3) or (..., 2)) – Observation points.
order (int) – Highest multipole kept: monopole, dipole, or quadrupole.
eps (float) – Permittivity of the medium.

Return type:

Returns:

ndarray, shape (…) – Potential in volts.

Examples

>>> m = multipole_moments([1.0, -1.0], [[0, 0, 0.5], [0, 0, -0.5]])
>>> far = [[0.0, 0.0, 100.0]]
>>> approx = multipole_potential(m, far, order=1, eps=1.0)
>>> exact = coulomb_potential([1.0, -1.0], [[0, 0, 0.5], [0, 0, -0.5]], far, eps=1.0)
>>> bool(abs(approx[0] / exact[0] - 1) < 1e-4)
True

physicskit.fields.electrostatics.solve_poisson(rho, spacing, *, phi0=None, fixed=None, method='sor', omega=None, tol=1e-08, max_iter=100000, eps=8.8541878188e-12)[source]#

Solve Poisson’s equation on a 2D or 3D grid by Jacobi or SOR relaxation.

Iterates the standard second-order finite-difference form of \(\nabla^2\phi = -\rho/\varepsilon\),

\[\phi_{ij}^{\text{new}} = \frac{(\phi_{i+1,j}+\phi_{i-1,j})/h_x^2 + (\phi_{i,j+1}+\phi_{i,j-1})/h_y^2 + \rho_{ij}/\varepsilon} {2/h_x^2 + 2/h_y^2},\]

either all at once from the previous sweep (Jacobi) or in place with over-relaxation \(\phi \leftarrow \phi + \omega(\phi^{\text{GS}} - \phi)\) (SOR; Press et al., Numerical Recipes, 3rd ed., §20.5). For an N x N grid SOR needs \(O(N)\) sweeps at the optimal \(\omega = 2/(1+\sin(\pi/N))\), against Jacobi’s \(O(N^2)\).

The outermost faces of the grid, and every cell where fixed is true, are Dirichlet cells held at their phi0 values – a grounded box by default, or electrodes/conductors at set voltages. In 2D the geometry is invariant along z, so rho is still a volume density and point sources are really line charges.

Parameters:

rho (array_like, shape (nx, ny) or (nx, ny, nz)) – Charge density in C/m^3.
spacing (float or sequence of float) – Grid spacing in meters, one value or one per axis.
phi0 (array_like, optional) – Initial guess and Dirichlet values (same shape as rho). Defaults to zero everywhere.
fixed (array_like of bool, optional) – Interior cells to hold at phi0 (conductors). The box faces are always held.
method (str) – Relaxation scheme.
omega (float | None) – SOR over-relaxation factor in (0, 2); defaults to 2 / (1 + sin(pi / max(shape))). Ignored for Jacobi.
tol (float) – Stop once the largest single-sweep update is below tol times the largest |phi|.
max_iter (int) – Maximum number of sweeps.
eps (float) – Permittivity of the medium.

Return type:

PoissonSolution

Returns:

PoissonSolution – (phi, n_iter, converged).

Examples

Two plates at +-1 V, grounded box: the potential is odd about the midplane.

>>> phi0 = np.zeros((41, 41))
>>> fixed = np.zeros((41, 41), dtype=bool)
>>> phi0[10:31, 15], phi0[10:31, 25] = 1.0, -1.0
>>> fixed[10:31, 15] = fixed[10:31, 25] = True
>>> sol = solve_poisson(np.zeros((41, 41)), 0.01, phi0=phi0, fixed=fixed)
>>> sol.converged, round(float(sol.phi[20, 20]), 6)
(True, 0.0)

physicskit.fields.electrostatics.solve_poisson_fft(rho, spacing, *, bc='periodic', eps=8.8541878188e-12)[source]#

Solve Poisson’s equation directly with a fast Fourier or sine transform.

bc="periodic" treats the grid as one period of an infinite lattice and solves \(\hat\phi(\mathbf k) = \hat\rho(\mathbf k) / (\varepsilon k^2)\) with the exact continuum \(k^2\) (spectral accuracy for smooth sources). The \(\mathbf k = 0\) mode is set to zero, which is the same as adding a uniform neutralizing background: only a net-neutral periodic system has a periodic potential.

bc="dirichlet" holds every face of the box at \(\phi = 0\) (the grid’s first and last nodes along each axis are the grounded walls) and diagonalizes the second-order finite-difference Laplacian with a type-I discrete sine transform, whose eigenvalues are \(\sum_a (2 - 2\cos(\pi m_a/(n_a-1)))/h_a^2\). It therefore returns the exact solution of the same linear system solve_poisson() iterates towards with a zero phi0 and no fixed cells, in one \(O(N\log N)\) step (Press et al., Numerical Recipes, 3rd ed., §20.4).

Parameters:

rho (array_like) – Charge density in C/m^3, any number of dimensions.
spacing (float or sequence of float) – Grid spacing in meters, one value or one per axis.
bc (str) – Boundary condition, as described above.
eps (float) – Permittivity of the medium.

Return type:

Returns:

ndarray – Potential in volts, same shape as rho.

Examples

A sinusoidal charge wave has potential \(\rho/(\varepsilon k^2)\):

>>> x = np.arange(64) * (2 * np.pi / 64)
>>> phi = solve_poisson_fft(np.sin(x)[:, None] * np.ones((64, 8)), 2 * np.pi / 64, eps=1.0)
>>> bool(np.allclose(phi[:, 0], np.sin(x)))
True

Magnetostatics: the Biot-Savart law for current-carrying wires, loops, and solenoids.

A steady current \(I\) along a wire produces the magnetic field

\[\mathbf B(\mathbf r) = \frac{\mu_0 I}{4\pi}\oint \frac{d\boldsymbol\ell \times (\mathbf r - \mathbf r')}{|\mathbf r - \mathbf r'|^3}\]

(Biot and Savart 1820; Jackson, Classical Electrodynamics, 3rd ed., eq. 5.4). biot_savart_field() evaluates it for any wire given as a polyline, integrating each straight segment in closed form rather than by quadrature, so the only discretization error is the polygonal approximation to a curved wire. Units are SI, like the rest of physicskit.fields (MU0 is the vacuum permeability).

physicskit.fields.magnetostatics.biot_savart_field(path, points, current=1.0, mu=1.25663706127e-06)[source]#

Magnetic field of a current-carrying polyline wire, by the Biot-Savart law.

Each straight segment from \(\mathbf a\) to \(\mathbf b\) is integrated exactly: with \(\mathbf r_1 = \mathbf r - \mathbf a\) and \(\mathbf r_2 = \mathbf r - \mathbf b\),

\[\mathbf B_{\text{seg}} = \frac{\mu I}{4\pi} \frac{(r_1 + r_2)\,(\mathbf r_1 \times \mathbf r_2)} {r_1 r_2\,(r_1 r_2 + \mathbf r_1\cdot\mathbf r_2)},\]

the vector form of the textbook finite-wire result \(B = \mu I(\sin\theta_2 - \sin\theta_1)/(4\pi s)\) (Griffiths, Introduction to Electrodynamics, 4th ed., Example 5.5; J. D. Hanson and S. P. Hirshman, Phys. Plasmas 9, 4410 (2002)). Points on a segment or its straight-line extension get no contribution from it (the field of a segment vanishes on its own axis).

Parameters:

path (array_like, shape (n_vertices, 3)) – Wire vertices, in meters, in the direction of current flow. Close a loop by repeating the first vertex at the end (as circular_loop_path() does).
points (array_like, shape (..., 3)) – Observation points.
current (float) – Current in amperes.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray, shape (…, 3) – Magnetic field in tesla.

Examples

A long straight wire recovers \(\mu_0 I / (2\pi s)\):

>>> wire = [[0, 0, -1e4], [0, 0, 1e4]]
>>> B = biot_savart_field(wire, [[0.1, 0.0, 0.0]], current=1.0)
>>> round(float(B[0, 1] / infinite_wire_field(1.0, 0.1)), 6)
1.0

physicskit.fields.magnetostatics.circular_loop_path(radius, n_segments=256, center=(0.0, 0.0, 0.0))[source]#

Vertices of a closed circular loop in a plane of constant z.

The loop runs counter-clockwise seen from +z, so a positive current produces a field along +z at its center.

Parameters:

radius (float) – Loop radius in meters.
n_segments (int) – Number of straight segments; the polygon’s on-axis field differs from a true circle’s by \(O(1/n^2)\).
center (array_like, shape (3,), default=(0, 0, 0)) – Loop center.

Return type:

Returns:

ndarray, shape (n_segments + 1, 3) – Vertices, with the first repeated at the end.

Examples

>>> path = circular_loop_path(1.0, n_segments=4)
>>> path.shape, bool(np.allclose(path[0], path[-1]))
((5, 3), True)

physicskit.fields.magnetostatics.infinite_wire_field(current, s, mu=1.25663706127e-06)[source]#

Field magnitude at distance s from an infinite straight wire, \(\mu I / (2\pi s)\).

(Griffiths eq. 5.38.) The field circles the wire in the right-hand sense about the current.

Parameters:

current (float) – Current in amperes.
s (array_like) – Perpendicular distance from the wire, in meters.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray – \(|\mathbf B|\) in tesla.

Examples

>>> round(float(infinite_wire_field(1.0, 1.0)), 10)  # 2e-7 T per ampere at 1 m
2e-07

physicskit.fields.magnetostatics.loop_axial_field(current, radius, z, mu=1.25663706127e-06)[source]#

On-axis field of a circular current loop.

\[B_z(z) = \frac{\mu I R^2}{2 (R^2 + z^2)^{3/2}}\]

(Griffiths Example 5.6; Jackson eq. 5.40 on the axis).

Parameters:

current (float) – Current in amperes.
radius (float) – Loop radius R in meters.
z (array_like) – Axial distance from the loop’s center.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray – \(B_z\) in tesla.

Examples

>>> round(float(loop_axial_field(1.0, 1.0, 0.0) / MU0), 12)  # mu0 I / (2R)
0.5

physicskit.fields.magnetostatics.magnetic_dipole_field(m, points, origin=(0.0, 0.0, 0.0), mu=1.25663706127e-06)[source]#

Magnetic field of an ideal point dipole.

\[\mathbf B(\mathbf x) = \frac{\mu}{4\pi} \frac{3\hat{\mathbf n}(\mathbf m\cdot\hat{\mathbf n}) - \mathbf m}{r^3}\]

(Jackson eq. 5.56, away from the origin). A small loop of area \(A\) carrying current \(I\) has \(m = IA\) along its normal, and its Biot-Savart field approaches this one far away.

Parameters:

m (array_like, shape (3,)) – Magnetic moment in A m^2.
points (array_like, shape (..., 3)) – Observation points.
origin (array_like, shape (3,), default=(0, 0, 0)) – Dipole location.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray, shape (…, 3) – Field in tesla.

Examples

>>> B = magnetic_dipole_field([0, 0, 1.0], [[0, 0, 1.0]], mu=4 * np.pi)
>>> B.round(12).tolist()
[[0.0, 0.0, 2.0]]

physicskit.fields.magnetostatics.solenoid_axial_field(current, n_turns, radius, length, z, mu=1.25663706127e-06)[source]#

On-axis field of a finite solenoid (uniform current sheet approximation).

\[B_z(z) = \frac{\mu n I}{2}\left[ \frac{z + L/2}{\sqrt{R^2 + (z + L/2)^2}} - \frac{z - L/2}{\sqrt{R^2 + (z - L/2)^2}}\right], \qquad n = N/L,\]

obtained by integrating loop_axial_field() over the winding (Griffiths Problem 5.11). Tends to \(\mu n I\) deep inside a long solenoid and to half that at either end.

Parameters:

current (float) – Current in amperes.
n_turns (int) – Total number of turns N.
radius (float) – Winding radius R in meters.
length (float) – Solenoid length L in meters, centered on z = 0.
z (array_like) – Axial position.
mu (float) – Permeability of the medium.

Return type:

Returns:

ndarray – \(B_z\) in tesla.

Examples

>>> ratio = solenoid_axial_field(1.0, 1000, 0.01, 1.0, [0.0, 0.5]) / (MU0 * 1000)
>>> ratio.round(4).tolist()  # center, and one end of a long solenoid
[0.9998, 0.5]

physicskit.fields.magnetostatics.solenoid_path(radius, length, n_turns, points_per_turn=64, center=(0.0, 0.0, 0.0))[source]#

Vertices of a finite helical solenoid wound along z.

Parameters:

radius (float) – Winding radius in meters.
length (float) – Axial length in meters, from z = -length/2 to +length/2 about center.
n_turns (int) – Number of turns.
points_per_turn (int) – Polyline vertices per turn.
center (array_like, shape (3,), default=(0, 0, 0)) – Solenoid center.

Return type:

Returns:

ndarray, shape (n_turns * points_per_turn + 1, 3) – Helix vertices, winding counter-clockwise seen from +z (field inside along +z for positive current).

Examples

>>> path = solenoid_path(0.01, 0.1, n_turns=10)
>>> round(float(path[-1, 2] - path[0, 2]), 12)
0.1

Finite-difference time-domain (FDTD) solutions of Maxwell’s equations on a Yee grid.

Implements the staggered-grid (Yee 1966) leapfrog update for 1D plane waves and 2D TMz-mode propagation, plus a graded-conductivity absorbing boundary (a simplified, non-split-field approximation to Berenger’s Perfectly Matched Layer) that damps outgoing waves instead of reflecting them off a hard electric wall.

Units are SI throughout; EPS0, MU0, C0 are the vacuum permittivity, permeability, and speed of light.

physicskit.fields.electrodynamics.courant_limit_1d(dx)[source]#

Maximum stable time step for the 1D FDTD update.

Parameters:: dx (float) – Spatial grid spacing in meters.
Return type:: float
Returns:: float – The Courant stability limit \(\Delta t = \Delta x / c_0\).

Examples

>>> round(float(courant_limit_1d(1e-3) * C0 / 1e-3), 6)
1.0

physicskit.fields.electrodynamics.courant_limit_2d(dx, dy)[source]#

Maximum stable time step for the 2D FDTD update.

Parameters:

dx (float) – Spatial grid spacing in meters along each axis.
dy (float) – Spatial grid spacing in meters along each axis.

Return type:

Returns:

float – The Courant stability limit \(\Delta t = 1 / (c_0 \sqrt{1/\Delta x^2 + 1/\Delta y^2})\).

Examples

>>> dt = courant_limit_2d(1e-3, 1e-3)
>>> round(float(dt * C0 / 1e-3), 6)
0.707107

physicskit.fields.electrodynamics.dielectric_slab(shape, i_start, i_end, eps_r_slab=4.0)[source]#

A planar dielectric slab: eps_r_slab between grid columns i_start and i_end, vacuum elsewhere.

A minimal demo permittivity map for fdtd_2d_tmz() / fdtd_2d_tmz_evolve(): a wave launched from one side crossing i_start partially reflects and partially refracts/transmits into the slab, and again at i_end on the way out, illustrating propagation through a medium interface.

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
i_start (int) – Grid-index bounds of the slab along the x axis (rows i_start:i_end).
i_end (int) – Grid-index bounds of the slab along the x axis (rows i_start:i_end).
eps_r_slab (float) – Relative permittivity inside the slab.

Return type:

Returns:

ndarray, shape (Nx, Ny) – eps_r map: eps_r_slab inside the slab, 1.0 outside.

Examples

>>> eps_r = dielectric_slab((40, 20), i_start=15, i_end=25, eps_r_slab=4.0)
>>> float(eps_r[10, 5]), float(eps_r[20, 5])
(1.0, 4.0)

physicskit.fields.electrodynamics.fdtd_1d(Ez0, Hy0, eps_r, mu_r, steps, dt, dx, sigma=None)[source]#

Evolve a 1D plane wave (Ez, Hy) on a Yee grid using leapfrog FDTD.

Discretizes Maxwell’s curl equations for fields varying only along x, \(\partial_t H_y = \partial_x E_z/\mu\) and \(\partial_t E_z = \partial_x H_y/\epsilon\) (the same signs as fdtd_2d_tmz()), so a wave moving toward +x has \(H_y = -E_z/\eta\) (Poynting vector \(-E_zH_y>0\)).

Parameters:

Ez0 (ndarray) – Initial electric field.
Hy0 (ndarray) – Initial magnetic field, staggered half a cell ahead of Ez0.
eps_r (ndarray) – Relative permittivity and permeability at each Ez grid point.
mu_r (ndarray) – Relative permittivity and permeability at each Ez grid point.
steps (int) – Number of time steps to advance.
dt (float) – Time step in seconds; must satisfy courant_limit_1d().
dx (float) – Spatial grid spacing in meters.
sigma (ndarray | None) – Electric conductivity (S/m) at each grid point, e.g. from pml_conductivity_profile(), for a lossy/absorbing boundary. Defaults to a lossless grid with hard (PEC) walls at both ends.

Return type:

Returns:

Ez, Hy (ndarray) – Fields at the final time step, same shapes as Ez0, Hy0.

See also

fdtd_2d_tmz: The 2D (TMz-mode) analog of this solver.

Examples

A smooth Gaussian pulse, launched with the impedance-matched initial condition of a purely right-moving wave (\(H_y=-E_z/\eta_0\)), propagates at exactly the vacuum speed of light:

>>> import numpy as np
>>> N = 800
>>> dx = 1e-3
>>> dt = 0.99 * courant_limit_1d(dx)
>>> eta0 = np.sqrt(MU0 / EPS0)
>>> x0, sigma_pulse = 100, 25
>>> Ez0 = np.exp(-((np.arange(N) - x0) ** 2) / (2 * sigma_pulse ** 2))
>>> xh = np.arange(N - 1) + 0.5
>>> Hy0 = -np.exp(-((xh - x0) ** 2) / (2 * sigma_pulse ** 2)) / eta0
>>> eps_r, mu_r = np.ones(N), np.ones(N)
>>> Ez, Hy = fdtd_1d(Ez0, Hy0, eps_r, mu_r, steps=360, dt=dt, dx=dx)
>>> peak = np.argmax(Ez)
>>> measured_speed = (peak - x0) * dx / (360 * dt)
>>> round(float(measured_speed / C0), 3)
1.002

physicskit.fields.electrodynamics.fdtd_2d_tmz(Ez0, Hx0, Hy0, eps_r, mu_r, steps, dt, dx, dy)[source]#

Evolve a 2D TMz-mode field (Ez, Hx, Hy) on a Yee grid using leapfrog FDTD.

Parameters:

Ez0 (ndarray) – Initial (out-of-plane) electric field.
Hx0 (ndarray) – Initial in-plane magnetic field components.
Hy0 (ndarray) – Initial in-plane magnetic field components.
eps_r (ndarray) – Relative permittivity and permeability maps.
mu_r (ndarray) – Relative permittivity and permeability maps.
steps (int) – Number of time steps to advance.
dt (float) – Time step in seconds; must satisfy courant_limit_2d().
dx (float) – Spatial grid spacing in meters along each axis.
dy (float) – Spatial grid spacing in meters along each axis.

Return type:

Returns:

Ez, Hx, Hy (ndarray, shape (Nx, Ny)) – Fields at the final time step.

See also

fdtd_1d: The 1D analog of this solver.
poynting_vector_tmz: Energy flux from the evolved fields.

Examples

>>> import numpy as np
>>> Nx, Ny = 50, 50
>>> Ez0 = np.zeros((Nx, Ny))
>>> Hx0 = np.zeros((Nx, Ny))
>>> Hy0 = np.zeros((Nx, Ny))
>>> eps_r, mu_r = np.ones((Nx, Ny)), np.ones((Nx, Ny))
>>> dx = dy = 1e-3
>>> dt = 0.5 * courant_limit_2d(dx, dy)
>>> Ez0[Nx // 2, Ny // 2] = 1.0
>>> Ez, Hx, Hy = fdtd_2d_tmz(Ez0, Hx0, Hy0, eps_r, mu_r, steps=10, dt=dt, dx=dx, dy=dy)
>>> Ez.shape
(50, 50)

physicskit.fields.electrodynamics.fdtd_2d_tmz_evolve(Ez0, Hx0, Hy0, eps_r, mu_r, steps, dt, dx, dy, sigma=None, source=None, snapshot_stride=1)[source]#

Time-step the 2D TMz FDTD update while recording Ez snapshots, for animation.

Unlike fdtd_2d_tmz(), which returns only the field at the final step, this driver records Ez every snapshot_stride steps and supports two features layered on top of the same base update: a lossy absorbing boundary (sigma, e.g. from pml_conductivity_profile_2d(), so an outgoing wave is damped instead of reflecting off the hard PEC walls) and a soft source callback injected into Ez after every field update (e.g. oscillating_dipole_source()).

Parameters:

Ez0 (ndarray) – As in fdtd_2d_tmz().
Hx0 (ndarray) – As in fdtd_2d_tmz().
Hy0 (ndarray) – As in fdtd_2d_tmz().
eps_r (ndarray) – As in fdtd_2d_tmz().
mu_r (ndarray) – As in fdtd_2d_tmz().
steps (int) – As in fdtd_2d_tmz().
dt (float) – As in fdtd_2d_tmz().
dx (float) – As in fdtd_2d_tmz().
dy (float) – As in fdtd_2d_tmz().
sigma (ndarray | None) – Electric conductivity map. None (default) reproduces fdtd_2d_tmz()’s lossless, hard-wall behavior exactly.
source (callable, optional) – source(Ez, t), called after each field update and expected to mutate Ez in place (a soft source: added to, not overwriting, the field already there).
snapshot_stride (int) – Record Ez every this many steps (plus the initial condition).

Return type:

Returns:

snapshots (ndarray, shape (n_recorded, Nx, Ny)) – Ez at t=0 and after every recorded step.
times (ndarray, shape (n_recorded,)) – Time of each recorded snapshot.

See also

fdtd_2d_tmz: The single-shot (final-state-only) solver this wraps.

Examples

>>> import numpy as np
>>> Nx, Ny = 40, 40
>>> Ez0 = np.zeros((Nx, Ny))
>>> Hx0 = Hy0 = np.zeros((Nx, Ny))
>>> eps_r = mu_r = np.ones((Nx, Ny))
>>> dx = dy = 1e-3
>>> dt = 0.5 * courant_limit_2d(dx, dy)
>>> source = oscillating_dipole_source(Nx // 2, Ny // 2, amplitude=1.0, freq=5e10)
>>> snaps, times = fdtd_2d_tmz_evolve(Ez0, Hx0, Hy0, eps_r, mu_r, steps=20, dt=dt, dx=dx, dy=dy, source=source, snapshot_stride=5)
>>> snaps.shape
(5, 40, 40)

physicskit.fields.electrodynamics.flux_tube_energy_density_2d(shape, dx, dy, separation, flux_quantum=1.0, tube_width=None)[source]#

2D energy-density map of the confined flux tube, for visualizing it as a stretching heatmap.

Extrudes flux_tube_field_1d() (evaluated along the axis joining the two charges) into 2D by giving it a fixed-width Gaussian transverse profile – the “fixed cross-section” that produces confinement – and forms the field energy density \(\tfrac{1}{2}\,\text{field}(x)^2\times\text{transverse}(y)^2\).

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
dx (float) – Grid spacing; the charges are placed on the centered grid at \(x=\pm\text{separation}/2\), \(y=0\).
dy (float) – Grid spacing; the charges are placed on the centered grid at \(x=\pm\text{separation}/2\), \(y=0\).
separation (float) – Distance between the two charges.
flux_quantum (float) – Charge magnitude, as in flux_tube_field_1d().
tube_width (float | None) – Transverse (y) width of the flux tube. Defaults to 5 * dy.

Return type:

Returns:

ndarray, shape (Nx, Ny) – Energy density, largest along the tube connecting the charges and decaying away from its axis.

Examples

>>> e = flux_tube_energy_density_2d((200, 40), dx=0.2, dy=0.2, separation=10.0)
>>> bool(e.sum() > 0)
True

physicskit.fields.electrodynamics.flux_tube_field_1d(x, separation, flux_quantum=1.0, wall_width=None)[source]#

Static 1D confined flux profile between two opposite point charges (toy confinement model).

Solves the 1D Gauss law \(d(\text{flux})/dx = \rho(x)\) for a field confined to a fixed-cross-section tube connecting two opposite charges at \(x=\pm\text{separation}/2\), each smoothed into a Gaussian of width wall_width (a numerical regularization of the point charge, not a physical effect). The field is computed by direct numerical integration (cumulative trapezoidal sum) of this charge density, so it is genuinely “solved for” at each separation rather than assumed.

Parameters:

x (ndarray) – 1D grid the charges live on.
separation (float) – Distance between the two charges.
flux_quantum (float) – Charge magnitude (equivalently, the plateau field strength inside the tube).
wall_width (float | None) – Smoothing width of each point charge. Defaults to 3 * (x[1] - x[0]), i.e. a few grid cells – small enough to look like a point charge but resolved on the grid.

Return type:

Returns:

ndarray, same shape as x – The confined field, approximately +flux_quantum between the charges and 0 outside, with smooth transitions of width wall_width.

See also

flux_tube_energy_density_2d: The corresponding 2D energy-density map, for animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-20, 20, 2000)
>>> field = flux_tube_field_1d(x, separation=10.0, flux_quantum=1.0)
>>> bool(abs(field[np.argmin(np.abs(x))] - 1.0) < 0.05)
True
>>> bool(abs(field[0]) < 0.05)
True

physicskit.fields.electrodynamics.oscillating_dipole_source(i0, j0, amplitude, freq)[source]#

Build a soft, sinusoidally-oscillating point source for fdtd_2d_tmz_evolve().

Models a simple oscillating dipole antenna: a source current concentrated at a single Yee cell, injected as a soft source (added to the existing field rather than overwriting it, so outgoing waves already at that cell are not clobbered).

Parameters:

i0 (int) – Grid indices of the source cell.
j0 (int) – Grid indices of the source cell.
amplitude (float) – Peak field amplitude injected each step.
freq (float) – Oscillation frequency in Hz.

Returns:

callable – source(Ez, t), suitable for fdtd_2d_tmz_evolve()’s source argument: adds amplitude * sin(2*pi*freq*t) to Ez[i0, j0] in place.

Examples

>>> import numpy as np
>>> Ez = np.zeros((10, 10))
>>> source = oscillating_dipole_source(5, 5, amplitude=2.0, freq=1.0)
>>> source(Ez, t=0.25)
>>> round(float(Ez[5, 5]), 6)
2.0

physicskit.fields.electrodynamics.pml_conductivity_profile(n_cells, pml_width, dx, order=3, sigma_max=None)[source]#

Build a polynomial-graded electric conductivity profile absorbing at both ends of a 1D grid.

A simplified (non-split-field) approximation to Berenger’s Perfectly Matched Layer: rather than solving auxiliary split-field equations, this grades the ordinary medium’s electric conductivity smoothly up from zero over pml_width cells at each boundary, damping outgoing waves via the standard lossy-dielectric FDTD update (see sigma in fdtd_1d()) instead of reflecting them off a hard wall.

Parameters:

n_cells (int) – Total number of grid points.
pml_width (int) – Number of cells over which the conductivity ramps from 0 to sigma_max.
dx (float) – Spatial grid spacing in meters, used to scale the default sigma_max.
order (int) – Polynomial grading order (Sadiku’s rule of thumb uses a cubic ramp).
sigma_max (float | None) – Peak conductivity at the outermost cell. Defaults to the standard empirical optimum (order + 1) / (150 * pi * dx).

Return type:

Returns:

ndarray, shape (n_cells,) – Conductivity profile, zero in the interior and ramping up at both edges.

Examples

>>> sigma = pml_conductivity_profile(200, pml_width=20, dx=1e-3)
>>> bool(np.all(sigma[20:-20] == 0))
True
>>> bool(sigma[0] > sigma[10] > 0)
True

physicskit.fields.electrodynamics.pml_conductivity_profile_2d(shape, pml_width, dx, dy, order=3, sigma_max=None)[source]#

Build a 2D graded-conductivity absorbing boundary by combining two 1D PML profiles.

Reuses pml_conductivity_profile() along each axis and takes the pointwise maximum, so a cell absorbs at whichever rate is larger – correct in the corners, where both the x and y ramps are active.

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
pml_width (int) – Number of cells over which each axis’s conductivity ramps up.
dx (float) – Spatial grid spacing along each axis.
dy (float) – Spatial grid spacing along each axis.
order (int) – Polynomial grading order, passed to pml_conductivity_profile().
sigma_max (float | None) – Peak conductivity; defaults (per axis) to the standard empirical optimum.

Return type:

Returns:

ndarray, shape (Nx, Ny) – Conductivity map, zero in the interior and ramping up at all four edges.

Examples

>>> sigma = pml_conductivity_profile_2d((100, 80), pml_width=10, dx=1e-3, dy=1e-3)
>>> bool(np.all(sigma[10:-10, 10:-10] == 0))
True
>>> bool(sigma[0, 0] > 0)
True

physicskit.fields.electrodynamics.poynting_vector_tmz(Ez, Hx, Hy)[source]#

Poynting vector \(\mathbf{S} = \mathbf{E} \times \mathbf{H}\) for a TMz-mode field.

Parameters:

Ez (ndarray) – Out-of-plane electric field.
Hx (ndarray) – In-plane magnetic field components, same shape as Ez.
Hy (ndarray) – In-plane magnetic field components, same shape as Ez.

Return type:

Returns:

Sx, Sy (ndarray) – In-plane energy flux density components, \(S_x = -E_z H_y\), \(S_y = E_z H_x\).

Examples

>>> import numpy as np
>>> Ez = np.array([[2.0]])
>>> Hx = np.array([[3.0]])
>>> Hy = np.array([[5.0]])
>>> poynting_vector_tmz(Ez, Hx, Hy)
(array([[-10.]]), array([[6.]]))

physicskit.fields.electrodynamics.tmz_cavity_mode(shape, dx, dy, m, n)[source]#

The analytic \(TM_{mn}\) standing-wave mode of an ideal rectangular PEC cavity.

fdtd_2d_tmz() (and fdtd_2d_tmz_evolve() without sigma) already enforce Ez=0 on all four edges every step – exactly a perfectly-conducting (PEC) cavity wall, no absorbing boundary needed. Seeding the grid with this analytic mode shape (which already vanishes on the boundary) launches a standing wave that should oscillate in place at omega_mn rather than propagate.

Parameters:

shape (tuple) – (Nx, Ny) grid shape.
dx (float) – Grid spacing; the cavity spans Lx=(Nx-1)*dx by Ly=(Ny-1)*dy.
dy (float) – Grid spacing; the cavity spans Lx=(Nx-1)*dx by Ly=(Ny-1)*dy.
m (int) – Mode indices (number of half-wavelengths along each axis).
n (int) – Mode indices (number of half-wavelengths along each axis).

Return type:

Returns:

Ez0 (ndarray, shape (Nx, Ny)) – \(\sin(m\pi x/L_x)\sin(n\pi y/L_y)\).
omega_mn (float) – The mode’s angular frequency, \(c_0\pi\sqrt{(m/L_x)^2+(n/L_y)^2}\).

Examples

>>> Ez0, omega = tmz_cavity_mode((41, 41), dx=1e-3, dy=1e-3, m=1, n=1)
>>> bool(np.all(Ez0[0, :] == 0) and np.all(Ez0[:, 0] == 0))
True
>>> bool(omega > 0)
True

Soliton-bearing nonlinear field equations: KdV, the nonlinear Schrodinger equation, and Sine-Gordon.

Each equation is integrated with a scheme suited to its structure:

Korteweg-de Vries – Strang-split pseudo-spectral: the stiff linear dispersion \(\partial_x^3\) is advanced exactly via the FFT, and the non-stiff advection \(6u\partial_x u\) via RK4.
Nonlinear Schrodinger – split-step Fourier: exact for the linear (kinetic) part and exact for the nonlinear (pointwise phase-rotation) part.
Sine-Gordon – explicit leapfrog finite differences on the wave equation.

physicskit.fields.solitons.kdv_evolve(u0, x, dt, steps)[source]#

Evolve a KdV initial condition forward in time on a periodic domain.

Parameters:

u0 (ndarray) – Initial field, sampled on the periodic grid x.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps (int) – Number of steps to advance.

Return type:

Returns:

ndarray – Field after steps * dt time units.

Notes

Because KdV solitons are true solitons, two colliding solitons of different speed pass through each other elastically: each emerges from the collision with its original amplitude and speed, shifted only in phase – unlike generic nonlinear waves, which would break up or change shape on collision.

Examples

A fast, tall soliton launched behind a slower, shorter one overtakes it; both retain their original amplitudes after the collision:

>>> import numpy as np
>>> from scipy.signal import find_peaks
>>> N, L = 512, 60.0
>>> x = np.linspace(-L / 2, L / 2, N, endpoint=False)
>>> u0 = kdv_soliton(x, c=9.0, x0=-20) + kdv_soliton(x, c=4.0, x0=-8)
>>> u = kdv_evolve(u0, x, dt=0.0005, steps=6000)
>>> peaks, _ = find_peaks(u, height=1.0)
>>> heights = sorted(u[peaks])
>>> [round(float(h), 1) for h in heights]
[2.0, 4.5]

physicskit.fields.solitons.kdv_evolve_frames(u0, x, dt, steps_per_frame, n_frames)[source]#

Evolve a KdV field and record a snapshot every steps_per_frame steps, for animation.

A thin looping wrapper around kdv_evolve() – no new physics, just repeated short calls with each chunk’s output fed back in as the next chunk’s initial condition (valid since kdv_evolve() is a deterministic one-step-at-a-time integrator), so the trajectory can be animated frame by frame.

Parameters:

u0 (ndarray) – Initial field.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps_per_frame (int) – Number of integration steps between recorded frames.
n_frames (int) – Number of frames to record after the initial condition.

Return type:

Returns:

frames (ndarray, shape (n_frames + 1, len(x))) – Field at t=0 and after each recorded chunk.
times (ndarray, shape (n_frames + 1,))

See also

physicskit.fields.visualizers.animate_field_1d: Render these frames as a line-plot animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-30, 30, 256, endpoint=False)
>>> u0 = kdv_soliton(x, c=4.0, x0=-15)
>>> frames, times = kdv_evolve_frames(u0, x, dt=0.001, steps_per_frame=200, n_frames=5)
>>> frames.shape
(6, 256)

physicskit.fields.solitons.kdv_soliton(x, c, x0=0.0)[source]#

Exact single-soliton solution of \(u_t + 6uu_x + u_{xxx} = 0\) at \(t=0\).

Parameters:

x (ndarray) – Spatial grid.
c (float) – Soliton speed (equal to twice its amplitude: amplitude \(=c/2\)).
x0 (float) – Initial center position.

Return type:

Returns:

ndarray – \(u(x, 0) = \tfrac{c}{2}\,\mathrm{sech}^2\!\big(\tfrac{\sqrt{c}}{2}(x-x_0)\big)\).

See also

kdv_evolve: Propagate this (or any) initial condition forward in time.

Examples

>>> import numpy as np
>>> round(float(kdv_soliton(0.0, c=4.0)), 6)
2.0

physicskit.fields.solitons.kdv_step(u_hat, k, dt)[source]#

Advance the KdV equation in Fourier space by one Strang-split step.

Parameters:

u_hat (ndarray) – Fourier coefficients of \(u\) (as from numpy.fft.fft).
k (ndarray) – Angular wavenumbers matching u_hat.
dt (float) – Time step.

Return type:

Returns:

ndarray – Fourier coefficients after one step of size dt.

physicskit.fields.solitons.nls_bright_soliton(x, t, A=1.0, v=0.0, x0=0.0)[source]#

Exact bright-soliton solution of the focusing NLS equation.

Solves \(i\psi_t + \tfrac{1}{2}\psi_{xx} + |\psi|^2\psi = 0\).

Parameters:

x (ndarray) – Spatial grid.
t (float) – Time.
A (float) – Soliton amplitude.
v (float) – Soliton velocity.
x0 (float) – Initial center position.

Return type:

Returns:

ndarray of complex – \(\psi(x,t) = A\,\mathrm{sech}(A(x-x_0-vt))\, e^{i[v(x-x_0) + (A^2-v^2)t/2]}\).

See also

nls_dark_soliton: The defocusing-equation counterpart.
nls_evolve: Propagate this (or any) initial condition numerically.

Examples

>>> import numpy as np
>>> x = np.array([0.0])
>>> abs(nls_bright_soliton(x, t=0.0, A=1.0))
array([1.])

physicskit.fields.solitons.nls_dark_soliton(x, t, rho0=1.0, x0=0.0)[source]#

Exact stationary dark (black) soliton solution of the defocusing NLS equation.

Solves \(i\psi_t + \tfrac{1}{2}\psi_{xx} - |\psi|^2\psi = 0\): a density notch that vanishes at its center, sitting on a uniform background of density rho0.

Parameters:

x (ndarray) – Spatial grid.
t (float) – Time.
rho0 (float) – Background density far from the soliton.
x0 (float) – Center position.

Return type:

Returns:

ndarray of complex – \(\psi(x,t) = \sqrt{\rho_0}\,\tanh\!\big(\sqrt{\rho_0}\,(x-x_0)\big)\, e^{-i\rho_0 t}\), with healing length \(\xi = 1/\sqrt{\rho_0}\).

See also

nls_bright_soliton: The focusing-equation counterpart.

Examples

>>> import numpy as np
>>> x = np.array([0.0])
>>> abs(nls_dark_soliton(x, t=0.0, rho0=1.0))
array([0.])

physicskit.fields.solitons.nls_evolve(psi0, x, dt, steps, g=1.0)[source]#

Evolve a nonlinear Schrodinger initial condition via split-step Fourier.

Solves \(i\psi_t + \tfrac{1}{2}\psi_{xx} + g|\psi|^2\psi = 0\) on a periodic domain: g > 0 is focusing (bright solitons), g < 0 is defocusing (dark solitons).

Parameters:

psi0 (ndarray) – Initial wavefunction, sampled on the periodic grid x.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps (int) – Number of steps to advance.
g (float) – Nonlinearity strength and sign.

Return type:

Returns:

ndarray of complex – Field after steps * dt time units.

Examples

A bright soliton propagates without dispersing – its envelope \(|\psi|\) is unchanged after 2 time units of evolution:

>>> import numpy as np
>>> N, L = 1024, 80.0
>>> x = np.linspace(-L / 2, L / 2, N, endpoint=False)
>>> psi0 = nls_bright_soliton(x, t=0.0, A=1.0)
>>> psi = nls_evolve(psi0, x, dt=0.001, steps=2000, g=1.0)
>>> shape_error = np.max(np.abs(np.abs(psi) - np.abs(psi0)))
>>> bool(shape_error < 1e-3)
True

physicskit.fields.solitons.nls_evolve_frames(psi0, x, dt, steps_per_frame, n_frames, g=1.0)[source]#

Evolve an NLS field and record a snapshot every steps_per_frame steps, for animation.

A thin looping wrapper around nls_evolve(), analogous to kdv_evolve_frames(): no new physics, just repeated chunked calls with the state carried forward, to produce an animatable sequence of frames.

Parameters:

psi0 (ndarray) – Initial wavefunction.
x (ndarray) – Uniformly spaced periodic spatial grid.
dt (float) – Time step.
steps_per_frame (int) – Number of integration steps between recorded frames.
n_frames (int) – Number of frames to record after the initial condition.
g (float) – Nonlinearity strength and sign.

Return type:

Returns:

frames (ndarray of complex, shape (n_frames + 1, len(x))) – Field at t=0 and after each recorded chunk.
times (ndarray, shape (n_frames + 1,))

See also

physicskit.fields.visualizers.animate_field_1d: Render these frames as a line-plot animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-40, 40, 512, endpoint=False)
>>> psi0 = nls_bright_soliton(x, t=0.0, A=1.0)
>>> frames, times = nls_evolve_frames(psi0, x, dt=0.001, steps_per_frame=200, n_frames=5, g=1.0)
>>> frames.shape
(6, 512)

physicskit.fields.solitons.sine_gordon_evolve(u0, u0_prev, x, dt, steps)[source]#

Evolve the Sine-Gordon equation with explicit leapfrog finite differences.

Parameters:

u0 (ndarray) – Field at the starting time.
u0_prev (ndarray) – Field one step before the starting time (needed to seed the two-level leapfrog scheme); for a traveling-wave initial condition, use the exact solution evaluated at t=-dt.
x (ndarray) – Uniformly spaced spatial grid, with fixed (Dirichlet) boundaries.
dt (float) – Time step; must satisfy the CFL condition dt <= dx (wave speed 1).
steps (int) – Number of steps to advance.

Return type:

Returns:

u, u_prev (ndarray) – The field at the final step and the step before it (the pair needed to continue the integration further).

Examples

A kink launched at speed \(v=0.5\) arrives at the expected Lorentz-contracted position after propagating:

>>> import numpy as np
>>> N, L, v, x0 = 4000, 200.0, 0.5, -50.0
>>> x = np.linspace(-L / 2, L / 2, N)
>>> dx = x[1] - x[0]
>>> dt = 0.4 * dx
>>> u_prev = sine_gordon_kink(x, -dt, v, x0)
>>> u0 = sine_gordon_kink(x, 0.0, v, x0)
>>> steps = 1000
>>> u, _ = sine_gordon_evolve(u0, u_prev, x, dt, steps)
>>> expected = sine_gordon_kink(x, steps * dt, v, x0)
>>> bool(np.max(np.abs(u - expected)) < 0.01)
True

physicskit.fields.solitons.sine_gordon_evolve_frames(u0, u0_prev, x, dt, steps_per_frame, n_frames)[source]#

Evolve a Sine-Gordon field and record a snapshot every steps_per_frame steps, for animation.

A thin looping wrapper around sine_gordon_evolve(), analogous to kdv_evolve_frames(): the (u, u_prev) leapfrog pair is carried forward chunk to chunk, with no change to the underlying physics.

Parameters:

u0 (ndarray) – Field at the starting time.
u0_prev (ndarray) – Field one step before the starting time (see sine_gordon_evolve()).
x (ndarray) – Uniformly spaced spatial grid.
dt (float) – Time step.
steps_per_frame (int) – Number of integration steps between recorded frames.
n_frames (int) – Number of frames to record after the initial condition.

Return type:

Returns:

frames (ndarray, shape (n_frames + 1, len(x))) – Field at t=0 and after each recorded chunk.
times (ndarray, shape (n_frames + 1,))

See also

physicskit.fields.visualizers.animate_field_1d: Render these frames as a line-plot animation.

Examples

>>> import numpy as np
>>> x = np.linspace(-50, 50, 800)
>>> dt = 0.4 * (x[1] - x[0])
>>> u_prev = sine_gordon_kink(x, -dt, v=0.5, x0=-20)
>>> u0 = sine_gordon_kink(x, 0.0, v=0.5, x0=-20)
>>> frames, times = sine_gordon_evolve_frames(u0, u_prev, x, dt, steps_per_frame=100, n_frames=4)
>>> frames.shape
(5, 800)

physicskit.fields.solitons.sine_gordon_kink(x, t, v=0.0, x0=0.0, polarity=1)[source]#

Exact kink (or antikink) solution of the Sine-Gordon equation \(u_{tt} - u_{xx} + \sin u = 0\).

Parameters:

x (ndarray) – Spatial grid.
t (float) – Time.
v (float) – Kink velocity (\(|v| < 1\), the wave speed of the linearized equation).
x0 (float) – Initial center position.
polarity (int) – 1 for a kink (field jumps by \(2\pi\)), -1 for an antikink.

Return type:

Returns:

ndarray – \(u(x,t) = 4\arctan\!\big(\exp[\text{polarity}\cdot\gamma(x-x_0-vt)]\big)\), with \(\gamma = 1/\sqrt{1-v^2}\).

See also

sine_gordon_evolve: Propagate this (or any) initial condition numerically.

Examples

>>> import numpy as np
>>> round(float(sine_gordon_kink(0.0, t=0.0)), 6)
3.141593

Gross-Pitaevskii mean-field theory for trapped, rotating Bose-Einstein condensates.

Solves the time-dependent Gross-Pitaevskii equation (GPE)

\[i\partial_t\psi = \Big[-\tfrac{1}{2}\nabla^2 + V(\mathbf{r}) + g|\psi|^2 - \Omega L_z\Big]\psi\]

in the rotating frame (units \(\hbar=m=1\)) by imaginary-time propagation (\(\tau=it\)), which turns the Schrodinger-like evolution into a gradient flow that relaxes any initial state toward a stationary point of the rotating-frame energy. The kinetic and potential/interaction terms are each handled exactly in their natural representation (spectral and real-space respectively, split-step); the angular-momentum term \(\Omega L_z\) is advanced explicitly via centered finite differences, which is non-stiff at the small imaginary-time steps used here.

A single quantized vortex, imprinted by multiplying the wavefunction by \((x-x_0)+i(y-y_0)\), is used throughout as the worked example of a “vortex lattice” building block: gpe_relax() finds it as a genuine stationary GPE solution, and comparing its rotating-frame energy against the vortex-free ground state reproduces the standard textbook criterion for the critical rotation frequency \(\Omega_c\) above which nucleating a vortex lowers the energy.

physicskit.fields.quantum_fields.casimir_energy_1d(d, c=1.0, cutoff=None)[source]#

Regularized zero-point (Casimir) energy of the discrete 1D cavity spectrum.

Uses the standard exponential-cutoff regularization: the sum \(\sum_n n x^n = x/(1-x)^2\) (with \(x=e^{-a}\), \(a=\pi c\,\text{cutoff}/d\)) has an exact closed form, so no series truncation is needed. Its small-cutoff (small-a) expansion is \(\tfrac{1}{4\sinh^2(a/2)} = 1/a^2 - 1/12 + O(a^2)\); the \(1/a^2\) piece is the (unphysical, cutoff-scheme-dependent) divergence that a continuum reference calculation regularized the same way would also produce, so it is subtracted exactly, leaving the finite remainder that survives as cutoff -> 0:

\[E(d) \to -\frac{\pi c}{24 d}\]

the standard 1D massless-field Casimir energy (equivalently, the Casimir energy of a CFT strip with central charge 1).

Parameters:

d (float) – Plate separation.
c (float) – Wave speed.
cutoff (float | None) – Regulator scale. Defaults to 1e-3 * d / c. This subtraction is a difference of two large, nearly-equal floating-point terms (\(\propto 1/\text{cutoff}^2\)), so making cutoff too small loses precision to cancellation rather than gaining it; 1e-4*d/c to 1e-2*d/c is the well-behaved range in double precision.

Return type:

Returns:

float – The regularized (finite, cutoff-subtracted) zero-point energy.

See also

casimir_mode_frequencies: The (bare, un-regularized) mode spectrum being summed.

Examples

The cutoff-and-subtract result agrees with the known closed form \(-\pi c/(24d)\), and (within the well-behaved cutoff range noted above) agrees more closely as the cutoff shrinks:

>>> d = 3.0
>>> exact = -np.pi * 1.0 / (24 * d)
>>> bool(abs(casimir_energy_1d(d, cutoff=3e-3 * d) - exact) < 1e-6)
True
>>> bool(abs(casimir_energy_1d(d, cutoff=7e-4 * d) - exact) < 1e-8)
True

The energy grows less negative (weaker confinement of vacuum energy) as the plates separate, giving an attractive force \(-dE/dd < 0\):

>>> bool(casimir_energy_1d(1.0) < casimir_energy_1d(2.0) < 0)
True

physicskit.fields.quantum_fields.casimir_mode_frequencies(d, c=1.0, n_max=200)[source]#

Discrete standing-wave mode frequencies of a 1D cavity of plate separation d.

Parameters:

d (float) – Plate separation.
c (float) – Wave speed (=1 in natural units; use the physical speed of light for SI-unit frequencies).
n_max (int) – Number of modes to return.

Return type:

Returns:

ndarray, shape (n_max,) – \(\omega_n = n\pi c/d\) for \(n=1,\dots,n_{max}\).

See also

casimir_energy_1d: The regularized zero-point energy of this mode spectrum.

Examples

>>> omega = casimir_mode_frequencies(d=2.0, c=1.0, n_max=3)
>>> [round(float(w), 4) for w in omega]
[1.5708, 3.1416, 4.7124]

physicskit.fields.quantum_fields.count_vortices(psi, density_threshold=0.05)[source]#

Locate quantized vortices by summing the phase winding around each grid plaquette.

Plaquettes in very low density regions are excluded, since the phase of a near-zero-amplitude wavefunction is dominated by numerical noise and produces spurious windings there.

A vortex core that sits exactly on a grid vertex (rather than strictly inside a plaquette) can go undetected, since its circulation is then split ambiguously between the four plaquettes touching that vertex; in practice this is avoided by seeding vortices at positions not exactly on the grid.

Parameters:

psi (ndarray) – Wavefunction.
density_threshold (float) – Plaquettes where \(|\psi|^2\) (relative to its maximum) falls below this fraction are excluded from the search.

Return type:

Returns:

ndarray of int, shape (n, n) – Integer winding number around each plaquette’s corner (i, j); nonzero entries mark a vortex core (+1 or -1 for a singly-quantized vortex/antivortex) enclosed by that plaquette.

Examples

>>> import numpy as np
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(48, 10.0)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi_vortex = gpe_imprint_vortex(psi0, X, Y, [(0.0, 0.0)])
>>> winding = count_vortices(psi_vortex)
>>> int(np.sum(np.abs(winding)))
1

physicskit.fields.quantum_fields.gpe_energy(psi, V, g, X, Y, K2)[source]#

Evaluate the energy and angular momentum of a GPE wavefunction.

Parameters:

psi (ndarray) – Wavefunction.
V (ndarray) – External trapping potential.
g (float) – Interaction strength.
X (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
Y (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
K2 (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).

Return type:

dict

Returns:

dict – {"kinetic", "potential", "interaction", "angular_momentum", "total"}, where total is the lab-frame energy (kinetic + potential + interaction) and angular_momentum is \(\langle L_z \rangle\). The rotating-frame energy at rotation rate Omega is total - Omega * angular_momentum.

Examples

>>> import numpy as np
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(48, 10.0)
>>> V = 0.5 * (X ** 2 + Y ** 2)
>>> psi = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi *= 1.0 / np.sqrt(np.sum(np.abs(psi) ** 2) * (X[1, 0] - X[0, 0]) ** 2)
>>> E = gpe_energy(psi, V, g=0.0, X=X, Y=Y, K2=K2)
>>> round(float(E["kinetic"] + E["potential"]), 4)
1.0

physicskit.fields.quantum_fields.gpe_evolve(psi0, V, g, dt, steps, K2, snapshot_stride=1)[source]#

Real-time propagation of the 2D Gross-Pitaevskii / cubic NLS equation via split-step Fourier.

Solves \(i\partial_t\psi = [-\tfrac{1}{2}\nabla^2 + V + g|\psi|^2]\psi\) forward in real time (unlike gpe_relax()’s imaginary-time relaxation, which only finds stationary states): the kinetic term is exact in Fourier space, and the potential-plus-nonlinear term is exact as a pointwise phase rotation in real space, the same split-step structure as physicskit.fields.solitons.nls_evolve() generalized to 2D and to an external potential. Real-time evolution is unitary and conserves the norm on its own (no renormalization needed, unlike imaginary time).

This one stepper serves two purposes, distinguished only by V and the sign of g:

With a harmonic trap (V = 0.5*(X**2+Y**2)) and repulsive interactions (g > 0), an off-center vortex genuinely precesses around the trap under the density gradient – real vortex dynamics, as opposed to the static relaxed state gpe_relax() finds.
With V = 0 (or weak) and attractive interactions (g < 0 here, opposite sign convention from nls_evolve’s g since this module’s equation carries a +g|psi|^2 term), a sufficiently tall, narrow initial packet self-focuses, concentrating into a narrower, taller peak – a numerical stand-in for wave collapse. True collapse is a singularity in finite time; this integrator (like any finite-grid scheme) cannot resolve it and the calculation should be stopped once the peak density is still visibly growing, not carried through the blow-up itself.

Parameters:

psi0 (ndarray) – Initial wavefunction.
V (ndarray) – External potential (use np.zeros_like for the free/collapse case).
g (float) – Interaction strength and sign (see above).
dt (float) – Time step.
steps (int) – Number of steps to advance.
K2 (ndarray) – Wavenumber-squared grid from harmonic_trap_grid().
snapshot_stride (int) – Record a snapshot every this many steps (plus the initial condition).

Return type:

Returns:

snapshots (ndarray of complex, shape (n_recorded, n, n)) – Wavefunction at t=0 and after every recorded step.
times (ndarray, shape (n_recorded,)) – Time of each recorded snapshot.

See also

gpe_relax: Imaginary-time relaxation to a stationary state.
physicskit.fields.solitons.nls_evolve: The 1D analog (opposite sign convention for g).

Examples

Norm is conserved under real-time evolution, unlike the imaginary-time propagation in gpe_relax() (which requires explicit renormalization):

>>> import numpy as np
>>> n, length = 48, 12.0
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(n, length)
>>> dxg = X[1, 0] - X[0, 0]
>>> V = 0.5 * (X ** 2 + Y ** 2)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi0 *= 1.0 / np.sqrt(np.sum(np.abs(psi0) ** 2) * dxg * dxg)
>>> snaps, times = gpe_evolve(psi0, V, g=2.0, dt=1e-3, steps=200, K2=K2, snapshot_stride=50)
>>> norms = np.sum(np.abs(snaps) ** 2, axis=(1, 2)) * dxg * dxg
>>> bool(np.max(np.abs(norms - norms[0])) < 1e-6)
True

physicskit.fields.quantum_fields.gpe_imprint_vortex(psi, X, Y, positions)[source]#

Imprint one singly-quantized vortex per position by multiplying in a phase winding.

Parameters:

psi (ndarray) – Wavefunction to imprint onto.
X (ndarray) – Real-space coordinate grids from harmonic_trap_grid().
Y (ndarray) – Real-space coordinate grids from harmonic_trap_grid().
positions (sequence of tuple(float, float)) – (x0, y0) center of each vortex to imprint.

Return type:

Returns:

ndarray of complex, shape (n, n) – psi multiplied by \(\prod_i [(x-x_{0,i}) + i(y-y_{0,i})]\), renormalized to preserve the input norm.

Examples

>>> import numpy as np
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(32, 12.0)
>>> psi0 = np.exp(-0.5 * (X ** 2 + Y ** 2)).astype(complex)
>>> psi = gpe_imprint_vortex(psi0, X, Y, [(0.0, 0.0)])
>>> bool(abs(psi[16, 16]) < 1e-10)
True

physicskit.fields.quantum_fields.gpe_relax(psi0, V, g, dtau, steps, X, Y, K2, Omega=0.0, n_particles=1.0)[source]#

Relax a GPE initial condition toward a stationary state via imaginary-time propagation.

Parameters:

psi0 (ndarray) – Initial wavefunction.
V (ndarray) – External trapping potential.
g (float) – Interaction (nonlinearity) strength.
dtau (float) – Imaginary-time step.
steps (int) – Number of steps to advance.
X (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
Y (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
K2 (ndarray) – Grids from harmonic_trap_grid() (KX, KY are not needed here).
Omega (float) – Rotation frequency of the trap.
n_particles (float) – Total particle number; the wavefunction is renormalized to this value after every step (imaginary-time propagation does not conserve norm on its own).

Return type:

Returns:

ndarray of complex, shape (n, n) – The relaxed (approximately stationary) wavefunction.

See also

gpe_energy: Evaluate the energy of a relaxed state.
gpe_imprint_vortex: Seed an initial condition with quantized vortices.

Notes

Energy decreases monotonically under this propagation whenever it is implemented correctly (a basic sanity check worth verifying on any new potential or parameter regime).

Examples

Comparing the rotating-frame energy of the vortex-free ground state against a state seeded with one centered vortex reproduces the standard vortex-nucleation criterion: \(\Omega_c = \Delta E/\Delta L_z\), above which the vortex state has lower energy in the rotating frame:

>>> import numpy as np
>>> n, length, g = 64, 12.0, 4.0
>>> X, Y, KX, KY, K2 = harmonic_trap_grid(n, length)
>>> V = 0.5 * (X ** 2 + Y ** 2)
>>> psi_vf = gpe_relax(np.exp(-0.5 * (X**2 + Y**2)).astype(complex), V, g, dtau=5e-4, steps=4000, X=X, Y=Y, K2=K2)
>>> psi0_v = gpe_imprint_vortex(np.exp(-0.5 * (X**2 + Y**2)).astype(complex), X, Y, [(0.0, 0.0)])
>>> psi_v = gpe_relax(psi0_v, V, g, dtau=5e-4, steps=4000, X=X, Y=Y, K2=K2)
>>> E0, E1 = gpe_energy(psi_vf, V, g, X, Y, K2), gpe_energy(psi_v, V, g, X, Y, K2)
>>> Omega_c = (E1["total"] - E0["total"]) / (E1["angular_momentum"] - E0["angular_momentum"])
>>> bool(0.0 < Omega_c < 1.0)
True

physicskit.fields.quantum_fields.harmonic_trap_grid(n, length)[source]#

Build a centered real-space and wavenumber grid for a 2D harmonic trap.

Parameters:

n (int) – Number of grid points along each axis.
length (float) – Physical domain size (the domain is [-length/2, length/2)).

Return type: