As a sequel of this question, I built a mathematical model for the Id current of the following (same) JFET circuit:

enter image description here

For a simplified model, I took LAMBDA=0 and set the internal resistances Rd=0 and Rs=0. I then proceeded to calculate the roots for Vg_limit (the value of the input voltage that takes the JFET to the limit between saturation and triode regions) and Id_sat (the drain current in the saturation region). For simplification, I didn't consider the triode region in this analysis.

After doing some math, the roots for the Vg and the equation for Id_sat are obtained as:

$$Vg_{limit} = \dfrac{2 \beta R_{ds} V_{dd} \pm \sqrt{4 \beta R_{ds} V_{dd} + 1} + 1}{2 \beta R_{ds}^2} $$

$$Id_{sat} = \dfrac{2 \beta R_{s} V_{gt} \pm \sqrt{2 \beta R_{S} V_{gt} + 1} + 1}{2 \beta R_{s}^2} $$

where $$ R_{ds} = R_d + R_s \quad and \quad V_{gt} = V_g - V_{to}$$

Out of the two roots for each equation above, only one makes sense physically.

I tested this for many NJF models from LTspice and, while it worked for quite a few, I faced some discrepancies for some other devices, in which the BETA was particularly low. As an example, if I take a simplified version of the 2N4118A, you can see from the plots below that Id_sat from the model is totally different than the one from SPICE. Moreover, the Vds - Vgst curve doesn't touch the origin, which should happen, as the circuit is fed with a sine wave having amplitude equal to Vg_limit.

enter image description here

enter image description here

Can you maybe spot what I'm doing wrong here? As I mentioned, this works for other devices, in which the BETA is larger. As a matter of fact, if I replace the value in this model with a higher one (0.4m or higher, instead of 0.095m), the plots start to match.

Another funny fact that I noticed was that, when increasing the temperatures to absurdly high values, for instance 250ºC, the plots also start to match. Strange indeed, as I forced BETATCE=0 and VTOTC=0.

Here is the PySpice code:

import numpy as np
import matplotlib.pyplot as plt
from PySpice.Spice.Netlist import Circuit
from prefixed import Float

