Integrator6581.h

/*
 * This file is part of libsidplayfp, a SID player engine.
 *
 * Copyright 2011-2021 Leandro Nini <drfiemost@users.sourceforge.net>
 * Copyright 2007-2010 Antti Lankila
 * Copyright 2004, 2010 Dag Lem <resid@nimrod.no>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 */
 
#ifndef INTEGRATOR6581_H
#define INTEGRATOR6581_H
 
#include <stdint.h>
#include <cassert>
#ifdef SLOPE_FACTOR
#  include <cmath>
#endif
 
#include "siddefs-fp.h"
 
namespace reSIDfp
{
 
class Integrator6581
{
private:
    const unsigned short* vcr_Vg;
    const unsigned short* vcr_n_Ids_term;
    const unsigned short* opamp_rev;
 
    unsigned int Vddt_Vw_2;
    mutable int vx;
    mutable int vc;
#ifdef SLOPE_FACTOR
    // Slope factor n = 1/k
    // where k is the gate coupling coefficient
    // k = Cox/(Cox+Cdep) ~ 0.7
    mutable double n;
#else
    const int n;
#endif
    const double N16;
    const unsigned short Vddt;
    const unsigned short nVt;
    const unsigned short nVmin;
    const unsigned short n_snake;
 
public:
    Integrator6581(const unsigned short* vcr_Vg, const unsigned short* vcr_n_Ids_term,
               const unsigned short* opamp_rev, unsigned short Vddt, unsigned short nVt,
               unsigned short nVmin, unsigned short n_snake, double N16) :
        vcr_Vg(vcr_Vg),
        vcr_n_Ids_term(vcr_n_Ids_term),
        opamp_rev(opamp_rev),
        Vddt_Vw_2(0),
        vx(0),
        vc(0),
#ifdef SLOPE_FACTOR
        n(1.4),
#else
        n(1),
#endif
        N16(N16),
        Vddt(Vddt),
        nVt(nVt),
        nVmin(nVmin),
        n_snake(n_snake) {}
 
    void setVw(unsigned short Vw) { Vddt_Vw_2 = ((Vddt - Vw) * (Vddt - Vw)) >> 1; }
 
    int solve(int vi) const;
};
 
} // namespace reSIDfp
 
#if RESID_INLINING || defined(INTEGRATOR_CPP)
 
namespace reSIDfp
{
 
RESID_INLINE
int Integrator6581::solve(int vi) const
{
    // Make sure Vgst>0 so we're not in subthreshold mode
    assert(vx < Vddt);
 
    // Check that transistor is actually in triode mode
    // Vds < Vgs - Vth
    assert(vi < Vddt);
 
    // "Snake" voltages for triode mode calculation.
    const unsigned int Vgst = Vddt - vx;
    const unsigned int Vgdt = Vddt - vi;
 
    const unsigned int Vgst_2 = Vgst * Vgst;
    const unsigned int Vgdt_2 = Vgdt * Vgdt;
 
    // "Snake" current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
    const int n_I_snake = n_snake * (static_cast<int>(Vgst_2 - Vgdt_2) >> 15);
 
    // VCR gate voltage.       // Scaled by m*2^16
    // Vg = Vddt - sqrt(((Vddt - Vw)^2 + Vgdt^2)/2)
    const int Vg = static_cast<int>(vcr_Vg[(Vddt_Vw_2 + (Vgdt_2 >> 1)) >> 16]);
    const int Vp = (Vg - nVt) / n; // Pinch-off voltage
    const int kVg = static_cast<int>(Vp) - nVmin;
 
    // VCR voltages for EKV model table lookup.
    const int Vgs = (vx < kVg) ? kVg - vx : 0;
    assert(Vgs < (1 << 16));
    const int Vgd = (vi < kVg) ? kVg - vi : 0;
    assert(Vgd < (1 << 16));
 
    // VCR current, scaled by m*2^15*2^15 = m*2^30
    const unsigned int If = static_cast<unsigned int>(vcr_n_Ids_term[Vgs]) << 15;
    const unsigned int Ir = static_cast<unsigned int>(vcr_n_Ids_term[Vgd]) << 15;
    const int n_I_vcr = (If - Ir) * n;
 
#ifdef SLOPE_FACTOR
    // estimate new slope factor based on gate voltage
    const double gamma = 1.0;   // body effect factor
    const double phi = 0.8;     // bulk Fermi potential
    const double Ut = 26.0e-3;  // Thermal voltage
    const double nVp = Vp / N16;
    n = 1. + (gamma / (2 * sqrt(nVp + phi + 4*Ut)));
    assert((n > 1.2) && (n < 1.8));
#endif
 
    // Change in capacitor charge.
    vc += n_I_snake + n_I_vcr;
 
    // vx = g(vc)
    const int tmp = (vc >> 15) + (1 << 15);
    assert(tmp < (1 << 16));
    vx = opamp_rev[tmp];
 
    // Return vo.
    return vx - (vc >> 14);
}
 
} // namespace reSIDfp
 
#endif
 
#endif