junction.hpp

Go to the documentation of this file.

#pragma once

 
#include <algorithm>
#include <cmath>
 
namespace pulp::signal::junction {

inline constexpr double kThermalVoltage = 0.02585;

inline constexpr double kMaxExponent = 80.0;

inline double safe_exp(double x) {
    return std::exp(std::clamp(x, -kMaxExponent, kMaxExponent));
}

inline double exprel(double x) {
    if (std::abs(x) < 1e-6) return 1.0 + x * 0.5 + x * x * (1.0 / 6.0);
    return std::expm1(std::clamp(x, -kMaxExponent, kMaxExponent)) / x;
}

struct JunctionPair {
    double saturation_current = 1e-14;

    double ideality = 1.0;

    double leg_a = 1.0;

    double leg_b = 1.0;

    double theta() const { return ideality * kThermalVoltage; }

    double current(double v) const {
        const double t = theta();
        return saturation_current * (safe_exp(leg_a * v / t) - safe_exp(-leg_b * v / t));
    }

    double conductance(double v) const {
        const double t = theta();
        return saturation_current / t *
               (leg_a * safe_exp(leg_a * v / t) + leg_b * safe_exp(-leg_b * v / t));
    }

    double antiderivative(double v) const {
        const double t = theta();
        return saturation_current * v * (exprel(leg_a * v / t) - exprel(-leg_b * v / t));
    }

    double conduction_estimate(double v_drive, double resistance) const {
        const double leg = v_drive >= 0.0 ? leg_a : leg_b;
        if (!(leg > 0.0)) return v_drive;
        const double ratio = std::abs(v_drive) / (resistance * saturation_current);
        if (!(ratio > 1.0)) return v_drive;
        const double magnitude = (theta() / leg) * std::log(ratio);
        return std::copysign(std::min(magnitude, std::abs(v_drive)), v_drive);
    }

    double knee_voltage(double target_current) const {
        if (!(leg_a > 0.0)) return 0.0;
        return (theta() / leg_a) * std::log(target_current / saturation_current + 1.0);
    }

    double adaa_current(double v, double previous, double epsilon) const {
        const double delta = v - previous;
        if (std::abs(delta) <= epsilon) return current(v);
        return (antiderivative(v) - antiderivative(previous)) / delta;
    }

    double adaa_conductance(double v, double previous, double epsilon) const {
        const double delta = v - previous;
        // HALF the conductance, not the conductance. As `v → p` the quotient
        // `(F1(v) − F1(p))/(v − p)` tends to `i(p)` but its SLOPE tends to
        // `i'(p)/2` — the mean of `i` over a collapsing interval moves at half
        // the rate of `i` itself. Returning the full conductance here overstates
        // the Jacobian by exactly 2x near the fallback, which is where a solve
        // spends most of its time on a slow-moving signal.
        if (std::abs(delta) <= epsilon) return 0.5 * conductance(v);
        const double integral = antiderivative(v) - antiderivative(previous);
        return (current(v) * delta - integral) / (delta * delta);
    }
};
 
}  // namespace pulp::signal::junction