Transforming the field equations

export function maxwellResidualsPlaneWave(input: {
  beta: number;
  wave: PlaneWaveKind;
  polarization: Polarization;
  events?: readonly Event4[];
  convention?: TransformConvention;
  E0?: number;
  omega?: number;
}): Readonly<{
  maxResidual: number;
  residuals: Float64Array;
  amplitudeFactor: number;
  frequencyFactor: number;
  gamma: number;
  passed: boolean;
  eventSeed: string;
}> {
  const v = input.beta;
  const g = gamma(v);
  if (g.status !== "value") {
    return Object.freeze({
      maxResidual: Number.POSITIVE_INFINITY,
      residuals: new Float64Array(6),
      amplitudeFactor: Number.NaN,
      frequencyFactor: Number.NaN,
      gamma: Number.NaN,
      passed: false,
      eventSeed: SR07_EVENT_SEED,
    });
  }
  const E0 = input.E0 ?? 1;
  const omega = input.omega ?? 1;
  const convention = input.convention ?? "printed";
  const n = waveNormal(input.wave);
  const eHat = wavePolarization(input.wave, input.polarization);
  const events = input.events ?? seededPlaneWaveEvents();
  const residuals = new Float64Array(6);
  let maxResidual = 0;
  const γ = g.value;
  const doppler = γ * (1 - n.x * v);
  const amplitudeFactor = input.wave === "plus-x" || input.wave === "minus-x" ? doppler : γ;
  const frequencyFactor = doppler;
  let allPassed = true;
  for (const event of events) {
    const w = planeWaveAt(n, eHat, E0, omega, event);
    const dEdt = applyPrintedOrFlipped(w.dEdt, w.dBdt, v, convention);
    const dEdx = applyPrintedOrFlipped(w.dEdx, w.dBdx, v, convention);
    const dEdy = applyPrintedOrFlipped(w.dEdy, w.dBdy, v, convention);
    const dEdz = applyPrintedOrFlipped(w.dEdz, w.dBdz, v, convention);
    const eTau = chainPrime(dEdt.E, dEdx.E, dEdy.E, dEdz.E, γ, v);
    const bTau = chainPrime(dEdt.B, dEdx.B, dEdy.B, dEdz.B, γ, v);
    const curlE: readonly [number, number, number] = [
      eTau.dEta.z - eTau.dZeta.y,
      eTau.dZeta.x - eTau.dXi.z,
      eTau.dXi.y - eTau.dEta.x,
    ];
    const curlB: readonly [number, number, number] = [
      bTau.dEta.z - bTau.dZeta.y,
      bTau.dZeta.x - bTau.dXi.z,
      bTau.dXi.y - bTau.dEta.x,
    ];
    const faraday: readonly [number, number, number] = [
      curlE[0] + bTau.dTau.x,
      curlE[1] + bTau.dTau.y,
      curlE[2] + bTau.dTau.z,
    ];
    const ampere: readonly [number, number, number] = [
      curlB[0] - eTau.dTau.x,
      curlB[1] - eTau.dTau.y,
      curlB[2] - eTau.dTau.z,
    ];
    const six = [...faraday, ...ampere];
    for (let i = 0; i < 6; i++) {
      const axis = i % 3 === 0 ? "x" : i % 3 === 1 ? "y" : "z";
      const terms =
        i < 3
          ? [curlE[i] ?? 0, bTau.dTau[axis], E0 * omega]
          : [curlB[i % 3] ?? 0, eTau.dTau[axis], E0 * omega];
      const tol = SR07_RESIDUAL_RELATIVE * maxAbs(terms) + SR07_RESIDUAL_FLOOR * E0 * omega;
      const r = Math.abs(six[i] ?? 0);
      if (r > (residuals[i] ?? 0)) residuals[i] = r;
      if (r > maxResidual) maxResidual = r;
      if (r > tol) allPassed = false;
    }
  }
  return Object.freeze({
    maxResidual,
    residuals,
    amplitudeFactor,
    frequencyFactor,
    gamma: γ,
    passed: allPassed,
    eventSeed: SR07_EVENT_SEED,
  });
}

export function transformGaussianHistorical(input: {
  E: Vec3;
  B: Vec3;
  boost: number;
  c?: number;
}): Readonly<{ E: Vec3; B: Vec3; beta: number }> {
  const c = input.c ?? C_SI;
  const v = input.boost;
  const g = gamma(v / c);
  if (g.status !== "value") {
    return Object.freeze({ E: input.E, B: input.B, beta: Number.NaN });
  }
  const β = g.value;
  const { E, B } = input;
  const vOverC = v / c;
  return Object.freeze({
    beta: β,
    E: Object.freeze({
      x: E.x,
      y: β * (E.y - vOverC * B.z),
      z: β * (E.z + vOverC * B.y),
    }),
    B: Object.freeze({
      x: B.x,
      y: β * (B.y + vOverC * E.z),
      z: β * (B.z - vOverC * E.y),
    }),
  });
}