// The felt of the stand as a viscoelastic material: a spring in parallel
// with a chain of Maxwell arms (a spring and a dashpot in series), whose
// relaxation times are spaced two per decade from 0.5 Hz to above the audio
// band and whose weights grow geometrically. Such a chain has a loss factor
// that is nearly flat over the band, which is what a single loss factor
// eta means, and unlike a dashpot it is causal and can be stepped in time:
// each arm carries one state per node, its extension p, which relaxes
// toward zero at its own rate while the node moves.
//
// Every quantity here is normalized to the stiffness at the reference
// frequency (the bounce of the plate on the stand), where the storage
// stiffness of the chain is 1, so that the felt modulus given by the user is
// the modulus at the frequency where it matters most. Below the reference
// the felt is softer and above it stiffer, as the Kramers-Kronig relations
// require of any material with a loss; at rest its stiffness is kInf.

export interface FeltChain {
  /** Relaxation time of each arm [s]. */
  tau: Float64Array;
  /** Stiffness of each arm, per unit of the stiffness at the reference frequency. */
  k: Float64Array;
  /** Relaxed (static) stiffness, per unit of the stiffness at the reference frequency. */
  kInf: number;
  /** The loss factor the chain was built for. */
  eta: number;
}

/** Lowest arm frequency [Hz]; the highest is FELT_ARMS_PER_DECADE arms per decade above it. */
export const FELT_LOWEST_HZ = 0.5;
export const FELT_ARMS_PER_DECADE = 2;
/** Number of arms: 0.5 Hz to 158 kHz. */
export const FELT_ARMS = 12;

/**
 * The chain for a loss factor eta, normalized at the frequency fRef [Hz]. A
 * chain of arms with frequencies in the ratio r and stiffnesses in the ratio
 * 1 + c approximates the power law K(i omega) ~ (i omega)^(2 delta / pi) of a
 * material with the constant loss angle delta = atan(eta), whose modulus
 * rises by the factor r^(2 delta / pi) per arm; c is that factor minus one.
 */
export function feltChain(eta: number, fRef: number): FeltChain {
  const r = Math.pow(10, 1 / FELT_ARMS_PER_DECADE);
  const c = eta > 0 ? Math.pow(r, (2 * Math.atan(eta)) / Math.PI) - 1 : 0;
  const nArms = c > 0 ? FELT_ARMS : 0;
  const tau = new Float64Array(nArms), k = new Float64Array(nArms);
  let S = 1;
  for (let i = 0; i < nArms; i++) {
    tau[i] = 1 / (2 * Math.PI * FELT_LOWEST_HZ * Math.pow(r, i));
    k[i] = c * S;
    S *= 1 + c;
  }
  const raw: FeltChain = { tau, k, kInf: 1, eta };
  const scale = 1 / feltModulus(raw, 2 * Math.PI * Math.max(fRef, 1e-6)).storage;
  for (let i = 0; i < nArms; i++) k[i] *= scale;
  return { tau, k, kInf: scale, eta };
}

/** Storage and loss stiffness of the chain at the angular frequency omega, per unit of the reference stiffness. */
export function feltModulus(ch: FeltChain, omega: number): { storage: number; loss: number } {
  let storage = ch.kInf, loss = 0;
  for (let i = 0; i < ch.tau.length; i++) {
    const x = omega * ch.tau[i], d = 1 + x * x;
    storage += (ch.k[i] * x * x) / d;
    loss += (ch.k[i] * x) / d;
  }
  return { storage, loss };
}

/**
 * The constants of the exact time step of the arms for a step dt: with the
 * node moving linearly by dw over the step, the extension of arm i goes from
 * p to p decay_i + dw growth_i, where decay_i = exp(-dt / tau_i) and
 * growth_i = (tau_i / dt)(1 - decay_i). The force of the arms on the node is
 * -k_node sum_i k_i p_i.
 */
export function feltStepConstants(ch: FeltChain, dt: number): { decay: Float64Array; growth: Float64Array } {
  const n = ch.tau.length;
  const decay = new Float64Array(n), growth = new Float64Array(n);
  for (let i = 0; i < n; i++) {
    const x = dt / ch.tau[i];
    decay[i] = Math.exp(-x);
    growth[i] = -Math.expm1(-x) / x;
  }
  return { decay, growth };
}