/*
 * PianoPhysicsData.h
 * 
 * Implementation-ready parameter values for piano physical modeling
 * Extracted from peer-reviewed acoustic research papers
 * 
 * Sources:
 *   - Chaigne et al., JASA 133(4), 2013 - Soundboard dynamics
 *   - Hall & Askenfelt, JASA 83(4), 1988 - Hammer-string interaction
 *   - Russell & Rossing, Acta Acustica 84, 1998 - Hammer nonlinearity
 *   - Giordano, JASA 103(4), 1998 - Bridge impedance
 *   - Conklin, JASA 99-100, 1996 - Piano design
 *   - Fletcher, JASA 36(1), 1964 - String inharmonicity
 * 
 * This data is scientific fact and not subject to copyright.
 */

#pragma once
#include <cmath>
#include <array>

namespace PianoPhysics {

//==============================================================================
// MATERIAL PROPERTIES
//==============================================================================

// European Spruce (Picea excelsa) - Soundboard
namespace Spruce {
    constexpr float density = 440.0f;           // kg/m³
    constexpr float E_longitudinal = 15.9e9f;   // Pa (along grain)
    constexpr float E_radial = 0.69e9f;         // Pa
    constexpr float E_tangential = 0.39e9f;     // Pa
    constexpr float G_LR = 0.62e9f;             // Pa (shear modulus)
    constexpr float G_LT = 0.77e9f;             // Pa
    constexpr float G_RT = 0.0036e9f;           // Pa
    constexpr float nu_LR = 0.44f;              // Poisson's ratio
    constexpr float nu_LT = 0.38f;
    constexpr float nu_RT = 0.47f;
}

// Beech - Bridge material
namespace Beech {
    constexpr float density = 674.0f;           // kg/m³
    constexpr float E_longitudinal = 14.0e9f;   // Pa
    constexpr float E_radial = 2.28e9f;         // Pa
    constexpr float E_tangential = 1.16e9f;     // Pa
}

// Steel - Strings
namespace Steel {
    constexpr float density = 7800.0f;          // kg/m³
    constexpr float youngsModulus = 210.0e9f;   // Pa
}

// Copper - Bass string wrapping
namespace Copper {
    constexpr float density = 8960.0f;          // kg/m³
}

//==============================================================================
// SOUNDBOARD PARAMETERS
//==============================================================================

namespace Soundboard {
    // Modal damping (loss factor eta)
    constexpr float lossFactor_low = 0.007f;    // < 300 Hz (0.7%)
    constexpr float lossFactor_mid = 0.012f;    // 300-1000 Hz (1.2%)
    constexpr float lossFactor_high = 0.020f;   // 1000-3000 Hz (2.0%)
    
    // Typical loss factor for good overall match
    constexpr float typicalLossFactor = 0.015f; // 1.5%
    
    // Waveguide transition frequency (ribs become waveguides above this)
    constexpr float waveguideTransitionHz = 1100.0f;
    
    // Mean rib spacing (Pleyel P131 measurements)
    constexpr float meanRibSpacing_m = 0.132f;  // 13.2 cm
    
