physicskit.fields #

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.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.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.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.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.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

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: