thermal_friction.hpp

Go to the documentation of this file.

// sbl/dsp/pm/thermal_friction.hpp — A bow whose friction remembers: the thermal model (Physical modeling — primitive)
//
// The friction-curve bow (table_friction.hpp) makes force a function of
// sliding speed alone, which rosin measurably does not do: during stick-slip
// it traces a hysteresis loop in the force–velocity plane that follows no
// steady-sliding curve (Smith & Woodhouse 2000). Woodhouse (2003) proposed
// the simplest law that does: the coefficient of friction depends on the
// temperature of the rosin at the contact, and on nothing else. Sliding
// heats the contact and the friction falls; sticking lets it cool and the
// friction recovers. The lag between speed and temperature is the loop, and
// with it come the results that matter for a player: Helmholtz motion is
// reached more reliably and more quickly, the Schelleng region is wider,
// and the string is less "twitchy" (Woodhouse 2003 §5). Maestre, Spa and
// Smith (2014) put the same idea in a digital waveguide.
//
// The model, in the waveguide's units (2Z = 1, as table_friction.hpp):
//
//     limit(T) = L(force) · μ(T)/μ_cold       the static limit, cold at L
//     stuck   : |dv_free| <= limit(T)  →  F = dv_free
//     sliding : F = limit(T) · sgn(dv_free)     Coulomb: no speed dependence
//     dT/dt   = heat · F · v_slip  −  T/τ  −  convect · v_slip · T
//
// with v_slip = |dv_free| − |F| the true sliding speed. Friction that does
// not depend on speed has a flat sliding branch, so the load line crosses it
// in exactly one place: no capture test, no rounding of the corner
// (Woodhouse 2003 §2.3). Release and capture happen at different limits
// because T differs — that is the hysteresis, and it comes from the state,
// not from a rule.
//
// T is kept as a fraction of the range over which μ falls, 0 cold to 1 fully
// softened, so the constants read as "how fast it heats" and "how fast it
// cools" in that range: Woodhouse's rosin falls from μ 1.2 to 0.35 over about
// 40 °C above ambient (his Figure 3), and a played string sits at 17–31 °C
// (his Figure 6c). Exact values matter less than they look: the system
// "self-buffers", warming until the falling friction stops the warming
// (§2.1). Calibrate with pictures, not by hand.
//
// Conduction into a solid is not one time constant: the Green's function
// falls as 1/√t, fast at first and slow after, which is what lets a contact
// cool within a two-millisecond stick at C5 and still keep a warm baseline
// across the period. Two exponentials stand in for it here — a fast one that
// takes `fast_share` of the heat and a slow one that takes the rest — and
// T is their sum.
//
// The same state can also drive a friction-curve bow (BowedString's Hybrid
// model): the table's speed-dependent force scaled by the rosin's warmth —
// f = N·μ(T) with the curve's shape kept, the combination Woodhouse's
// conclusions point to (2003 §5). advance() feeds the state from a force and
// sliding speed the curve found.
//
// L(force) is the friction-curve bow's per-slice static limit, scaled by
// `cold_gain`: with the default 1/μ_ratio the *hot* limit equals the table's
// limit, so a bow pressure means "the friction of warm rosin equals the table
// bow's stick limit", and the two bows can be drawn side by side on one
// Schelleng diagram. (Calibrating the cold limit to the table's instead
// starves the string: Coulomb friction has no falling slope to self-excite,
// so the corner alone must recapture the string, and it needs the hot
// friction to be about what the table bow sticks with. Woodhouse's cello
// string sits at a cold limit near twenty bow speeds and a hot one near six.)
//
// Energy: a bow is a source, and this one has no passivity proof either.
// |F| never exceeds |dv_free| — sliding force is at most the limit, and
// sliding means the limit is below |dv_free| — so it cannot drive the string
// faster than the bow; the SoftLimiter downstream stays load-bearing.
//
// A law, not an exciter: a Bow<ThermalFriction> owns the speed, the weight
// and the junction. process() is the law solved against a string's line at
// unit impedance; limit() and advance() are what a junction needs to solve
// it against any other line (AP-041 Phase 6, the rubbed bowl).
//
// Usage:
//   sbl::dsp::pm::ThermalFriction law;
//   law.set_limits(lut::bow_stk_limits, 8);
//   law.set_force_su(pressure);
//   float injection = law.process(v_bow - v_incoming);
 
#ifndef SBL_DSP_PM_THERMAL_FRICTION_HPP_
#define SBL_DSP_PM_THERMAL_FRICTION_HPP_
 
#include <cmath>
#include <cstdint>
 
#include <sbl/dsp/diagram/graph.hpp>
#include <sbl/dsp/math/clamp.hpp>
#include <sbl/dsp/pm/contracts.hpp>
#include <sbl/dsp/types/fixed.hpp>
 
namespace sbl::dsp::pm {
 
class ThermalFriction {
public:
    /// @note All public methods are ISR-safe — bounded computation, no I/O.
 
    /// Rosin, as Woodhouse measured it: friction at full heat is this fraction of cold.
    static constexpr float DEFAULT_MU_RATIO = 0.35f / 1.2f;
    /// Fraction of the range per unit of (force · sliding speed) per second.
    /// Calibrated 2026-09-15 on the Davis Jr. string: Helmholtz at C3 with the body,
    /// C3–C5 bare (AP-040 Phase 4 addenda).
    static constexpr float DEFAULT_HEAT = 2000.0f;
    /// Conduction into the string and the bow: the slow component's time to cool by 1/e.
    static constexpr float DEFAULT_TAU_MS = 10.0f;
    /// Convection: cold rosin flowing into the contact, per unit sliding speed per second.
    static constexpr float DEFAULT_CONVECT = 0.0f;
    /// The cold limit over the table's. 1: cold equals the table's limit. (Any
    /// more and this waveguide's stuck force, capped by loop loss at a few bow
    /// speeds, never reaches it.)
    static constexpr float DEFAULT_COLD_GAIN = 1.0f;
    /// The fast cooling component: time constant and the share of heat it takes.
    static constexpr float DEFAULT_FAST_MS = 0.5f;
    static constexpr float DEFAULT_FAST_SHARE = 0.6f;
 
    /**
     * @brief The cold static limit per force slice — the friction-curve bow's `limits`
     *
     * @param limits  one static-friction limit per slice, light bow first
     * @param slices  number of slices
     */
    void set_limits(const float* limits, uint8_t slices) {
        limits_ = limits;
        slices_ = slices;
        set_force_su(force_su_);
    }
 
    /// Bow force across the slices [0, 1]; the cold limit interpolates between them.
    void set_force_su(float su) {
        force_su_ = su < 0.0f ? 0.0f : (su > 1.0f ? 1.0f : su);
        if (limits_ == nullptr || slices_ == 0) {
            cold_limit_ = 0.0f;
            return;
        }
        const float pos = force_su_ * static_cast<float>(slices_ - 1);
        uint8_t i = static_cast<uint8_t>(pos);
        if (i >= slices_ - 1 && slices_ > 1) i = static_cast<uint8_t>(slices_ - 2);
        const float blend = slices_ > 1 ? pos - static_cast<float>(i) : 0.0f;
        cold_limit_ = cold_gain_ * (slices_ > 1 ? limits_[i] + blend * (limits_[i + 1] - limits_[i]) : limits_[0]);
    }
 
    /**
     * @brief The rosin: how it heats, cools and softens (configuration, not signals)
     *
     * @param heat     temperature fraction per unit of force·speed per second
     * @param tau_ms   conduction time constant
     * @param convect  convection coefficient per unit sliding speed per second
     * @param mu_ratio friction at full heat relative to cold, in (0, 1]
     * @param cold_gain the cold limit as a multiple of the table's limit
     */
    void set_rosin(float heat, float tau_ms, float convect, float mu_ratio,
                   float cold_gain = DEFAULT_COLD_GAIN) {
        heat_ = heat > 0.0f ? heat : 0.0f;
        cool_ = tau_ms > 0.0f ? 1.0f / (tau_ms * 0.001f) : 0.0f;
        convect_ = convect > 0.0f ? convect : 0.0f;
        mu_ratio_ = (mu_ratio > 0.0f && mu_ratio <= 1.0f) ? mu_ratio : DEFAULT_MU_RATIO;
        cold_gain_ = cold_gain > 0.0f ? cold_gain : DEFAULT_COLD_GAIN;
        set_force_su(force_su_);
    }
 
