Vector Field:
Stochastic Flows & Real-time Integration

A vector field is a mapping that associates a vector to every point in a coordinate space. Mathematically, in an $n$-dimensional Euclidean space, a vector field is represented as a function:

Vector fields are the core mathematical framework for modeling physical phenomena where quantities vary in both direction and magnitude across space. Examples include the velocity of fluid particles, electrostatic and gravitational forces, atmospheric wind currents, and the phase portraits of non-linear dynamical systems.

Visualizing these fields is a primary challenge in scientific computing. Rather than drawing static arrays of discrete arrows—which clutter regions of high intensity and fail to convey global structure—we compute continuous **streamlines** (integral curves) that follow the field's tangents.

Our implementation leverages **Macroquad** for fast rendering, **Egui** for immediate-mode GUI control, and **Glam** for vector operations. The engine includes:

• 7 Mathematical Presets: Incorporates classic fields (Vortex, Electric Dipole, Saddle) alongside noise-based flow models (Perlin Curl, Perlin Gradient, Trig Flow) and interactive N-Body orbital mechanics.
• Dual Physics Models: Supports **Velocity Flow** ($v = \mathbf{F}(\mathbf{x})$) where particles trace fields directly, and **Force Field** ($a = \mathbf{F}(\mathbf{x})/m$) enabling inertial motion, orbits, and slingshots.
• Dynamic Viewport Interactivity: Allows users to perturb field coordinates locally via attraction, repulsion, vortices, or by spawning, dragging, and deleting massive gravitational bodies.
• Tapered Alpha-Faded Particle Trails: Implements alpha-decay buffers mapping historical coordinates onto screen-space vectors for high-fidelity visualization.

Analyzing a vector field $\mathbf{F} = (P, Q, R)$ requires examining its local derivatives. The behavior of a vector field at any given point is governed by two fundamental operators in vector calculus: **Divergence** and **Curl**.

Divergence ($\nabla \cdot \mathbf{F}$)

Divergence measures the net outward flux of a vector field per unit volume from an infinitesimal boundary around a point. It quantifies whether a point acts as a **source** ($\nabla \cdot \mathbf{F} > 0$) or a **sink** ($\nabla \cdot \mathbf{F} < 0$). In Cartesian coordinates:

Curl ($\nabla \times \mathbf{F}$)

Curl measures the rotation or vorticity of a vector field about a point. It yields a vector describing the axis and rate of rotation. In 3D:

In a two-dimensional vector field $\mathbf{F}(x, y) = (P, Q)$, the curl is simplified to a scalar component acting orthogonal to the coordinate plane:

Our simulation environment supports two fundamentally distinct physical models to translate vector values into coordinate displacements:

1. Velocity Flow (First-order system): Particle velocity maps directly to the vector field ($\mathbf{v} = \mathbf{F}(\mathbf{x})$). The pathlines follow streamlines exactly: $$\frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x})$$
2. Force Field (Second-order system with Inertia): The vector field acts as an accelerating force ($\mathbf{a} = \mathbf{F}(\mathbf{x})/m$). Particles carry inertia, enabling orbits, chaotic loops, and slingshot physics: $$\frac{d^2\mathbf{x}}{dt^2} + \gamma \frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x})$$ *(where $\gamma$ is a damping coefficient to prevent infinite energy accumulation).*

Numerical Integrators: Euler vs. Runge-Kutta 4 (RK4)

To approximate these paths, the engine evaluates coordinates step-by-step. Standard **Euler integration** ($\mathbf{x}_{n+1} = \mathbf{x}_n + \mathbf{F}(\mathbf{x}_n)\Delta t$) is computationally trivial but accumulates massive truncation errors $\mathcal{O}(\Delta t)$.

For stable orbits and high-vorticity zones, the sandbox implements a **Runge-Kutta 4th Order (RK4)** scheme:

Below is a complete Python implementation for calculating vector field streamlines and numerical Jacobians using `numpy` and `matplotlib`:

import numpy as np
import matplotlib.pyplot as plt

def compute_jacobian(F, x, y, dx=1e-5):
    """
    Computes the numerical Jacobian matrix of a 2D vector field F at (x, y)
    using central differences.
    """
    Fx_x, Fy_x = F(x + dx, y)
    Fx_x_neg, Fy_x_neg = F(x - dx, y)
    
    Fx_y, Fy_y = F(x, y + dx)
    Fx_y_neg, Fy_y_neg = F(x, y - dx)
    
    dFx_dx = (Fx_x - Fx_x_neg) / (2 * dx)
    dFy_dx = (Fy_x - Fy_x_neg) / (2 * dx)
    dFx_dy = (Fx_y - Fx_y_neg) / (2 * dx)
    dFy_dy = (Fy_y - Fy_y_neg) / (2 * dx)
    
    jacobian = np.array([
        [dFx_dx, dFx_dy],
        [dFy_dx, dFy_dy]
    ])
    return jacobian

# Define a vortex preset field
def vortex_field(x, y):
    return np.array([-y, x])

# Plotting Streamlines
x, y = np.meshgrid(np.linspace(-3, 3, 20), np.linspace(-3, 3, 20))
u, v = vortex_field(x, y)

plt.figure(figsize=(6, 6), facecolor='#09090b')
ax = plt.axes()
ax.set_facecolor('#09090b')
ax.streamplot(x, y, u, v, color='#38bdf8', linewidth=1, arrowsize=1.2)
ax.tick_params(colors='#71717a')
plt.title("Vector Field Streamlines", color='#fafafa')
plt.show()