springPhysics

fn springPhysics(t: f32, targetPosition: f32, initialPos: f32, initialVel: f32, stiffness: f32, damping: f32, mass: f32) -> vec2<f32> {
  // Ensure positive values for stiffness, damping, and mass
  let k = max(0.0001, stiffness);
  let d = max(0.0, damping);
  let m = max(0.0001, mass);
  
  // Calculate the angular frequency and damping ratio
  let omega = sqrt(k / m);
  let zeta = d / (2.0 * sqrt(k * m));
  
  // Initial displacement from targetPosition position
  let x0 = initialPos - targetPosition;
  let v0 = initialVel;
  
  var position: f32 = 0.0;
  var velocity: f32 = 0.0;
  
  if (zeta < 1.0) {
    // Underdamped case
    let omega_d = omega * sqrt(1.0 - zeta * zeta);
    let A = x0;
    let B = (v0 + zeta * omega * x0) / omega_d;
    
    // Calculate exponential decay term
    let expTerm = exp(-zeta * omega * t);
    
    // Calculate position and velocity
    position = targetPosition + expTerm * (A * cos(omega_d * t) + B * sin(omega_d * t));
    velocity = expTerm * (
                -zeta * omega * A * cos(omega_d * t) - omega_d * A * sin(omega_d * t) +
                -zeta * omega * B * sin(omega_d * t) + omega_d * B * cos(omega_d * t)
               );
  } else if (zeta == 1.0) {
    // Critically damped case
    let A = x0;
    let B = v0 + omega * x0;
    
    // Calculate exponential decay term
    let expTerm = exp(-omega * t);
    
    // Calculate position and velocity
    position = targetPosition + expTerm * (A + B * t);
    velocity = expTerm * (B - omega * (A + B * t));
  } else {
    // Overdamped case
    let omega1 = -omega * (zeta + sqrt(zeta * zeta - 1.0));
    let omega2 = -omega * (zeta - sqrt(zeta * zeta - 1.0));
    
    let A = (v0 - omega2 * x0) / (omega1 - omega2);
    let B = x0 - A;
    
    // Calculate position and velocity
    position = targetPosition + A * exp(omega1 * t) + B * exp(omega2 * t);
    velocity = A * omega1 * exp(omega1 * t) + B * omega2 * exp(omega2 * t);
  }
  
  return vec2<f32>(position, velocity);
}