JFET Id calculation mismatch in SPICE

Question

I'm simulating the following circuit in PySpice:

^{simulate this circuit – Schematic created using CircuitLab}

In this configuration, the N-JFET is always in saturation mode, as you can see by the plot of Vds - (Vgs - Vto) = Vds - Vgt, which is always positive:

Hence, the value of Id should be given by the following equation, according to the S&H model:

$$ I_d = 0.5 \beta (V_{gs} - V_t)^2 (1 + \lambda V_{ds}) $$

However, when I compare the value of the measured Id = (Vdd - Vd) / Rd to the one calculated from the formula above, I get the following mismatch:

Could anybody help me on what I missed here? Below you can find the Python 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 = 20  # [celsius]
    A = 1.0  # [V] amplitude, linear
    T = 0.1  # [seconds]

    # Spice circuit, based on the LSK189A from LTSpice
    circuit = Circuit('JFET')
    LSK189A = circuit.model("LSK189A", "NJF",
                            Beta="2.2m",
                            Betatce="-0.5",
                            Vto="-1.13",
                            Vtotc="-2.5m",
                            Lambda="4.3m",
                            Is="3f",
                            Xti="0",
                            Isr="0",
                            Alpha="30u",
                            N="1",
                            Rd="11",
                            Rs="30",
                            Cgd="3.19p",
                            Cgs="2.92p",
                            Fc="0.5",
                            Vk="120",
                            M="320m",
                            Pb="0.8",
                            Kf="0.0009f",
                            Af="1",
                            Gdsnoi="2.15",
                            Nlev="1",  # changed from original 3
                            Mfg="Linear_Systems",
                            )
    # component values
    Vdd = 9.0
    Rg = 1E6
    Rd = 4.4E3
    Rs = 1E3

    # Netlist
    circuit.V('Vdd', 'vdd', circuit.gnd, Vdd)
    circuit.SinusoidalVoltageSource(
        'in', 'gate', circuit.gnd, amplitude=A, 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="LSK189A")


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

    plt.figure()
    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)
    y_id_spice = (np.array(analysis["vdd"]) - np.array(analysis["drain"])) / Rd
    y_vgt_spice = np.array(
        analysis["gate"]) - np.array(analysis["source"]) - Float(jfet_spice_params["VTO"])
    y_vds_spice = np.array(analysis["drain"]) - np.array(analysis["source"])
    # if > 0 then saturation
    y_limit_spice = y_vds_spice - y_vgt_spice
    t = np.array(analysis.time)

    plt.plot(t, y_id_spice, label="Id_spice")
    Id_spice_calc = (1 + Float(jfet_spice_params["LAMBDA"]) * y_vds_spice) * 0.5 * Float(
        jfet_spice_params["BETA"]) * y_vgt_spice**2
    plt.plot(t, Id_spice_calc, label="Id_spice_calc")

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

    plt.figure()
    plt.title('V x t')
    plt.plot(t, y_vgt_spice, label=f"Vgt")
    plt.plot(t, y_limit_spice, label=f"Vds-Vgt")
    plt.xlabel('T [s]')
    plt.ylabel('V [V]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.legend()

    plt.show()


if __name__ == "__main__":
    main()

Answer 1

The 0.5 factor should not be there in your calculation. See 3-22 page 143 of 2nd edition Semiconductor Device Modeling With Spice (Guiseppe et al.) which is the same equation as you have but the 0.5 has been rolled into \$\beta\$.

The remaining discrepancy in the static model is mostly accounted for by Rs (30Ω in the 'A' bin model), which reduces Vgs.

I presume the 'B' model is intended to be closer to typical with 'A' and 'C' representing outliers..