    // T60 values by frequency band (seconds)
    inline float getT60(float freqHz) {
        if (freqHz < 120.0f)  return 1.8f - (freqHz - 60.0f) * 0.012f;
        if (freqHz < 500.0f)  return 1.1f - (freqHz - 120.0f) * 0.002f;
        if (freqHz < 2000.0f) return 0.7f - (freqHz - 500.0f) * 0.0003f;
        return 0.25f - (freqHz - 2000.0f) * 0.00003f;
    }
}

//==============================================================================
// STRING PARAMETERS (Broadwood Grand measurements)
//==============================================================================

struct StringData {
    int midiNote;
    float frequency;     // Hz (equal temperament)
    float length;        // m
    float coreDiameter;  // mm
    float wrapDiameter;  // mm (0 = unwrapped)
    float mass;          // g (total vibrating mass)
    float tension;       // N
    float strikingRatio; // hammer position / string length
};

// Measured string parameters for C notes
constexpr std::array<StringData, 7> stringDataC = {{
    // MIDI, freq,   length, core, wrap,  mass,  tension, strike
    {  24,   32.7f,  1.945f, 1.00f, 4.55f, 142.0f, 680.0f, 0.117f }, // C1
    {  36,   65.4f,  1.425f, 1.00f, 3.60f,  52.0f, 710.0f, 0.102f }, // C2
    {  48,  130.8f,  0.975f, 1.025f,1.75f,  21.0f, 740.0f, 0.109f }, // C3
    {  60,  261.6f,  0.620f, 1.025f,0.0f,   15.0f, 770.0f, 0.125f }, // C4
    {  72,  523.3f,  0.337f, 1.00f, 0.0f,    8.1f, 810.0f, 0.134f }, // C5
    {  84, 1046.5f,  0.172f, 0.925f,0.0f,    3.9f, 680.0f, 0.151f }, // C6
    {  96, 2093.0f,  0.089f, 0.850f,0.0f,    1.7f, 560.0f, 0.180f }, // C7
}};

//==============================================================================
// INHARMONICITY
//==============================================================================

namespace Inharmonicity {
    // Measured B values at specific points
    // f_n = n * f1 * sqrt(1 + B * n²)
    constexpr float B_A0  = 0.00020f;   // MIDI 21
    constexpr float B_C2  = 0.00030f;   // MIDI 36 (wound)
    constexpr float B_C3  = 0.00070f;   // MIDI 48
    constexpr float B_C4  = 0.0015f;    // MIDI 60
    constexpr float B_C5  = 0.0050f;    // MIDI 72
    constexpr float B_C6  = 0.018f;     // MIDI 84
    constexpr float B_C7  = 0.080f;     // MIDI 96
    constexpr float B_C8  = 0.40f;      // MIDI 108
    
    // Interpolation function for any MIDI note
    inline float getB(int midiNote) {
        // Bass wound strings (MIDI 21-52): slow exponential growth
        if (midiNote <= 52) {
            float norm = static_cast<float>(midiNote - 21) / 31.0f;
            return 0.00018f * std::pow(10.0f, norm * 0.7f);
        }
        // Plain steel strings (MIDI 53-108): faster growth
        float norm = static_cast<float>(midiNote - 53) / 55.0f;
        return 0.0008f * std::pow(10.0f, norm * 2.7f);
    }
    
    // Compute partial frequency with inharmonicity
    inline float partialFreq(float f1, int n, float B) {
        return static_cast<float>(n) * f1 * std::sqrt(1.0f + B * n * n);
    }
}

//==============================================================================
// HAMMER PARAMETERS
//==============================================================================

namespace Hammer {
    // Nonlinear felt model: F = K * C^p
    
    // Measured values for specific notes (Chaigne & Askenfelt)
    struct HammerData {
        int midiNote;
        float K;        // Stiffness (SI units, varies with p)
        float p;        // Nonlinearity exponent
        float mass_kg;  // Mass in kg
    };
    
    constexpr std::array<HammerData, 3> measuredHammers = {{
        { 36, 4.0e8f, 2.3f, 0.0110f },   // C2
        { 60, 4.5e9f, 2.8f, 0.0092f },   // C4  
        { 84, 1.0e10f, 3.0f, 0.0052f },  // C6
    }};
    
    // Exponent p varies smoothly from bass to treble
    // Source: Russell & Rossing (1998)
    inline float getExponent(int midiNote) {
        // Bass: p ≈ 2.0-2.3
        // Middle: p ≈ 2.6-3.0  
        // Treble: p ≈ 3.0-4.0
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        return 2.0f + norm * 2.0f;  // 2.0 to 4.0
    }
    
    // Stiffness coefficient K (log-interpolated)
    inline float getStiffness(int midiNote) {
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        return 4.0e8f * std::pow(10.0f, norm * 1.4f);
    }
    
    // Hammer mass (kg)
    inline float getMass(int midiNote) {
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        return 0.012f - norm * 0.008f;  // 12g to 4g
    }
    
