1 files changed, 112 insertions, 0 deletions
diff --git a/src/physics/models.ts b/src/physics/models.ts
new file mode 100644
index 0000000..8fb121f
--- /dev/null
+++ b/src/physics/models.ts
@@ -0,0 +1,112 @@
+import { COULOMB_CONSTANT, electricField, type Charge, type Point } from './electricField';
+
+export const VACUUM_PERMITTIVITY = 8.8541878128e-12;
+export const ELEMENTARY_CHARGE = 1.602176634e-19;
+const THERMAL_ENERGY_EV = 8.617333262e-5 * 300;
+
+function positive(...values: number[]) {
+  if (values.some(value => !Number.isFinite(value) || value <= 0)) {
+    throw new RangeError('Physical dimensions and material constants must be finite and positive.');
+  }
+}
+
+/** Vacuum potential in volts, referenced to infinity. Undefined inside the
+ * same 0.18 m visualization exclusion used by electricField. */
+export function potential(point: Point, charges: readonly Charge[]): number {
+  let volts = 0;
+  for (const charge of charges) {
+    if (charge.coulombs === 0) continue;
+    const r = Math.hypot(point.x - charge.x, point.y - charge.y, point.z - charge.z);
+    if (r < 0.18) return NaN;
+    volts += COULOMB_CONSTANT * charge.coulombs / r;
+  }
+  return volts;
+}
+
+/** Force on a test charge, excluding that charge from the supplied sources. */
+export function forceOn(test: Charge, sources: readonly Charge[]): Point {
+  const field = electricField(test, sources);
+  return { x: field.x * test.coulombs, y: field.y * test.coulombs, z: field.z * test.coulombs };
+}
+
+export function pairEnergy(q1: number, q2: number, separation: number): number {
+  positive(separation);
+  return COULOMB_CONSTANT * q1 * q2 / separation;
+}
+
+/** Exact signed flux for an origin-centred spherical Gaussian surface.
+ * A point charge exactly on the surface makes the integral undefined. */
+export function gaussianFlux(charges: readonly Charge[], radius: number): number {
+  positive(radius);
+  let enclosed = 0;
+  for (const charge of charges) {
+    if (charge.coulombs === 0) continue;
+    const r = Math.hypot(charge.x, charge.y, charge.z);
+    if (Math.abs(r - radius) < 1e-9) return NaN;
+    if (r < radius) enclosed += charge.coulombs;
+  }
+  return enclosed / VACUUM_PERMITTIVITY;
+}
+
+export interface CapacitorParameters {
+  area: number;
+  gap: number;
+  relativePermittivity: number;
+  voltage: number;
+  /** If isolated, voltage means the initial vacuum charging voltage. */
+  isolated?: boolean;
+}
+export function capacitor(p: CapacitorParameters) {
+  positive(p.area, p.gap, p.relativePermittivity);
+  const vacuumCapacitance = VACUUM_PERMITTIVITY * p.area / p.gap;
+  const capacitance = vacuumCapacitance * p.relativePermittivity;
+  const voltage = p.isolated ? p.voltage / p.relativePermittivity : p.voltage;
+  const charge = capacitance * voltage;
+  const field = voltage / p.gap;
+  return {
+    capacitance, voltage, charge, field,
+    energy: 0.5 * capacitance * voltage ** 2,
+    polarization: VACUUM_PERMITTIVITY * (p.relativePermittivity - 1) * field,
+  };
+}
+
+export interface ConductorParameters {
+  voltage: number;
+  resistivity: number;
+  length: number;
+  area: number;
+  carrierDensity: number;
+}
+export function conductor(p: ConductorParameters) {
+  positive(p.resistivity, p.length, p.area, p.carrierDensity);
+  const resistance = p.resistivity * p.length / p.area;
+  const current = p.voltage / resistance;
+  return {
+    resistance, current,
+    power: p.voltage * current,
+    density: current / p.area,
+    electronDrift: -current / (p.carrierDensity * ELEMENTARY_CHARGE * p.area),
+  };
+}
+
+export interface SemiconductorParameters {
+  intrinsicDensity: number;
+  /** Fully ionized donors minus acceptors, in m^-3. */
+  netDoping: number;
+  electronMobility: number;
+  holeMobility: number;
+}
+/** Non-degenerate thermal equilibrium at 300 K. Solving for the minority
+ * through mass action avoids catastrophic subtraction at large doping. */
+export function semiconductor(p: SemiconductorParameters) {
+  positive(p.intrinsicDensity, p.electronMobility, p.holeMobility);
+  const majority = (Math.abs(p.netDoping) + Math.hypot(p.netDoping, 2 * p.intrinsicDensity)) / 2;
+  const minority = p.intrinsicDensity ** 2 / majority;
+  const electrons = p.netDoping >= 0 ? majority : minority;
+  const holes = p.netDoping >= 0 ? minority : majority;
+  return {
+    electrons, holes,
+    fermiShift: THERMAL_ENERGY_EV * Math.log(electrons / p.intrinsicDensity),
+    conductivity: ELEMENTARY_CHARGE * (electrons * p.electronMobility + holes * p.holeMobility),
+  };
+}