高槻・北摂の街に響く音楽文化メセナ：
室内楽ホールの音響物理学（初期反射音・残響時間・IACC）と地域の音楽芸術団体40年の歩み

#![forbid(unsafe_code)]

/// Computes Norris-Eyring reverberation time RT60 (seconds) with air absorption 4mV.
#[inline]
pub fn compute_eyring_rt60(
    volume_m3: f64,
    surface_area_m2: f64,
    mean_absorption_alpha: f64,
    air_attenuation_m: f64,
) -> f64 {
    let alpha = mean_absorption_alpha.clamp(0.01, 0.99);
    let denom = -surface_area_m2 * (1.0 - alpha).ln() + 4.0 * air_attenuation_m * volume_m3;
    if denom <= 1e-9 {
        return 0.0;
    }
    (0.161 * volume_m3) / denom
}

/// Computes ISO 3382-1 Clarity Index C80 (in dB) from a discrete Room Impulse Response (RIR).
pub fn compute_c80_db(impulse_response: &[f64], sample_rate_hz: usize) -> f64 {
    let split_idx = ((sample_rate_hz as f64) * 0.080).round() as usize;
    let split = split_idx.min(impulse_response.len());

    let mut early_energy = 0.0f64;
    for &x in &impulse_response[..split] {
        early_energy += x * x;
    }

    let mut late_energy = 0.0f64;
    for &x in &impulse_response[split..] {
        late_energy += x * x;
    }

    if late_energy <= 1e-18 {
        return 60.0;
    }
    10.0 * (early_energy / late_energy).max(1e-12).log10()
}

/// Computes the Interaural Cross-Correlation Coefficient (IACC) over |tau| <= 1 ms.
pub fn compute_binaural_iacc(
    left_ir: &[f64],
    right_ir: &[f64],
    sample_rate_hz: usize,
) -> f64 {
    let n = left_ir.len().min(right_ir.len());
    if n == 0 {
        return 1.0;
    }

    let mut e_left = 0.0f64;
    let mut e_right = 0.0f64;
    for i in 0..n {
        e_left += left_ir[i] * left_ir[i];
        e_right += right_ir[i] * right_ir[i];
    }
    let norm = (e_left * e_right).sqrt();
    if norm <= 1e-18 {
        return 0.0;
    }

    let max_lag = ((sample_rate_hz as f64) * 0.001).round() as isize;
    let mut max_corr = 0.0f64;

    for lag in -max_lag..=max_lag {
        let mut cross = 0.0f64;
        for i in 0..n {
            let j = (i as isize) + lag;
            if j >= 0 && (j as usize) < n {
                cross += left_ir[i] * right_ir[j as usize];
            }
        }
        let c = (cross / norm).abs();
        if c > max_corr {
            max_corr = c;
        }
    }

    max_corr.clamp(0.0, 1.0)
}