    // Hysteresis (reduces force on return stroke)
    // Source: Stulov
    inline float getHysteresis(int midiNote) {
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        return 0.18f - norm * 0.10f;  // 0.18 to 0.08
    }
    
    // Striking position as fraction of string length
    inline float getStrikingRatio(int midiNote) {
        // Varies from ~1/8.5 in bass to ~1/5.5 in treble
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        return 0.10f + norm * 0.08f;  // 0.10 to 0.18
    }
    
    // Contact duration scaling with velocity
    // τ ∝ F_max^((1-p)/(2p))
    // At higher velocity, contact is shorter
    inline float getContactDurationScale(float velocity01, float p) {
        // velocity01: 0 = pp, 1 = ff
        // Returns multiplier relative to pp duration
        float Fscale = 0.2f + 0.8f * velocity01;  // Force scales with velocity
        float exponent = (1.0f - p) / (2.0f * p);
        return std::pow(Fscale, exponent);
    }
}

//==============================================================================
// BRIDGE IMPEDANCE
//==============================================================================

namespace Bridge {
    // Mean impedance values (kg/s)
    // Source: Wogram, Conklin, Giordano
    
    constexpr float Z_low = 1000.0f;     // < 1 kHz
    constexpr float Z_mid = 700.0f;      // 1-2 kHz
    constexpr float Z_high = 400.0f;     // 2-5 kHz
    constexpr float Z_veryHigh = 250.0f; // > 5 kHz
    
    // Typical fluctuation around mean: ±12 dB
    
    // High-pass cutoff frequencies
    constexpr float bassBridgeHP_Hz = 70.0f;
    constexpr float trebleBridgeHP_Hz = 140.0f;
    
    // Bridge/soundboard coupling strength
    constexpr float bassCoupling = 0.06f;
    constexpr float trebleCoupling = 0.04f;
    
    // Bass/treble split point
    constexpr int splitMidi = 52;  // E3
    
    inline float getImpedance(float freqHz) {
        if (freqHz < 1000.0f) return Z_low;
        if (freqHz < 2000.0f) return Z_low - (freqHz - 1000.0f) * 0.0003f;
        if (freqHz < 5000.0f) return Z_mid - (freqHz - 2000.0f) * 0.0001f;
        return Z_veryHigh;
    }
}

//==============================================================================
// DAMPER PARAMETERS
//==============================================================================

namespace Damper {
    // Damping time constant when damper engages (approximate)
    constexpr float engageTime_s = 0.015f;    // 15 ms
    
    // Final damping coefficient (very high = quick stop)
    constexpr float maxDamping = 0.90f;
    
    // Bass strings have more mass, damp slower
    inline float getDampingRate(int midiNote) {
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        return 0.85f + norm * 0.10f;  // 0.85 to 0.95
    }
}

//==============================================================================
// OVERALL STRING DECAY (T60)
//==============================================================================

namespace StringDecay {
    // Frequency-dependent T60 for string modes
    // Higher partials decay faster
    
    inline float getT60(float freqHz, int midiNote) {
        // Base T60 varies with register
        float norm = static_cast<float>(midiNote - 21) / 87.0f;
        float baseT60 = 20.0f - norm * 16.5f;  // 20s bass to 3.5s treble
        
        // High frequencies decay faster (loss filter effect)
        float freqFactor = 1.0f / (1.0f + freqHz / 5000.0f);
        
        return baseT60 * freqFactor;
    }
}

//==============================================================================
// VELOCITY MAPPING
//==============================================================================

namespace Velocity {
    // MIDI velocity to hammer velocity (m/s)
    // Typical range: 0.5 m/s (pp) to 5.5 m/s (fff)
    
    inline float midiToHammerVelocity(int midiVel) {
        float norm = static_cast<float>(midiVel) / 127.0f;
        // Logarithmic-ish mapping for musical response
        return 0.5f + 5.0f * norm * norm;
    }
    
    // Brightness tilt with velocity
    // Louder = brighter (more high frequency content)
    inline float getBrightnessTilt(float hammerVel) {
        // 0 to 1 range
        return std::min(1.0f, hammerVel / 5.0f);
    }
}

} // namespace PianoPhysics