    /// The fast cooling component: its time constant and the share of the heat it takes.
    void set_fast_cooling(float fast_ms, float fast_share) {
        cool_fast_ = fast_ms > 0.0f ? 1.0f / (fast_ms * 0.001f) : 0.0f;
        fast_share_ = fast_share < 0.0f ? 0.0f : (fast_share > 1.0f ? 1.0f : fast_share);
    }
 
    /// Friction relative to cold, from the temperature: 1 cold down to mu_ratio hot.
    float mu_rel() const { return 1.0f - (1.0f - mu_ratio_) * temperature(); }
 
    /**
     * @brief Advance the contact one sample from a force and sliding speed found elsewhere
     *
     * For the hybrid bow: the friction curve finds the force; this keeps the
     * temperature. Sliding speed 0 while stuck.
     */
    void advance(float force, float v_slip) {
        const float f = std::fabs(force);
        const float v = v_slip < 0.0f ? 0.0f : v_slip;
        const float heating = heat_ * f * v;
        t_fast_ += dt_ * (fast_share_ * heating - (cool_fast_ + convect_ * v) * t_fast_);
        t_slow_ += dt_ * ((1.0f - fast_share_) * heating - (cool_ + convect_ * v) * t_slow_);
        t_fast_ = math::clamp01(t_fast_);
        t_slow_ = math::clamp01(t_slow_);
    }
 
    /**
     * @brief One sample: relative velocity in, injection force out
     *
     * @param dv_free v_bow - v_incoming, in wave variables
     */
    float process(float dv_free) {
        if (limits_ == nullptr) return 0.0f;
        const float a = std::fabs(dv_free);
        const float limit = cold_limit_ * mu_rel();
 
        float f;
        float v_slip;
        if (a <= limit || !(a == a)) {   // holds, or NaN: treat as at rest
            stuck_ = true;
            f = (a == a) ? dv_free : 0.0f;
            v_slip = 0.0f;
        } else {
            stuck_ = false;
            f = dv_free < 0.0f ? -limit : limit;
            v_slip = a - limit;
        }
 
        // The contact: heated by the work of sliding, cooled by conduction
        // always and by convection while sliding. Forward Euler; at audio
        // rate the step is a fortieth of the shortest time constant.
        advance(limit, v_slip);
        return f;
    }
 
    /// True while the string travels with the bow.
    bool sticking() const { return stuck_; }
 
    /// The static limit right now, warmed: where the string would release.
    float limit() const { return cold_limit_ * mu_rel(); }
 
    /// Contact temperature as a fraction of the range over which the rosin softens.
    float temperature() const {
        const float t = t_fast_ + t_slow_;
        return t > 1.0f ? 1.0f : t;
    }
 
    /// Cold and stuck.
    void reset() {
        stuck_ = true;
        t_fast_ = 0.0f;
        t_slow_ = 0.0f;
        dt_ = 1.0f / types::SAMPLE_RATE_F;   // the one sample-rate rule (FDP-055 Addendum A)
    }
 
    /// One Junction node: the friction nonlinearity. Never called from audio code (AP-037).
    diagram::Ports describe(diagram::Graph& g, uint8_t group, const char* name = "bow (thermal)") const {
        diagram::Ports p;
        p.in = p.out = g.add_node(diagram::Kind::Junction, name, group);
        p.group = group;
        return p;
    }
 
private:
    const float* limits_ = nullptr;
    uint8_t slices_ = 0;
    float force_su_ = 1.0f;
    float cold_limit_ = 0.0f;
 
    float heat_ = DEFAULT_HEAT;
    float cool_ = 1.0f / (DEFAULT_TAU_MS * 0.001f);
    float convect_ = DEFAULT_CONVECT;
    float mu_ratio_ = DEFAULT_MU_RATIO;
    float cold_gain_ = DEFAULT_COLD_GAIN;
    float cool_fast_ = 1.0f / (DEFAULT_FAST_MS * 0.001f);
    float fast_share_ = DEFAULT_FAST_SHARE;
    float dt_ = 1.0f / types::SAMPLE_RATE_F;
 
    float t_fast_ = 0.0f;
    float t_slow_ = 0.0f;
    bool stuck_ = true;
};
 
static_assert(is_friction_law_v<ThermalFriction>, "ThermalFriction is a friction law (physical-modeling.md §4.5)");
 
} // namespace sbl::dsp::pm
 
#endif // SBL_DSP_PM_THERMAL_FRICTION_HPP_