Charge and current density

export function evaluateSr12(input: Sr12Input): Sr12Snapshot {
  const c = C_SI;
  const v = input.boost;
  const eventId = "event-charge-current-0";

  const packScalar = (
    key: keyof typeof SR12_ROWS,
    make: (id: string, unit: string, kind: string, owner: string) => ScientificResult,
  ) => {
    const r = SR12_ROWS[key];
    return make(r[0], r[1], r[2], r[3]);
  };

  const valScalar = (key: keyof typeof SR12_ROWS, num: number) =>
    packScalar(key, (id, unit, kind, owner) => asValue(id, unit, kind, owner, num));

  const valVector = (key: keyof typeof SR12_ROWS, vec: Vec3) =>
    packScalar(key, (id, unit, kind, owner) => asVectorValue(id, unit, kind, owner, vec));

  const refuseAll = (reason: string, condition: string): Sr12Snapshot => {
    const out = (key: keyof typeof SR12_ROWS) =>
      packScalar(key, (id, unit, kind, owner) =>
        asOutside(id, unit, kind, owner, condition, reason),
      );
    return Object.freeze({
      beta: v / c,
      gamma: Number.NaN,
      chargeDensityStationary: out("rhoM"),
      chargeDensityMoving: out("rhoC"),
      currentDensityStationary: out("JM"),
      currentDensityMoving: out("JC"),
      fourCurrentInvariant: out("inv"),
      fourCurrentInvariantNormalized: out("invNorm"),
      lorentzFactor: out("gamma"),
      continuityResidualStationary: out("contM"),
      continuityResidualMoving: out("contC"),
      loopLegChargePositive: out("legPos"),
      loopLegChargeNegative: out("legNeg"),
      loopTotalCharge: out("loopTotal"),
      sphereTotalChargeStationary: out("sphereM"),
      sphereTotalChargeMoving: out("sphereC"),
      eventId,
    });
  };

  if (
    ![
      input.boost,
      input.chargeDensity,
      input.currentDensity.x,
      input.currentDensity.y,
      input.currentDensity.z,
      input.carrierVelocity.x,
      input.carrierVelocity.y,
      input.carrierVelocity.z,
      input.sphereRadius,
      input.sphereCharge,
      input.loopCurrent,
      input.loopLengthX,
      input.loopLengthY,
      input.pulseWidth,
      input.pulseAmplitude,
    ].every(Number.isFinite)
  ) {
    return refuseAll(
      "All density, current, velocity and boost inputs must be finite numbers.",
      "finite inputs",
    );
  }

  const beta = v / c;
  if (Math.abs(beta) > 0.95) {
    return refuseAll("Observer boost speed must satisfy |v| <= 0.95c.", "|v| <= 0.95c");
  }

  const g = gamma(beta);
  if (g.status !== "value") {
    return refuseAll("Lorentz factor could not be evaluated.", "|v| < c");
  }
  const γ = g.value;

  const rho0 = input.chargeDensity;
  let J0: Vec3 = input.currentDensity;

  if (input.mode === "convection") {
    const uMagSq =
      input.carrierVelocity.x ** 2 + input.carrierVelocity.y ** 2 + input.carrierVelocity.z ** 2;
    if (uMagSq >= c * c) {
      return refuseAll("Carrier velocity must satisfy |u| < c.", "|u| < c");
    }
    J0 = Object.freeze({
      x: rho0 * input.carrierVelocity.x,
      y: rho0 * input.carrierVelocity.y,
      z: rho0 * input.carrierVelocity.z,
    });
  }

  const transformed = transformChargeCurrent({ rho: rho0, J: J0, boost: v, c });
  const invs = fourCurrentInvariants(rho0, J0, c);

  // Mode-specific calculations
  const loop = currentLoopCharges(input.loopCurrent, input.loopLengthX, input.loopLengthY, v, c);

  const sphereQ0 = sphereTotalCharge(input.sphereRadius, rho0);
  const sphereQPrime = sphereQ0; // Exact charge invariance

  const pulse = gaussianPulseContinuity(
    input.pulseAmplitude,
    input.carrierVelocity.x,
    input.pulseWidth,
    0,
    0,
    v,
    c,
  );

  return Object.freeze({
    beta,
    gamma: γ,
    chargeDensityStationary: valScalar("rhoM", rho0),
    chargeDensityMoving: valScalar("rhoC", transformed.rho),
    currentDensityStationary: valVector("JM", J0),
    currentDensityMoving: valVector("JC", transformed.J),
    fourCurrentInvariant: valScalar("inv", invs.si),
    fourCurrentInvariantNormalized: valScalar(
      "invNorm",
      (c * rho0) ** 2 / (c * c) - (J0.x ** 2 + J0.y ** 2 + J0.z ** 2) / (c * c),
    ),
    lorentzFactor: valScalar("gamma", γ),
    continuityResidualStationary: valScalar("contM", pulse.stationaryResidual || 0),
    continuityResidualMoving: valScalar("contC", pulse.movingResidual || 0),
    loopLegChargePositive: valScalar("legPos", loop.legChargePositive),
    loopLegChargeNegative: valScalar("legNeg", loop.legChargeNegative),
    loopTotalCharge: valScalar("loopTotal", loop.totalCharge),
    sphereTotalChargeStationary: valScalar("sphereM", sphereQ0),
    sphereTotalChargeMoving: valScalar("sphereC", sphereQPrime),
    eventId,
  });
}

export function transformChargeCurrent(input: {
  rho: number;
  J: Vec3;
  boost: number;
  c?: number;
}): Readonly<{
  rho: number;
  J: Vec3;
  beta: number;
  gamma: number;
}> {
  const c = input.c ?? C_SI;
  const v = input.boost;
  const beta = v / c;
  const g = gamma(beta);
  if (g.status !== "value") {
    return Object.freeze({
      rho: Number.NaN,
      J: Object.freeze({ x: Number.NaN, y: Number.NaN, z: Number.NaN }),
      beta,
      gamma: Number.NaN,
    });
  }
  const γ = g.value;
  const { rho, J } = input;
  const rhoPrime = γ * (rho - (v * J.x) / (c * c));
  const JprimeX = γ * (J.x - v * rho);
  return Object.freeze({
    rho: rhoPrime,
    J: Object.freeze({
      x: JprimeX,
      y: J.y,
      z: J.z,
    }),
    beta,
    gamma: γ,
  });
}