def main():
    F0 = 20  # [Hz]
    FS = 48000  # [Hz]
    temp_c = 25  # [celsius]
    T = 0.1  # [seconds]

    # Ground truth: spice model
    circuit = Circuit('JFET')
    my2N4118A = circuit.model("my2N4118A", "NJF",
    # component values
    VDD = 9.0
    RG = 1E6
    RD = 4.4E3
    RS = 1E3
    RDS = RD+RS

    jfet_spice_params = {
        p.upper(): my2N4118A._parameters[p]
        for p in my2N4118A._parameters

    BETA = Float(jfet_spice_params["BETA"])
    VTO = Float(jfet_spice_params["VTO"])

    # Roots for ID at the limit region, where
    # ID - BETA * (VDS)**2 = ID - BETA * (VDD - ID * RDS)**2 = 0
    id_limit = np.array([
        (2*BETA*RDS*VDD - np.sqrt(4*BETA*RDS*VDD + 1) + 1)/(2*BETA*(RDS**2)),
        (2*BETA*RDS*VDD + np.sqrt(4*BETA*RDS*VDD + 1) + 1)/(2*BETA*(RDS**2))
    vd = VDD - id_limit * RD
    vg = vd + VTO
    vs = id_limit * RS
    vds = vd - vs
    if (id_limit[0] > 0 and vds[0] > 0) and (id_limit[1] <= 0 or vds[1] <= 0):
        vg_limit = vg[0]
    elif (id_limit[1] > 0 and vds[1] > 0) and (id_limit[0] <= 0 or vds[0] <= 0):
        vg_limit = vg[1]
        vg_limit = None
        print("Error! VG limit sanity check failed")

    # Netlist
    circuit.V('Vdd', 'vdd', circuit.gnd, VDD)
        'in', 'gate', circuit.gnd, amplitude=vg_limit, frequency=F0)
    circuit.R('Rg', 'gate', circuit.gnd, RG)
    circuit.R('Rd', 'vdd', 'drain', RD)
    circuit.R('Rs', 'source', circuit.gnd, RS)
    circuit.J('JFET', 'drain', 'gate', 'source', model="my2N4118A")

    plt.title('Id x t')

    # Spice model
    simulator = circuit.simulator(
        temperature=temp_c, nominal_temperature=temp_c)
    # analysis = simulator.dc(Vin=Vsl)
    analysis = simulator.transient(step_time=1/FS, end_time=T)
    vin_spice = np.array(analysis["gate"])
    vgs_spice = np.array(analysis["gate"]) - np.array(analysis["source"])
    id_spice = (np.array(analysis["vdd"]) - np.array(analysis["drain"])) / RD
    vgst_spice = vgs_spice - VTO
    vds_spice = np.array(analysis["drain"]) - np.array(analysis["source"])
    # if > 0 then saturation, else triode
    v_limit_spice = vds_spice - vgst_spice
    t = np.array(analysis.time)

    # Roots for IDsat, when LAMBDA = 0
    vgt = vin_spice - VTO
    BETA_2 = 2 * BETA
    BETA_2_RS_VGT = BETA_2 * RS * vgt
    # Out of the 2 existent roots, only the one with "- np.sqrt(...)" makes sense physically
    id_sat = (BETA_2_RS_VGT - np.sqrt(2 * BETA_2_RS_VGT + 1) + 1) / \
    id_sat[vgt <= 0] = 0

    plt.plot(t, id_spice, label="Id_spice")
    plt.plot(t, id_sat, label="Id_model")

    plt.xlabel('T [s]')
    plt.ylabel('Id [A]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')

    plt.title('V x t')
    plt.plot(t, vgst_spice, label=f"Vgst")
    plt.plot(t, v_limit_spice, label=f"Vds-Vgst")
    plt.xlabel('T [s]')
    plt.ylabel('V [V]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')


if __name__ == "__main__":

  • \$\begingroup\$ Garbage in, garbage out. It is meaningless to set the temperature to 250C when the physical device cannot tolerate it. You know the result at 250C should be that the device no longer functions and if the simulation indicates otherwise you know you are looking at garbage. A more complex model will have indicated it fails, and an even more complex model may attempt to simulate the malfunction but a simpler model will not and just continue as usual. And if you are expecting it simulate the malfunction then you should obviously not be expecting anything "normal" looking. \$\endgroup\$
    – DKNguyen
    Apr 30, 2023 at 17:06
  • \$\begingroup\$ Could be. The 250ºC simulation was a side-note of something I noticed when playing around with the temperature. I found it curious that it "matches" the model for such an extreme condition. Fact is, under normal temperatures, there is a mismatch. \$\endgroup\$
    – edwillys
    Apr 30, 2023 at 18:09
  • \$\begingroup\$ My hunch is that you might still be caught in the same pattern of thinking that got you stuck on your previous question. Are you sure you're accurately duplicating what SPICE is doing? For example, I don't see the gate-source and gate-drain PN junctions implemented anywhere in your Python model. Did you take a look at that SPICE2 thesis I linked you in the other question (in this case: Figure A2.12)? I would also familiarize yourself with GMIN if you can. \$\endgroup\$
    – Ste Kulov
    Apr 30, 2023 at 18:41
  • 1
    \$\begingroup\$ @SteKulov You were right on the hint about the PN junction diodes. I dissected the thesis and the Semiconductor Device Modeling With Spice, 2nd ed. and indeed for the devices in which I saw the mismatch, the gate-source junction was starting to conduct, as Vgs>0. Bottom line for the model: I need to calculate another threshold input voltage, above which Vgs>0. For my needs, this is actually not a valid use case, so I'll just limit the input range so that it doesn't reach that point. \$\endgroup\$
    – edwillys
    May 1, 2023 at 13:52
  • 1
    \$\begingroup\$ @edwillys Sorry, I'll clarify. What I meant by the "pattern of thinking" comment was that you're applying semiconductor models generally, but not what's specifically implemented in SPICE. To know what it's SPICE, you have to see what's actually implemented in SPICE. The Massobrio book and Nagel thesis are great starting points. But keep in mind, depending on how deep you go, you might have to actually dig down into ngspice source code (which I've had to do before) to see what's going on find where discrepancies originate from. \$\endgroup\$
    – Ste Kulov
    May 1, 2023 at 21:16

1 Answer 1


It turns out that I didn't consider the Id behavior when Vgs > 0. In this region, the gate-source PN junction (as hinted by Ste Kulov) starts to conduct and the Igs starts to play a major role in the circuit, effectively clamping Id as Vin further increases.

For my needs, the use-case when Vgs > 0 is not a valid one. Hence, I just need to limit my analysis for an input voltage below a threshold Vin_limit_pn, where Vgs=0. Hence, as the JFET is in the saturation region:

$$ I_{d\_limit\_pn} = \beta \times (0+V_{TO}^2) $$ $$ V_{in\_limit\_pn} = V_{gs} + V_s = 0 + R_s \times I_{d\_limit\_pn} = R_S \times \beta \times V_{TO}^2 $$

  • 1
    \$\begingroup\$ Good workaround. Nice job! \$\endgroup\$
    – Ste Kulov
    May 1, 2023 at 21:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.