diff options
author | 2019-04-02 15:11:11 +0200 | |
---|---|---|
committer | 2019-04-07 19:08:33 +0200 | |
commit | 140dc2237eed85e3451ab95a4bcdae84ee8600c5 (patch) | |
tree | f25eb489043d3cb008b061f45233df17ca977abc | |
parent | 24c53ebe1b59c04ebc7a4719c799d8381770633f (diff) |
netlist: TYPE=2 opamp model and MB3614 fixes.
- Untie diode code
- Fix some typos
- add TYPE=2 opamp model which omits output voltage limitation
Useful for determining causes of non-convergence
- Fix MB3614 opamp parameters to match datasheet
-rw-r--r-- | nl_examples/cmos_inverter_clk.cpp | 10 | ||||
-rw-r--r-- | nl_examples/congo_bongo.c | 34 | ||||
-rw-r--r-- | nl_examples/opamp.c | 4 | ||||
-rw-r--r-- | src/lib/netlist/analog/nld_bjt.cpp | 7 | ||||
-rw-r--r-- | src/lib/netlist/analog/nld_generic_models.h | 89 | ||||
-rw-r--r-- | src/lib/netlist/analog/nld_opamps.cpp | 64 | ||||
-rw-r--r-- | src/lib/netlist/analog/nlid_twoterm.h | 3 | ||||
-rw-r--r-- | src/lib/netlist/macro/nlm_opamp.cpp | 2 | ||||
-rw-r--r-- | src/lib/netlist/plib/putil.h | 7 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.cpp | 5 |
10 files changed, 140 insertions, 85 deletions
diff --git a/nl_examples/cmos_inverter_clk.cpp b/nl_examples/cmos_inverter_clk.cpp index 09064a59d19..95dba92bb79 100644 --- a/nl_examples/cmos_inverter_clk.cpp +++ b/nl_examples/cmos_inverter_clk.cpp @@ -17,18 +17,18 @@ NETLIST_START(cmos_inverter_clk) //SOLVER(Solver, 480000) SOLVER(Solver, 1e5) //SOLVER(Solver, 100000000000) - PARAM(Solver.ACCURACY, 1e-7 ) + PARAM(Solver.ACCURACY, 1e-6 ) PARAM(Solver.NR_LOOPS, 500000) - PARAM(Solver.DYNAMIC_TS, 0) + PARAM(Solver.DYNAMIC_TS, 1) PARAM(Solver.DYNAMIC_LTE, 1e-5) - PARAM(Solver.DYNAMIC_MIN_TIMESTEP, 2e-8) + PARAM(Solver.DYNAMIC_MIN_TIMESTEP, 2e-9) ANALOG_INPUT(V5, 5) // CLOCK(clk, 0.5e6) #if (USE_CLOCK) CLOCK(V, 5000) - PARAM(NETLIST.DEFAULT_MOS_CAPMODEL, 0) // Disable capacitance modeling + PARAM(NETLIST.DEFAULT_MOS_CAPMODEL, 2) // Disable capacitance modeling //CLOCK(V, 500000) #else VS(V, 5) @@ -56,7 +56,7 @@ NETLIST_START(cmos_inverter_clk) NET_C(M.D, C.1) NET_C(M.S, C.2) #endif -#if 0 +#if 1 LOG(log_G, M.G) LOG(log_D, M.D) LOGD(log_X, RG.1, RG.2) diff --git a/nl_examples/congo_bongo.c b/nl_examples/congo_bongo.c index 492865c79c4..c31258b0bb2 100644 --- a/nl_examples/congo_bongo.c +++ b/nl_examples/congo_bongo.c @@ -43,27 +43,28 @@ NETLIST_START(dummy) // IGNORED O_AUDIO0: O_AUDIO0 64 0 // .END - PARAM(Solver.ACCURACY, 1e-8) - PARAM(Solver.NR_LOOPS, 90) - PARAM(Solver.SOR_FACTOR, 1.01) - PARAM(Solver.GS_LOOPS, 4) + PARAM(Solver.ACCURACY, 1e-7) + PARAM(Solver.NR_LOOPS, 90000) + PARAM(Solver.SOR_FACTOR, 1.0) + PARAM(Solver.GS_LOOPS, 99) //PARAM(Solver.METHOD, "GMRES") PARAM(Solver.METHOD, "MAT_CR") //PARAM(Solver.METHOD, "SOR") #if USE_OPTMIZATIONS #if USE_FRONTIERS - SOLVER(Solver, 24000) + SOLVER(Solver, 48000) #else SOLVER(Solver, 48000) #endif - PARAM(Solver.DYNAMIC_TS, 0 ) + PARAM(Solver.DYNAMIC_TS, 0) PARAM(Solver.PARALLEL, 1) + PARAM(Solver.DYNAMIC_MIN_TIMESTEP, 2e-6) #else - SOLVER(Solver, 24000) + SOLVER(Solver, 48000) PARAM(Solver.DYNAMIC_TS, 1) PARAM(Solver.DYNAMIC_LTE, 1e-4) - PARAM(Solver.DYNAMIC_MIN_TIMESTEP, 5e-7) + PARAM(Solver.DYNAMIC_MIN_TIMESTEP, 1e-7) PARAM(Solver.PARALLEL, 0) PARAM(Solver.PIVOT, 0) #endif @@ -91,8 +92,9 @@ NETLIST_START(dummy) NET_C(CO.2, RO.1) NET_C(RO.2, GND) - // FIXME: Same as 1N4148 - NET_MODEL("1S2075 D(Is=2.52n Rs=.568 N=1.752 Cjo=4p M=.4 tt=20n Iave=200m Vpk=75)") + // Found here: https://hamesspam.sakura.ne.jp/hes2016/160521.html + NET_MODEL("1S2075 D(IS=1.387E-9 N=1.702 RS=1.53 CJO=1.92pf VJ=0.4996 M=0.0605 TT=5ns BV=75 IBV=100E-15)") + NET_MODEL("2SC1941 NPN(IS=46.416f BF=210 NF=1.0022 VAF=600 IKF=500m ISE=60f NE=1.5 BR=2.0122 NR=1.0022 VAR=10G IKR=10G ISC=300p NC=2 RB=13.22 IRB=10G RBM=13.22 RE=100m RC=790m CJE=26.52p VJE=900m MJE=518m TF=1.25n XTF=10 VTF=10 ITF=500m PTF=0 CJC=4.89p VJC=750m MJC=237m XCJC=500m TR=100n CJS=0 VJS=750m MJS=500m XTB=1.5 EG=1.11 XTI=3 KF=0 AF=1 FC=500m)") INCLUDE(CongoBongo_schematics) @@ -103,16 +105,18 @@ NETLIST_START(dummy) * to the feedback loop consisting of D9 and Q2. */ +#if 0 RES(RX1, 1e100) NET_C(RX1.1, Q2.C) NET_C(RX1.2, XU16.7) +#endif - /* The opamp actually has an FPF of about 500k. This doesn't work here and causes oscillations. + /* The opamp actually has an FPF of about 200k. This doesn't work here and causes oscillations. * FPF here therefore about half the Solver clock. */ - PARAM(XU16.B.MODEL, "MB3614(UGF=11k)") - PARAM(XU17.C.MODEL, "MB3614(UGF=11k)") - + PARAM(XU16.B.MODEL, "MB3614(TYPE=3 UGF=22k)") + PARAM(XU17.C.MODEL, "MB3614(TYPE=3 UGF=22k)") +#if 0 PARAM(XU17.A.MODEL, "MB3614(TYPE=1)") PARAM(XU17.B.MODEL, "MB3614(TYPE=1)") PARAM(XU17.D.MODEL, "MB3614(TYPE=1)") @@ -120,7 +124,7 @@ NETLIST_START(dummy) //PARAM(XU16.A.MODEL, "MB3614(TYPE=1)") PARAM(XU16.C.MODEL, "MB3614(TYPE=1)") PARAM(XU16.D.MODEL, "MB3614(TYPE=1)") - +#endif #if 0 PARAM(XU13.A.MODEL, "MB3614(TYPE=1)") PARAM(XU13.B.MODEL, "MB3614(TYPE=1)") diff --git a/nl_examples/opamp.c b/nl_examples/opamp.c index d3dc80f76a9..b13d92dab67 100644 --- a/nl_examples/opamp.c +++ b/nl_examples/opamp.c @@ -50,8 +50,8 @@ NETLIST_START(main) NET_C(RL.2, GND) NET_C(RL.1, op1.OUT) - LOG(logX, op1.OUT) - LOG(logY, clk) + LOG(log_X, op1.OUT) + LOG(log_Y, clk) NETLIST_END() NETLIST_START(opamp) diff --git a/src/lib/netlist/analog/nld_bjt.cpp b/src/lib/netlist/analog/nld_bjt.cpp index a174b3876b2..98a0d80836f 100644 --- a/src/lib/netlist/analog/nld_bjt.cpp +++ b/src/lib/netlist/analog/nld_bjt.cpp @@ -49,7 +49,7 @@ namespace analog // nld_Q - Base classes // ----------------------------------------------------------------------------- - /*! Class representing the bjt model paramers. + /*! Class representing the bjt model parameters. * * This is the model representation of the bjt model. Typically, SPICE uses * the following parameters. A "Y" in the first column indicates that the @@ -71,7 +71,7 @@ namespace analog * | | ISC | leakage saturation current | A | 0 | 8 | | * | | NC | leakage emission coefficient | - | 2 | 1.5 | | * | | RB | zero bias base resistance | | 0 | 100 | * | - * | | IRB | current where base resistance falls halfway to its min value | A | infinte | 0.1 | * | + * | | IRB | current where base resistance falls halfway to its min value | A | infinite | 0.1 | * | * | | RBM | minimum base resistance at high currents | | RB | 10 | * | * | | RE | emitter resistance | | 0 | 1 | * | * | | RC | collector resistance | | 0 | 10 | * | @@ -271,6 +271,7 @@ namespace analog connect("B", "m_CJC.1"); connect("C", "m_CJC.2"); } + } protected: @@ -428,6 +429,8 @@ namespace analog const nl_double Ie = (sIe + gee * m_gD_BE.Vd() - gec * m_gD_BC.Vd()) * polarity; const nl_double Ic = (sIc - gce * m_gD_BE.Vd() + gcc * m_gD_BC.Vd()) * polarity; + // "Circuit Design", page 174 + m_D_EB.set_mat( gee, gec - gee, -Ie, gce - gee, gee - gec, Ie); m_D_CB.set_mat( gcc, gce - gcc, -Ic, diff --git a/src/lib/netlist/analog/nld_generic_models.h b/src/lib/netlist/analog/nld_generic_models.h index 6bb8b815baf..c18452e1109 100644 --- a/src/lib/netlist/analog/nld_generic_models.h +++ b/src/lib/netlist/analog/nld_generic_models.h @@ -47,17 +47,25 @@ namespace analog capacitor_e type() const { return capacitor_e::VARIABLE_CAPACITY; } + // Circuit Simulation, page 284, 5.360 + // q(un+1) - q(un) = int(un, un+1, C(U)) = (C0+C1)/2 * (un+1-un) + // The direct application of formulas 5.359 and 5.360 has + // issues with pulses. Therefore G and Ieq are expressed differently + // so that G depends on un+1 only and Ieq on un only. + // In both cases, i = G * un+1 + Ieq + nl_double G(nl_double cap) const { - //return m_h * (2.0 * cap - m_c) + m_gmin; + //return m_h * cap + m_gmin; return m_h * 0.5 * (cap + m_c) + m_gmin; + //return m_h * cap + m_gmin; } nl_double Ieq(nl_double cap, nl_double v) const { - plib::unused_var(v); - // return m_h * (cap * v - m_charge) - G(cap) * v; + //return -m_h * 0.5 * ((cap + m_c) * m_v + (cap - m_c) * v) ; return -m_h * 0.5 * (cap + m_c) * m_v; + //return -m_h * cap * m_v; } void timestep(nl_double cap, nl_double v, nl_double step) @@ -76,28 +84,31 @@ namespace analog nl_double m_gmin; }; + // "Circuit simulation", page 274 template <> class generic_capacitor<capacitor_e::CONSTANT_CAPACITY> { public: generic_capacitor(device_t &dev, const pstring &name) : m_h(dev, name + ".m_h", 0.0) + , m_v(dev, name + ".m_v", 0.0) , m_gmin(0.0) { } capacitor_e type() const { return capacitor_e::CONSTANT_CAPACITY; } nl_double G(nl_double cap) const { return cap * m_h + m_gmin; } - nl_double Ieq(nl_double cap, nl_double v) const { return - G(cap) * v; } + nl_double Ieq(nl_double cap, nl_double v) const { return - G(cap) * m_v; } void timestep(nl_double cap, nl_double v, nl_double step) { - plib::unused_var(cap, v); m_h = 1.0 / step; + m_v = v; } void setparams(nl_double gmin) { m_gmin = gmin; } private: state_var<nl_double> m_h; + state_var<double> m_v; nl_double m_gmin; }; @@ -127,6 +138,7 @@ namespace analog , m_gmin(1e-15) , m_VtInv(0.0) , m_Vcrit(0.0) + , m_name(name) { set_param(1e-15, 1, 1e-15, 300.0); } @@ -135,37 +147,46 @@ namespace analog { nl_double IseVDVt(0.0); - if (TYPE == diode_e::BIPOLAR && nVd < m_Vmin) + if (TYPE == diode_e::BIPOLAR) { - m_Vd = nVd; - m_G = m_gmin; - m_Id = - m_Is; + //printf("%s: %g %g\n", m_name.c_str(), nVd, (double) m_Vd); + if (nVd > m_Vcrit) + { + const nl_double d = std::min(1e100, nVd - m_Vd); + const nl_double a = std::abs(d) * m_VtInv; + m_Vd = m_Vd + (d < 0 ? -1.0 : 1.0) * std::log1p(a) * m_Vt; + } + else + m_Vd = std::max(-1e100, nVd); + //m_Vd = nVd; + + if (m_Vd < m_Vmin) + { + m_G = m_gmin; + m_Id = - m_Is; + } + else + { + IseVDVt = std::exp(m_logIs + m_Vd * m_VtInv); + m_Id = IseVDVt - m_Is; + m_G = IseVDVt * m_VtInv + m_gmin; + } } - else if (TYPE == diode_e::MOS && nVd < constants::zero()) + else if (TYPE == diode_e::MOS) { - m_Vd = nVd; - m_G = m_Is * m_VtInv + m_gmin; - m_Id = m_G * m_Vd; - } - else if (/*TYPE == diode_e::MOS || */nVd < m_Vcrit) - { - m_Vd = nVd; - IseVDVt = std::exp(std::min(300.0, m_logIs + m_Vd * m_VtInv)); - m_Id = IseVDVt - m_Is; - m_G = IseVDVt * m_VtInv + m_gmin; - } - else - { - if (TYPE == diode_e::MOS && m_Vd < constants::zero()) - m_Vd = std::min(m_Vmin, nVd); - - const nl_double d = (nVd - m_Vd); - const nl_double a = std::abs(nVd - m_Vd) * m_VtInv; - m_Vd = m_Vd + (d < 0 ? -1.0 : 1.0) * std::log1p(a) * m_Vt; - IseVDVt = std::exp(m_logIs + m_Vd * m_VtInv); - //const double IseVDVt = m_Is * std::exp(m_Vd * m_VtInv); - m_Id = IseVDVt - m_Is; - m_G = IseVDVt * m_VtInv + m_gmin; + if (nVd < constants::zero()) + { + m_Vd = nVd; + m_G = m_Is * m_VtInv + m_gmin; + m_Id = m_G * m_Vd; + } + else /* log stepping should already be done in mosfet */ + { + m_Vd = nVd; + IseVDVt = std::exp(std::min(300.0, m_logIs + m_Vd * m_VtInv)); + m_Id = IseVDVt - m_Is; + m_G = IseVDVt * m_VtInv + m_gmin; + } } } @@ -207,6 +228,8 @@ namespace analog nl_double m_VtInv; nl_double m_Vcrit; + + pstring m_name; }; diff --git a/src/lib/netlist/analog/nld_opamps.cpp b/src/lib/netlist/analog/nld_opamps.cpp index 8c3c655bfa9..0539a5dea51 100644 --- a/src/lib/netlist/analog/nld_opamps.cpp +++ b/src/lib/netlist/analog/nld_opamps.cpp @@ -115,6 +115,8 @@ namespace netlist , m_VREF(*this, "VREF") { m_type = static_cast<int>(m_model.m_TYPE); + if (m_type < 1 || m_type > 3) + log().fatal(MF_1_UNKNOWN_OPAMP_TYPE, m_type); if (m_type == 1) { @@ -127,38 +129,42 @@ namespace netlist connect("RP1.1", "G1.OP"); } - else if (m_type == 3) + if (m_type == 2 || m_type == 3) { create_and_register_subdevice("CP1", m_CP); create_and_register_subdevice("EBUF", m_EBUF); - create_and_register_subdevice("DN", m_DN, "D(IS=1e-15 N=1)"); - create_and_register_subdevice("DP", m_DP, "D(IS=1e-15 N=1)"); - - //m_DP->m_model.setTo("D(IS=1e-15 N=1)"); - //m_DN->m_model.setTo("D(IS=1e-15 N=1)"); register_subalias("PLUS", "G1.IP"); register_subalias("MINUS", "G1.IN"); - register_subalias("OUT", "EBUF.OP"); - - connect("EBUF.ON", "VREF"); connect("G1.ON", "VREF"); connect("RP1.2", "VREF"); connect("CP1.2", "VREF"); + connect("EBUF.ON", "VREF"); connect("EBUF.IN", "VREF"); connect("RP1.1", "G1.OP"); connect("CP1.1", "RP1.1"); + connect("EBUF.IP", "RP1.1"); + } + if (m_type == 2) + { + register_subalias("OUT", "EBUF.OP"); + } + if (m_type == 3) + { + + create_and_register_subdevice("DN", m_DN, "D(IS=1e-15 N=1)"); + create_and_register_subdevice("DP", m_DP, "D(IS=1e-15 N=1)"); + connect("DP.K", "VH"); connect("VL", "DN.A"); connect("DP.A", "DN.K"); connect("DN.K", "RP1.1"); - connect("EBUF.IP", "RP1.1"); + + register_subalias("OUT", "EBUF.OP"); } - else - log().fatal(MF_1_UNKNOWN_OPAMP_TYPE, m_type); } @@ -194,6 +200,7 @@ namespace netlist const double cVt = 0.0258 * 1.0; // * m_n; const double cId = m_model.m_DAB; // 3 mA const double cVd = cVt * std::log(cId / 1e-15 + 1.0); + m_VH.push(m_VCC() - m_model.m_VLH - cVd); m_VL.push(m_GND() + m_model.m_VLL + cVd); m_VREF.push((m_VCC() + m_GND()) / 2.0); @@ -211,28 +218,41 @@ namespace netlist m_RP.set_R(RO); m_G1.m_G.setTo(G); } - else if (m_type == 3) + if (m_type == 3 || m_type == 2) { m_EBUF->reset(); - m_DP->reset(); - m_DN->reset(); m_CP->reset(); m_RP.reset(); - - m_EBUF->m_G.setTo(1.0); - m_EBUF->m_RO.setTo(m_model.m_RO); - +#if 0 + double CP = 0.0; + double RP = 1000; + double G = m_model.m_UGF / m_model.m_FPF / RP; +#else double CP = m_model.m_DAB / m_model.m_SLEW; - double RP = 0.5 / 3.1459 / CP / m_model.m_FPF; + double RP = 0.5 / constants::pi() / CP / m_model.m_FPF; double G = m_model.m_UGF / m_model.m_FPF / RP; - - //printf("Min Freq %s: %f\n", name().c_str(), 1.0 / (CP*RP / (G*RP))); +#endif + printf("OPAMP %s: %g %g %g\n", name().c_str(), CP, RP, G); + if (m_model.m_SLEW / (4.0 * constants::pi() * 0.0258) < m_model.m_UGF) + printf("failed!\n"); m_CP->m_C.setTo(CP); m_RP.set_R(RP); m_G1.m_G.setTo(G); } + if (m_type == 2) + { + m_EBUF->m_G.setTo(1.0); + m_EBUF->m_RO.setTo(m_model.m_RO); + } + if (m_type == 3) + { + m_DP->reset(); + m_DN->reset(); + m_EBUF->m_G.setTo(1.0); + m_EBUF->m_RO.setTo(m_model.m_RO); + } } diff --git a/src/lib/netlist/analog/nlid_twoterm.h b/src/lib/netlist/analog/nlid_twoterm.h index 8e93d1a6fba..63ab43db591 100644 --- a/src/lib/netlist/analog/nlid_twoterm.h +++ b/src/lib/netlist/analog/nlid_twoterm.h @@ -313,7 +313,7 @@ namespace analog * |NL? |name |parameter |units|default| example|area | * |:--:|:-----|:--------------------------------|:----|------:|-------:|:----:| * | Y |IS |saturation current |A |1.0e-14| 1.0e-14| * | - * | |RS |ohmic resistanc |Ohm | 0| 10| * | + * | |RS |ohmic resistance |Ohm | 0| 10| * | * | Y |N |emission coefficient |- | 1| 1| | * | |TT |transit-time |sec | 0| 0.1ns| | * | |CJO |zero-bias junction capacitance |F | 0| 2pF| * | @@ -359,6 +359,7 @@ namespace analog register_subalias("K", m_N); } + /* needed by nld_opamp ... */ template <class CLASS> NETLIB_NAME(D)(CLASS &owner, const pstring &name, const pstring &model) : NETLIB_NAME(twoterm)(owner, name) diff --git a/src/lib/netlist/macro/nlm_opamp.cpp b/src/lib/netlist/macro/nlm_opamp.cpp index 849f7820046..2c61430c935 100644 --- a/src/lib/netlist/macro/nlm_opamp.cpp +++ b/src/lib/netlist/macro/nlm_opamp.cpp @@ -225,7 +225,7 @@ NETLIST_START(OPAMP_lib) NET_MODEL("LM324 OPAMP(TYPE=3 VLH=2.0 VLL=0.2 FPF=5 UGF=500k SLEW=0.3M RI=1000k RO=50 DAB=0.00075)") NET_MODEL("LM358 OPAMP(TYPE=3 VLH=2.0 VLL=0.2 FPF=5 UGF=500k SLEW=0.3M RI=1000k RO=50 DAB=0.001)") - NET_MODEL("MB3614 OPAMP(TYPE=3 VLH=1.4 VLL=0.02 FPF=2 UGF=500k SLEW=0.6M RI=1000k RO=50 DAB=0.0002)") + NET_MODEL("MB3614 OPAMP(TYPE=3 VLH=1.4 VLL=0.02 FPF=2 UGF=200k SLEW=0.6M RI=1000k RO=50 DAB=0.002)") NET_MODEL("UA741 OPAMP(TYPE=3 VLH=1.0 VLL=1.0 FPF=5 UGF=1000k SLEW=0.5M RI=2000k RO=75 DAB=0.0017)") NET_MODEL("LM747 OPAMP(TYPE=3 VLH=1.0 VLL=1.0 FPF=5 UGF=1000k SLEW=0.5M RI=2000k RO=50 DAB=0.0017)") NET_MODEL("LM747A OPAMP(TYPE=3 VLH=2.0 VLL=2.0 FPF=5 UGF=1000k SLEW=0.7M RI=6000k RO=50 DAB=0.0015)") diff --git a/src/lib/netlist/plib/putil.h b/src/lib/netlist/plib/putil.h index 5cc4fc53187..7ec668a1f8f 100644 --- a/src/lib/netlist/plib/putil.h +++ b/src/lib/netlist/plib/putil.h @@ -62,10 +62,11 @@ namespace plib template <typename T> struct constants { - static constexpr T zero() noexcept { return static_cast<T>(0); } - static constexpr T one() noexcept { return static_cast<T>(1); } - static constexpr T two() noexcept { return static_cast<T>(2); } + static constexpr T zero() noexcept { return static_cast<T>(0); } + static constexpr T one() noexcept { return static_cast<T>(1); } + static constexpr T two() noexcept { return static_cast<T>(2); } static constexpr T sqrt2() noexcept { return static_cast<T>(1.414213562373095048801688724209); } + static constexpr T pi() noexcept { return static_cast<T>(3.14159265358979323846264338327950); } /*! * \brief Electric constant of vacuum diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp index 0a7bf5ef3c3..1b92017c284 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.cpp +++ b/src/lib/netlist/solver/nld_matrix_solver.cpp @@ -570,9 +570,12 @@ namespace devices analog_net_t *n = m_nets[k]; terms_for_net_t *t = m_terms[k].get(); - const nl_double DD_n = (n->Q_Analog() - t->m_last_V); + //const nl_double DD_n = (n->Q_Analog() - t->m_last_V); + // avoid floating point exceptions + const nl_double DD_n = std::max(-1e100, std::min(1e100,(n->Q_Analog() - t->m_last_V))); const nl_double hn = cur_ts; + //printf("%g %g %g %g\n", DD_n, hn, t->m_DD_n_m_1, t->m_h_n_m_1); nl_double DD2 = (DD_n / hn - t->m_DD_n_m_1 / t->m_h_n_m_1) / (hn + t->m_h_n_m_1); nl_double new_net_timestep; |