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),
  };
}