//=============================================================================
// Description: IBIS macro library implemented in Verilog-A(MS)
// Created: September 2005 by Arpad Muranyi (Intel Corporation)
//=============================================================================
// List of building blocks contained in this library:
//
// IBIS_R             (p, n);
// IBIS_VCR           (p, n, ps, ns);
// IBIS_CCR           (p, n, ps, ns);
//
// IBIS_C             (p, n);
// IBIS_VCC           (p, n, ps, ns);       // Commented out, need to fix it
// IBIS_CCC           (p, n, ps, ns);       // Commented out, need to fix it
//
// IBIS_L             (p, n);
// IBIS_VCL           (p, n, ps, ns);       // Commented out, need to fix it
// IBIS_CCL           (p, n, ps, ns);       // Commented out, need to fix it
// IBIS_K             (p1, n1, p2, n2);
//
// IBIS_V             (p, n);
// IBIS_VCVS          (p, n, ps, ns);
// IBIS_CCVS          (p, n, ps, ns);
// IBIS_VCVS_DELAY    (p, n, ps, ns);
// IBIS_CCVS_DELAY    (p, n, ps, ns);
// IBIS_VCVS_SUM      (p, n, ps1, ns1, ps2, ns2);
// IBIS_CCVS_SUM      (p, n, ps1, ns1, ps2, ns2);
// IBIS_VCVS_MULT     (p, n, ps1, ns1, ps2, ns2);
// IBIS_CCVS_MULT     (p, n, ps1, ns1, ps2, ns2);
// IBIS_VCVS_DIV      (p, n, ps1, ns1, ps2, ns2);
// IBIS_CCVS_DIV      (p, n, ps1, ns1, ps2, ns2);
//
// IBIS_I             (p, n);
// IBIS_VCCS          (p, n, ps, ns);
// IBIS_CCCS          (p, n, ps, ns);
// IBIS_VCCS_DELAY    (p, n, ps, ns);
// IBIS_CCCS_DELAY    (p, n, ps, ns);
// IBIS_VCCS_SUM      (p, n, ps1, ns1, ps2, ns2);
// IBIS_CCCS_SUM      (p, n, ps1, ns1, ps2, ns2);
// IBIS_VCCS_MULT     (p, n, ps1, ns1, ps2, ns2);
// IBIS_CCCS_MULT     (p, n, ps1, ns1, ps2, ns2);
// IBIS_VCCS_DIV      (p, n, ps1, ns1, ps2, ns2);
// IBIS_CCCS_DIV      (p, n, ps1, ns1, ps2, ns2);
//
// IBIS_T             (n1, ref1, n2, ref2);
//
// IBIS_Input         (PC_ref, GC_ref, Input, Rcv_D);
// IBIS_Output        (PU_ref, PD_ref, Output, In_D, PC_ref, GC_ref);
// IBIS_IO            (PU_ref, PD_ref, IO, In_D, En_D, Rcv_D, PC_ref, GC_ref);
// IBIS_3state        (PU_ref, PD_ref, IO, In_D, En_D, PC_ref, GC_ref);
// IBIS_OpenSink      (PU_ref, PD_ref, IO, In_D, En_D, PC_ref, GC_ref);
// IBIS_IO_OpenSink   (PU_ref, PD_ref, IO, In_D, En_D, Rcv_D, PC_ref, GC_ref);
// IBIS_OpenSource    (PU_ref, PD_ref, IO, In_D, En_D, PC_ref, GC_ref);
// IBIS_IO_OpenSource (PU_ref, PD_ref, IO, In_D, En_D, Rcv_D, PC_ref, GC_ref);
//-----------------------------------------------------------------------------
// List of building blocks to be implemented:
//
// IBIS_BUFFER_SERIES
// IBIS_SPARAM
//
//|                 Input, Output, I/O, 3-state, Open_drain, I/O_open_drain,
//|                 Open_sink, I/O_open_sink, Open_source, I/O_open_source,
//|                 Input_ECL, Output_ECL, I/O_ECL, 3-state_ECL, Terminator,
//|                 Series, and Series_switch.
//-----------------------------------------------------------------------------
// Other ideas:
//
// More than 2 element mutual inductor(s)?
// More than 2 input adder, multiplier?
// Min(x,y) and Max(x,y) function?
// Other math functions (power, root, exp, log, etc...)?
//
//-----------------------------------------------------------------------------
// List of features to be implemented from Don Telian's list:
//
// Equation-based sources (voltage, current, numeric)
// AM:  This cannot be done due to language limitations.  Even
//      though string parameters do exist, they cannot be
//      interpreted as expression (i.e. converted to "code").
//      However, it would be easy to write equation based source
//      statements in the templates directly assuming some
//      familiarity with the language.  This could prevent or
//      make the substitution of native SPICE elements impossible,
//      though.
//
// B element with c_comp on/off and table scaling (static/dynamic)
//  - Hspice can’t do c_comp or dynamic
// AM:  Both of these features can be added easily to the existing
//      *-AMS IBIS buffer implementations.  Dynamic scaling would
//      require addition control ports, as discussed with the
//      controlled R, L, and C building blocks.  We need to decide
//      how many and what type of control ports we want to add.
//
// Expressions in parameters
// AM:  Same answer as 1st entry above.
//
// Table-based PWL sources
// AM:  This can be done using the $table_model keyword.  The
//      problem is that currently the table or the file name
//      reference to the table is hard coded, which prevents
//      mutiple instances of the same exact module name to be
//      used with different data tables.  There is some hope 
//      that in the LRM v2.3 this will be fixed.  Until then
//      manual cutting and pasting will be required.
//
//***   The PWL tables can be passed into the module as numerical
//***   arrays from the calling statement.  This would eliminate
//***   the need for cutting and pasting the module itself for
//***   multiple instantiation with different data.
//
// Time-controlled sources
//  - Hspice can not do
// AM:  These can be implemented in *-AMS easily.  We just need to
//      decide on what they should do.
//
// Use of previously computed value (prev(x)) in equation
//  - Hspice can not
// AM:  Need to check into this how it can be done in *-AMS.
//
//=============================================================================
//=============================================================================
module IBIS_R (p, n);
  electrical p, n;
  branch    (p, n) Out;
  parameter real Rval  = 1.0;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * Rval * I(Out);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCR (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * V(In) * I(Out);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCR (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real Scale = 1.0;

  analog begin
    V(In)  <+ 0.0;
    V(Out) <+ Scale * I(In) * I(Out);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_C (p, n);
  electrical p, n;
  branch    (p, n) Out;
  parameter real Cval  = 1.0;
  parameter real V0    = 0.0;
  parameter real Scale = 1.0;

  analog begin
    if (Scale * Cval > 0.0) begin
      V(Out) <+ idt(I(Out), V0 * Scale * Cval) / (Scale * Cval);
    end else begin
      I(Out) <+ 0.0;
    end
  end
endmodule
//=============================================================================
//=============================================================================
//module IBIS_VCC (p, n, ps, ns);
//  electrical p, n, ps, ns;
//  branch    (ps, ns) In;
//  branch    (p,  n)  Out;
//  parameter real V0    = 0.0;
//  parameter real Scale = 1.0;
//
//  analog begin
//    if (Scale * V(In) > 0.0) begin
//      V(Out) <+ idt(I(Out), V0 * Scale * V(In)) / (Scale * V(In));
//   end else begin
//      I(Out) <+ 0.0;
//    end
//  end
//endmodule
//=============================================================================
//=============================================================================
//module IBIS_CCC (p, n, ps, ns);
//  electrical p, n, ps, ns;
//  branch    (ps, ns) In;
//  branch    (p,  n)  Out;
//  parameter real V0    = 0.0;
//  parameter real Scale = 1.0;
//
//  analog begin
//    V(In)  <+ 0.0;
//    if (Scale * I(In) > 0.0) begin
//      V(Out) <+ idt(I(Out), V0 * Scale * I(In)) / (Scale * I(In));
//   end else begin
//      I(Out) <+ 0.0;
//    end
//  end
//endmodule
//=============================================================================
//=============================================================================
module IBIS_L (p, n);
  electrical p, n;
  branch    (p, n) Out;
  parameter real Lval  = 1.0;
  parameter real I0    = 0.0;
  parameter real Scale = 1.0;

  analog begin
    if (Scale * Lval > 0.0) begin
      I(Out) <+ idt(V(Out), I0 * Scale * Lval) / (Scale * Lval);
    end else begin
      V(Out) <+ 0.0;
    end
  end
endmodule
//=============================================================================
//=============================================================================
//module IBIS_VCL (p, n, ps, ns);
//  electrical p, n, ps, ns;
//  branch    (ps, ns) In;
//  branch    (p,  n)  Out;
//  parameter real I0    = 0.0;
//  parameter real Scale = 1.0;
//
//  analog begin
//    if (Scale * V(In) > 0.0) begin
//      I(Out) <+ idt(V(Out), I0 * Scale * V(In)) / (Scale * V(In));
//    end else begin
//      V(Out) <+ 0.0;
//    end
//  end
//endmodule
//=============================================================================
//=============================================================================
//module IBIS_CCL (p, n, ps, ns);
//  electrical p, n, ps, ns;
//  branch    (ps, ns) In;
//  branch    (p,  n)  Out;
//  parameter real I0    = 0.0;
//  parameter real Scale = 1.0;
//
//  analog begin
//    V(In)  <+ 0.0;
//    if (Scale * I(In) > 0.0) begin
//      I(Out) <+ idt(V(Out), I0 * Scale * I(In)) / (Scale * I(In));
//    end else begin
//      V(Out) <+ 0.0;
//    end
//  end
//endmodule
//=============================================================================
//=============================================================================
module IBIS_K (p1, n1, p2, n2);
  electrical p1, n1, p2, n2;
  branch    (p1, n1) Out1;
  branch    (p2, n2) Out2;
  parameter real Lval_1 = 1.0;
  parameter real Lval_2 = 1.0;
  parameter real Kval   = 0.0;
  parameter real I0_1   = 0.0;
  parameter real I0_2   = 0.0;
  parameter real Scale  = 1.0;

  // This should really be a "localparam" to prevent it from being
  // changed from the outside.
  parameter real M = Kval * sqrt(Lval_1 * Lval_2);

  analog begin
    if (Scale * Lval_1 > 0.0) begin
      I(Out1) <+ (idt(V(Out1), I0_1*Scale*Lval_1) - M*(I(Out2)-I0_2)) / (Scale*Lval_1);
    end else begin
      V(Out1) <+ 0.0;
    end
    if (Scale * Lval_2 > 0.0) begin
      I(Out2) <+ (idt(V(Out2), I0_2*Scale*Lval_2) - M*(I(Out1)-I0_1)) / (Scale*Lval_2);
    end else begin
      V(Out2) <+ 0.0;
    end
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_V (p, n);
  electrical p, n;
  branch    (p, n) Out;
  parameter real Vdc   = 1.0;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * Vdc;
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCVS (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * V(In);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCVS (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real Scale = 1.0;

  analog begin
    V(In)  <+ 0.0;
    V(Out) <+ Scale * I(In);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCVS_DELAY (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real TD    = 0.0;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * absdelay(V(In), TD);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCVS_DELAY (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real TD    = 0.0;
  parameter real Scale = 1.0;

  analog begin
    V(In)  <+ 0.0;
    V(Out) <+ Scale * absdelay(I(In), TD);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCVS_SUM (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * (V(In1) + V(In2));
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCVS_SUM (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(In1) <+ 0.0;
    V(In2) <+ 0.0;
    V(Out) <+ Scale * (I(In1) + I(In2));
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCVS_MULT (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(Out) <+ Scale * V(In1) * V(In2);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCVS_MULT (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(In1) <+ 0.0;
    V(In2) <+ 0.0;
    V(Out) <+ Scale * I(In1) * I(In2);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCVS_DIV (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    if (V(In2) != 0.0) begin
      V(Out) <+ Scale * V(In1) / V(In2);
    end else begin
      V(Out) <+ 1.0e+15;          // instead of inf;
    end
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCVS_DIV (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(In1) <+ 0.0;
    V(In2) <+ 0.0;
    if (I(In2) != 0.0) begin
      V(Out) <+ Scale * I(In1) / I(In2);
    end else begin
      V(Out) <+ 1.0e+15;          // instead of inf;
    end
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_I (p, n);
  electrical p, n;
  branch    (p, n) Out;
  parameter real Idc   = 1.0;
  parameter real Scale = 1.0;

  analog begin
    I(Out) <+ Scale * Idc;
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCCS (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real Scale = 1.0;

  analog begin
    I(Out) <+ Scale * V(In);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCCS (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real Scale = 1.0;

  analog begin
    V(In)  <+ 0.0;
    I(Out) <+ Scale * I(In);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCCS_DELAY (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real TD    = 0.0;
  parameter real Scale = 1.0;

  analog begin
    I(Out) <+ Scale * absdelay(V(In), TD);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCCS_DELAY (p, n, ps, ns);
  electrical p, n, ps, ns;
  branch    (ps, ns) In;
  branch    (p,  n)  Out;
  parameter real TD    = 0.0;
  parameter real Scale = 1.0;

  analog begin
    V(In)  <+ 0.0;
    I(Out) <+ Scale * absdelay(I(In), TD);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCCS_SUM (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    I(Out) <+ Scale * (V(In1) + V(In2));
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCCS_SUM (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(In1) <+ 0.0;
    V(In2) <+ 0.0;
    I(Out) <+ Scale * (I(In1) + I(In2));
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCCS_MULT (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    I(Out) <+ Scale * V(In1) * V(In2);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCCS_MULT (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(In1) <+ 0.0;
    V(In2) <+ 0.0;
    I(Out) <+ Scale * I(In1) * I(In2);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_VCCS_DIV (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    if (V(In2) != 0.0) begin
      I(Out) <+ Scale * V(In1) / V(In2);
    end else begin
      I(Out) <+ 1.0e+15;          // instead of inf;
    end
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_CCCS_DIV (p, n, ps1, ns1, ps2, ns2);
  electrical p, n, ps1, ns1, ps2, ns2;
  branch    (ps1, ns1) In1;
  branch    (ps2, ns2) In2;
  branch    (p,   n)   Out;
  parameter real Scale = 1.0;

  analog begin
    V(In1) <+ 0.0;
    V(In2) <+ 0.0;
    if (I(In2) != 0.0) begin
      I(Out) <+ Scale * I(In1) / I(In2);
    end else begin
      I(Out) <+ 1.0e+15;          // instead of inf;
    end
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_T (n1, ref1, n2, ref2);
  electrical n1, ref1, n2, ref2;
  branch    (n1, ref1) End1;
  branch    (n2, ref2) End2;
  parameter real Z0 = 50.0;
  parameter real TD =  1.0e-9;

  analog begin
    V(End1) <+ absdelay(V(End2) + Z0*I(End2), TD) + Z0*I(End1);
    V(End2) <+ absdelay(V(End1) + Z0*I(End1), TD) + Z0*I(End2);
  end
endmodule
//=============================================================================
//=============================================================================
module IBIS_Input (PC_ref, GC_ref, Input, Rcv_D);
  output     Rcv_D;
  electrical Rcv_D;
  inout      Input, PC_ref, GC_ref;
  electrical Input, PC_ref, GC_ref;

  branch (PC_ref, Input) pc;      // branches driven in model
  branch (Input, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.50,      // C_comp splitting coefficients
    k_C_comp_gc = 0.50,

    Vinh        = 2.0,       // Default Vinh value
    Vinl        = 0.8;       // Default Vinl value
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //===========================================================================
  analog begin

    //=========================================================================
    // This section contains the analog expressions which calculate
    // the clamping currents from the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+ -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(gc) <+  $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
    // A simple receiver logic
    //-------------------------------------------------------------------------
    if      (V(gc) > Vinh)   V(Rcv_D) <+ 1.0;
    else if (V(gc) < Vinl)   V(Rcv_D) <+ 0.0;
    else                     V(Rcv_D) <+ 0.5;
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_Input
//=============================================================================
//=============================================================================
module IBIS_Output (PU_ref, PD_ref, Output, In_D, PC_ref, GC_ref);
  input      In_D;
  electrical In_D;
  inout      Output, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical Output, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, Output) pc;     // branches driven in model
  branch (PU_ref, Output) pu;
  branch (Output, PD_ref) pd;
  branch (Output, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pu_ref    = 5.0,       // Pullup reference voltage
    V_pd_ref    = 0.0,       // Pulldown reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpu_length = 4,
    IVpd_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pu[1:IVpu_length] = { 0.10,  0.00, -0.10, -0.20},
    V_pu[1:IVpu_length] = {-5.00,  0.00,  5.00, 10.00},
    I_pd[1:IVpd_length] = {-0.10,  0.00,  0.10,  0.20},
    V_pd[1:IVpd_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTr2_length = 4,
    VTf1_length = 4,
    VTf2_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {0.00,  0.00,     2.50,     2.50},
    Tr1[1:VTr1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vr2[1:VTr2_length] = {2.50,  2.50,     5.00,     5.00},
    Tr2[1:VTr2_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {5.00,  5.00,     2.50,     2.50},
    Tf1[1:VTf1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vf2[1:VTf2_length] = {2.50,  2.50,     0.00,     0.00},
    Tf2[1:VTf2_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 0.0,                 // V_fixture values
    Vfx_r2 = 5.0,
    Vfx_f1 = 5.0,
    Vfx_f2 = 0.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_r2 = 50.0,
    Rfx_f1 = 50.0,
    Rfx_f2 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kr2[1:Max_length];
  real Kf1[1:Max_length];
  real Kf2[1:Max_length];

  integer pu_on;                  // Logic state variables
  integer pu_off;
  integer pd_on;
  integer pd_off;

  real Tpu_on_event;              // Variables to store event time
  real Tpu_off_event;
  real Tpd_on_event;
  real Tpd_off_event;

  real Kpu;                       // Output current scaling coefficients
  real Kpd;
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_r2;
  integer  index_f1;
  integer  index_f2;
  integer  new_index_r1;
  integer  new_index_r2;
  integer  new_index_f1;
  integer  new_index_f2;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vr2_common[1:3];
  real     Vf1_common[1:3];
  real     Vf2_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_r2;
  real     dVwfm_f1;
  real     dVwfm_f2;

  real     Iout1;
  real     Iout2;
  real     Ipd1;
  real     Ipd2;
  real     Ipu1;
  real     Ipu2;
  real     Igc1;
  real     Igc2;
  real     Ipc1;
  real     Ipc2;
  real     den;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_r2      = 1;
      index_f1      = 1;
      index_f2      = 1;
      new_index_r1  = VTr1_length;
      new_index_r2  = VTr2_length;
      new_index_f1  = VTf1_length;
      new_index_f2  = VTf2_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tr2[1]);
      old_t = min(old_t,  Tf1[1]);
      old_t = min(old_t,  Tf2[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tr2[VTr2_length]);
      last_t = max(new_t, Tf1[VTf1_length]);
      last_t = max(new_t, Tf2[VTf2_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tr2[index_r2] <= old_t) new_index_r2 = index_r2 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        if (Tf2[index_f2] <= old_t) new_index_f2 = index_f2 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_r2 <= VTr2_length) begin
          if (Tr2[new_index_r2] < new_t) new_t = Tr2[new_index_r2];
          index_r2 = new_index_r2;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        if (new_index_f2 <= VTf2_length) begin
          if (Tf2[new_index_f2] < new_t) new_t = Tf2[new_index_f2];
          index_f2 = new_index_f2;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vr2_common[2] = $table_model(T_common[1], Tr2, Vr2, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");
      Vf2_common[2] = $table_model(T_common[1], Tf2, Vf2, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vr2_common[3] = $table_model(T_common[i+1], Tr2, Vr2, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
          Vf2_common[3] = $table_model(T_common[i+1], Tf2, Vf2, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tR2  = %e\tF1  = %e\tF2  = %e", T_common[i], Vr1_common[2], Vr2_common[2], Vf1_common[2], Vf2_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_r2 = 0.0;
          dVwfm_f1 = 0.0;
          dVwfm_f2 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_r2 = ((Vr2_common[2] - Vr2_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr2_common[3] - Vr2_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f2 = ((Vf2_common[2] - Vf2_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf2_common[3] - Vf2_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdR2  = %e\tdF1  = %e\tdF2 = %e", dVwfm_r1, dVwfm_r2, dVwfm_f1, dVwfm_f2);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 Kr2 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;
        Iout2 = ((Vr2_common[2] - Vfx_r2) / Rfx_r2) + C_comp * dVwfm_r2;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipc2  = -1.0 * $table_model(V_pc_ref - Vr2_common[2], V_pc, I_pc, "LL");
        Ipu1  = -1.0 * $table_model(V_pu_ref - Vr1_common[2], V_pu, I_pu, "LL");
        Ipu2  = -1.0 * $table_model(V_pu_ref - Vr2_common[2], V_pu, I_pu, "LL");
        Ipd1  =        $table_model(Vr1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Ipd2  =        $table_model(Vr2_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");
        Igc2  =        $table_model(Vr2_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        den   = (Ipu1 * Ipd2) - (Ipd1 * Ipu2);

        // Since these are rising waveforms
        // Kr1 == Kpu_on and Kr2 == Kpd_off
        Kr1[i] = (Ipd2*(Iout1+Igc1-Ipc1) + Ipd1*(Ipc2-Igc2-Iout2)) / den;
        Kr2[i] = (Ipu2*(Iout1+Igc1-Ipc1) + Ipu1*(Ipc2-Igc2-Iout2)) / den;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 Kf2 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;
        Iout2 = ((Vf2_common[2] - Vfx_f2) / Rfx_f2) + C_comp * dVwfm_f2;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipc2  = -1.0 * $table_model(V_pc_ref - Vf2_common[2], V_pc, I_pc, "LL");
        Ipu1  = -1.0 * $table_model(V_pu_ref - Vf1_common[2], V_pu, I_pu, "LL");
        Ipu2  = -1.0 * $table_model(V_pu_ref - Vf2_common[2], V_pu, I_pu, "LL");
        Ipd1  =        $table_model(Vf1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Ipd2  =        $table_model(Vf2_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");
        Igc2  =        $table_model(Vf2_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        den   = (Ipu1 * Ipd2) - (Ipd1 * Ipu2);

        // Since these are falling waveforms
        // Kf1 == Kpd_on and Kf2 == Kpu_off
        Kf1[i] = (Ipu2*(Iout1+Igc1-Ipc1) + Ipu1*(Ipc2-Igc2-Iout2)) / den;
        Kf2[i] = (Ipd2*(Iout1+Igc1-Ipc1) + Ipd1*(Ipc2-Igc2-Iout2)) / den;

        //$strobe("\nTcom = %e\tKr1 = %e\tKr2 = %e\tKf1 = %e\tKf2 = %e", T_common[i], Kr1[i], Kr2[i], Kf1[i], Kf2[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vr2_common[1] = Vr2_common[2];
        Vf1_common[1] = Vf1_common[2];
        Vf2_common[1] = Vf2_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vr2_common[2] = Vr2_common[3];
        Vf1_common[2] = Vf1_common[3];
        Vf2_common[2] = Vf2_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kr2[i] = Kr2[Common_length];
        Kf1[i] = Kf1[Common_length];
        Kf2[i] = Kf2[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if (V(In_D) > Vth) begin         // Find initial logic state
        pu_on  = 1;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 0;
      end else if (V(In_D) < Vth) begin
        pu_on  = 0;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 1;
      end else begin
        pu_on  = 1;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 0;
      end
      //=======================================================================
    end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5)) begin
      if (V(In_D) > Vth) begin         // Find logic state
        pu_on  = 1;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 0;
      end else if (V(In_D) < Vth) begin
        pu_on  = 0;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 1;
      end else begin
        pu_on  = 1;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 0;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D changed
    //-------------------------------------------------------------------------
    @(cross(pu_on-0.5,  0)) begin
      Tpu_on_event  = $abstime;
      //$strobe("\nPU_ON  = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pu_on, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pu_off-0.5, 0)) begin
      Tpu_off_event = $abstime;
      //$strobe("\nPU_OFF = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pu_off, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_on-0.5,  0)) begin
      Tpd_on_event  = $abstime;
      //$strobe("\nPD_ON  = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pd_on, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_off-0.5, 0)) begin
      Tpd_off_event = $abstime;
      //$strobe("\nPD_OFF = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pd_off, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpu_on_event == 0.0 && Tpu_off_event == 0.0 &&
        Tpd_on_event == 0.0 && Tpd_off_event == 0.0) begin
                                                 // Initialization
      if (pu_on == 1) begin                      // Start with the end of the
        Kpu = Kr1[Common_length];                // Vt curves for those which
      end else if (pu_off == 1) begin            // are fully on initially
        Kpu = Kf2[Common_length];
      end else begin
        Kpu = 0.0;
      end

      if (pd_on == 1) begin
        Kpd = Kf1[Common_length];
      end else if (pd_off == 1) begin
        Kpd = Kr2[Common_length];
      end else begin
        Kpd = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pu_on  == 1) Kpu = $table_model($abstime - Tpu_on_event,  T_common, Kr1, "LL");
      else if (pu_off == 1) Kpu = $table_model($abstime - Tpu_off_event, T_common, Kf2, "LL");
      else                  Kpu = Kr1[1];

      if      (pd_on  == 1) Kpd = $table_model($abstime - Tpd_on_event,  T_common, Kf1, "LL");
      else if (pd_off == 1) Kpd = $table_model($abstime - Tpd_off_event, T_common, Kr2, "LL");
      else                  Kpd = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+ Kpu * -$table_model(V(pu), V_pu, I_pu, "LL") + k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+ Kpd *  $table_model(V(pd), V_pd, I_pd, "LL") + k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_Output
//=============================================================================
//=============================================================================
module IBIS_IO (PU_ref, PD_ref, IO, In_D, En_D, Rcv_D, PC_ref, GC_ref);
  input      In_D, En_D;
  electrical In_D, En_D;
  output     Rcv_D;
  electrical Rcv_D;
  inout      IO, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical IO, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, IO) pc;    // branches driven in model
  branch (PU_ref, IO) pu;
  branch (IO, PD_ref) pd;
  branch (IO, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    Vinh        = 2.0,       // Default Vinh value
    Vinl        = 0.8,       // Default Vinl value

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pu_ref    = 5.0,       // Pullup reference voltage
    V_pd_ref    = 0.0,       // Pulldown reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpu_length = 4,
    IVpd_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pu[1:IVpu_length] = { 0.10,  0.00, -0.10, -0.20},
    V_pu[1:IVpu_length] = {-5.00,  0.00,  5.00, 10.00},
    I_pd[1:IVpd_length] = {-0.10,  0.00,  0.10,  0.20},
    V_pd[1:IVpd_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTr2_length = 4,
    VTf1_length = 4,
    VTf2_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {0.00,  0.00,     2.50,     2.50},
    Tr1[1:VTr1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vr2[1:VTr2_length] = {2.50,  2.50,     5.00,     5.00},
    Tr2[1:VTr2_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {5.00,  5.00,     2.50,     2.50},
    Tf1[1:VTf1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vf2[1:VTf2_length] = {2.50,  2.50,     0.00,     0.00},
    Tf2[1:VTf2_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 0.0,                 // V_fixture values
    Vfx_r2 = 5.0,
    Vfx_f1 = 5.0,
    Vfx_f2 = 0.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_r2 = 50.0,
    Rfx_f1 = 50.0,
    Rfx_f2 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kr2[1:Max_length];
  real Kf1[1:Max_length];
  real Kf2[1:Max_length];

  integer pu_on;                  // Logic state variables
  integer pu_off;
  integer pd_on;
  integer pd_off;

  real Tpu_on_event;              // Variables to store event time
  real Tpu_off_event;
  real Tpd_on_event;
  real Tpd_off_event;

  real Kpu;                       // Output current scaling coefficients
  real Kpd;
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_r2;
  integer  index_f1;
  integer  index_f2;
  integer  new_index_r1;
  integer  new_index_r2;
  integer  new_index_f1;
  integer  new_index_f2;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vr2_common[1:3];
  real     Vf1_common[1:3];
  real     Vf2_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_r2;
  real     dVwfm_f1;
  real     dVwfm_f2;

  real     Iout1;
  real     Iout2;
  real     Ipd1;
  real     Ipd2;
  real     Ipu1;
  real     Ipu2;
  real     Igc1;
  real     Igc2;
  real     Ipc1;
  real     Ipc2;
  real     den;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_r2      = 1;
      index_f1      = 1;
      index_f2      = 1;
      new_index_r1  = VTr1_length;
      new_index_r2  = VTr2_length;
      new_index_f1  = VTf1_length;
      new_index_f2  = VTf2_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tr2[1]);
      old_t = min(old_t,  Tf1[1]);
      old_t = min(old_t,  Tf2[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tr2[VTr2_length]);
      last_t = max(new_t, Tf1[VTf1_length]);
      last_t = max(new_t, Tf2[VTf2_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tr2[index_r2] <= old_t) new_index_r2 = index_r2 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        if (Tf2[index_f2] <= old_t) new_index_f2 = index_f2 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_r2 <= VTr2_length) begin
          if (Tr2[new_index_r2] < new_t) new_t = Tr2[new_index_r2];
          index_r2 = new_index_r2;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        if (new_index_f2 <= VTf2_length) begin
          if (Tf2[new_index_f2] < new_t) new_t = Tf2[new_index_f2];
          index_f2 = new_index_f2;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vr2_common[2] = $table_model(T_common[1], Tr2, Vr2, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");
      Vf2_common[2] = $table_model(T_common[1], Tf2, Vf2, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vr2_common[3] = $table_model(T_common[i+1], Tr2, Vr2, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
          Vf2_common[3] = $table_model(T_common[i+1], Tf2, Vf2, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tR2  = %e\tF1  = %e\tF2  = %e", T_common[i], Vr1_common[2], Vr2_common[2], Vf1_common[2], Vf2_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_r2 = 0.0;
          dVwfm_f1 = 0.0;
          dVwfm_f2 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_r2 = ((Vr2_common[2] - Vr2_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr2_common[3] - Vr2_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f2 = ((Vf2_common[2] - Vf2_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf2_common[3] - Vf2_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdR2  = %e\tdF1  = %e\tdF2 = %e", dVwfm_r1, dVwfm_r2, dVwfm_f1, dVwfm_f2);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 Kr2 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;
        Iout2 = ((Vr2_common[2] - Vfx_r2) / Rfx_r2) + C_comp * dVwfm_r2;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipc2  = -1.0 * $table_model(V_pc_ref - Vr2_common[2], V_pc, I_pc, "LL");
        Ipu1  = -1.0 * $table_model(V_pu_ref - Vr1_common[2], V_pu, I_pu, "LL");
        Ipu2  = -1.0 * $table_model(V_pu_ref - Vr2_common[2], V_pu, I_pu, "LL");
        Ipd1  =        $table_model(Vr1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Ipd2  =        $table_model(Vr2_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");
        Igc2  =        $table_model(Vr2_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        den   = (Ipu1 * Ipd2) - (Ipd1 * Ipu2);

        // Since these are rising waveforms
        // Kr1 == Kpu_on and Kr2 == Kpd_off
        Kr1[i] = (Ipd2*(Iout1+Igc1-Ipc1) + Ipd1*(Ipc2-Igc2-Iout2)) / den;
        Kr2[i] = (Ipu2*(Iout1+Igc1-Ipc1) + Ipu1*(Ipc2-Igc2-Iout2)) / den;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 Kf2 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;
        Iout2 = ((Vf2_common[2] - Vfx_f2) / Rfx_f2) + C_comp * dVwfm_f2;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipc2  = -1.0 * $table_model(V_pc_ref - Vf2_common[2], V_pc, I_pc, "LL");
        Ipu1  = -1.0 * $table_model(V_pu_ref - Vf1_common[2], V_pu, I_pu, "LL");
        Ipu2  = -1.0 * $table_model(V_pu_ref - Vf2_common[2], V_pu, I_pu, "LL");
        Ipd1  =        $table_model(Vf1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Ipd2  =        $table_model(Vf2_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");
        Igc2  =        $table_model(Vf2_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        den   = (Ipu1 * Ipd2) - (Ipd1 * Ipu2);

        // Since these are falling waveforms
        // Kf1 == Kpd_on and Kf2 == Kpu_off
        Kf1[i] = (Ipu2*(Iout1+Igc1-Ipc1) + Ipu1*(Ipc2-Igc2-Iout2)) / den;
        Kf2[i] = (Ipd2*(Iout1+Igc1-Ipc1) + Ipd1*(Ipc2-Igc2-Iout2)) / den;

        //$strobe("\nTcom = %e\tKr1 = %e\tKr2 = %e\tKf1 = %e\tKf2 = %e", T_common[i], Kr1[i], Kr2[i], Kf1[i], Kf2[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vr2_common[1] = Vr2_common[2];
        Vf1_common[1] = Vf1_common[2];
        Vf2_common[1] = Vf2_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vr2_common[2] = Vr2_common[3];
        Vf1_common[2] = Vf1_common[3];
        Vf2_common[2] = Vf2_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kr2[i] = Kr2[Common_length];
        Kf1[i] = Kf1[Common_length];
        Kf2[i] = Kf2[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin   // Find initial logic state
        pu_on  = 1;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 1;
      end else begin
        pu_on  = 1;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 0;
      end
      //=======================================================================
    end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D, En_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5) or cross(V(En_D)-Vth, 0, 100p, 0.5)) begin
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin        // Find logic state
        pu_on  = 1;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 1;
      end else begin
        pu_on  = 1;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 0;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D or En_D changed
    //-------------------------------------------------------------------------
    @(cross(pu_on-0.5,  0)) begin
      Tpu_on_event  = $abstime;
      //$strobe("\nPU_ON  = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pu_on, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pu_off-0.5, 0)) begin
      Tpu_off_event = $abstime;
      //$strobe("\nPU_OFF = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pu_off, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_on-0.5,  0)) begin
      Tpd_on_event  = $abstime;
      //$strobe("\nPD_ON  = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pd_on, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_off-0.5, 0)) begin
      Tpd_off_event = $abstime;
      //$strobe("\nPD_OFF = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pd_off, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpu_on_event == 0.0 && Tpu_off_event == 0.0 &&
        Tpd_on_event == 0.0 && Tpd_off_event == 0.0) begin
                                                 // Initialization
      if (pu_on == 1) begin                      // Start with the end of the
        Kpu = Kr1[Common_length];                // Vt curves for those which
      end else if (pu_off == 1) begin            // are fully on initially
        Kpu = Kf2[Common_length];
      end else begin
        Kpu = 0.0;
      end

      if (pd_on == 1) begin
        Kpd = Kf1[Common_length];
      end else if (pd_off == 1) begin
        Kpd = Kr2[Common_length];
      end else begin
        Kpd = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pu_on  == 1) Kpu = $table_model($abstime - Tpu_on_event,  T_common, Kr1, "LL");
      else if (pu_off == 1) Kpu = $table_model($abstime - Tpu_off_event, T_common, Kf2, "LL");
      else                  Kpu = Kr1[1];

      if      (pd_on  == 1) Kpd = $table_model($abstime - Tpd_on_event,  T_common, Kf1, "LL");
      else if (pd_off == 1) Kpd = $table_model($abstime - Tpd_off_event, T_common, Kr2, "LL");
      else                  Kpd = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+ Kpu * -$table_model(V(pu), V_pu, I_pu, "LL") + k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+ Kpd *  $table_model(V(pd), V_pd, I_pd, "LL") + k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
    // A simple receiver logic
    //-------------------------------------------------------------------------
    if      (V(pd) > Vinh)   V(Rcv_D) <+ 1.0;
    else if (V(pd) < Vinl)   V(Rcv_D) <+ 0.0;
    else                     V(Rcv_D) <+ 0.5;
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_IO
//=============================================================================
//=============================================================================
module IBIS_3state (PU_ref, PD_ref, IO, In_D, En_D, PC_ref, GC_ref);
  input      In_D, En_D;
  electrical In_D, En_D;
  inout      IO, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical IO, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, IO) pc;    // branches driven in model
  branch (PU_ref, IO) pu;
  branch (IO, PD_ref) pd;
  branch (IO, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    Vinh        = 2.0,       // Default Vinh value
    Vinl        = 0.8,       // Default Vinl value

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pu_ref    = 5.0,       // Pullup reference voltage
    V_pd_ref    = 0.0,       // Pulldown reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpu_length = 4,
    IVpd_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pu[1:IVpu_length] = { 0.10,  0.00, -0.10, -0.20},
    V_pu[1:IVpu_length] = {-5.00,  0.00,  5.00, 10.00},
    I_pd[1:IVpd_length] = {-0.10,  0.00,  0.10,  0.20},
    V_pd[1:IVpd_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTr2_length = 4,
    VTf1_length = 4,
    VTf2_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {0.00,  0.00,     2.50,     2.50},
    Tr1[1:VTr1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vr2[1:VTr2_length] = {2.50,  2.50,     5.00,     5.00},
    Tr2[1:VTr2_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {5.00,  5.00,     2.50,     2.50},
    Tf1[1:VTf1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vf2[1:VTf2_length] = {2.50,  2.50,     0.00,     0.00},
    Tf2[1:VTf2_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 0.0,                 // V_fixture values
    Vfx_r2 = 5.0,
    Vfx_f1 = 5.0,
    Vfx_f2 = 0.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_r2 = 50.0,
    Rfx_f1 = 50.0,
    Rfx_f2 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kr2[1:Max_length];
  real Kf1[1:Max_length];
  real Kf2[1:Max_length];

  integer pu_on;                  // Logic state variables
  integer pu_off;
  integer pd_on;
  integer pd_off;

  real Tpu_on_event;              // Variables to store event time
  real Tpu_off_event;
  real Tpd_on_event;
  real Tpd_off_event;

  real Kpu;                       // Output current scaling coefficients
  real Kpd;
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_r2;
  integer  index_f1;
  integer  index_f2;
  integer  new_index_r1;
  integer  new_index_r2;
  integer  new_index_f1;
  integer  new_index_f2;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vr2_common[1:3];
  real     Vf1_common[1:3];
  real     Vf2_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_r2;
  real     dVwfm_f1;
  real     dVwfm_f2;

  real     Iout1;
  real     Iout2;
  real     Ipd1;
  real     Ipd2;
  real     Ipu1;
  real     Ipu2;
  real     Igc1;
  real     Igc2;
  real     Ipc1;
  real     Ipc2;
  real     den;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_r2      = 1;
      index_f1      = 1;
      index_f2      = 1;
      new_index_r1  = VTr1_length;
      new_index_r2  = VTr2_length;
      new_index_f1  = VTf1_length;
      new_index_f2  = VTf2_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tr2[1]);
      old_t = min(old_t,  Tf1[1]);
      old_t = min(old_t,  Tf2[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tr2[VTr2_length]);
      last_t = max(new_t, Tf1[VTf1_length]);
      last_t = max(new_t, Tf2[VTf2_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tr2[index_r2] <= old_t) new_index_r2 = index_r2 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        if (Tf2[index_f2] <= old_t) new_index_f2 = index_f2 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_r2 <= VTr2_length) begin
          if (Tr2[new_index_r2] < new_t) new_t = Tr2[new_index_r2];
          index_r2 = new_index_r2;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        if (new_index_f2 <= VTf2_length) begin
          if (Tf2[new_index_f2] < new_t) new_t = Tf2[new_index_f2];
          index_f2 = new_index_f2;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vr2_common[2] = $table_model(T_common[1], Tr2, Vr2, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");
      Vf2_common[2] = $table_model(T_common[1], Tf2, Vf2, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vr2_common[3] = $table_model(T_common[i+1], Tr2, Vr2, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
          Vf2_common[3] = $table_model(T_common[i+1], Tf2, Vf2, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tR2  = %e\tF1  = %e\tF2  = %e", T_common[i], Vr1_common[2], Vr2_common[2], Vf1_common[2], Vf2_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_r2 = 0.0;
          dVwfm_f1 = 0.0;
          dVwfm_f2 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_r2 = ((Vr2_common[2] - Vr2_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr2_common[3] - Vr2_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f2 = ((Vf2_common[2] - Vf2_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf2_common[3] - Vf2_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdR2  = %e\tdF1  = %e\tdF2 = %e", dVwfm_r1, dVwfm_r2, dVwfm_f1, dVwfm_f2);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 Kr2 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;
        Iout2 = ((Vr2_common[2] - Vfx_r2) / Rfx_r2) + C_comp * dVwfm_r2;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipc2  = -1.0 * $table_model(V_pc_ref - Vr2_common[2], V_pc, I_pc, "LL");
        Ipu1  = -1.0 * $table_model(V_pu_ref - Vr1_common[2], V_pu, I_pu, "LL");
        Ipu2  = -1.0 * $table_model(V_pu_ref - Vr2_common[2], V_pu, I_pu, "LL");
        Ipd1  =        $table_model(Vr1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Ipd2  =        $table_model(Vr2_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");
        Igc2  =        $table_model(Vr2_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        den   = (Ipu1 * Ipd2) - (Ipd1 * Ipu2);

        // Since these are rising waveforms
        // Kr1 == Kpu_on and Kr2 == Kpd_off
        Kr1[i] = (Ipd2*(Iout1+Igc1-Ipc1) + Ipd1*(Ipc2-Igc2-Iout2)) / den;
        Kr2[i] = (Ipu2*(Iout1+Igc1-Ipc1) + Ipu1*(Ipc2-Igc2-Iout2)) / den;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 Kf2 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;
        Iout2 = ((Vf2_common[2] - Vfx_f2) / Rfx_f2) + C_comp * dVwfm_f2;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipc2  = -1.0 * $table_model(V_pc_ref - Vf2_common[2], V_pc, I_pc, "LL");
        Ipu1  = -1.0 * $table_model(V_pu_ref - Vf1_common[2], V_pu, I_pu, "LL");
        Ipu2  = -1.0 * $table_model(V_pu_ref - Vf2_common[2], V_pu, I_pu, "LL");
        Ipd1  =        $table_model(Vf1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Ipd2  =        $table_model(Vf2_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");
        Igc2  =        $table_model(Vf2_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        den   = (Ipu1 * Ipd2) - (Ipd1 * Ipu2);

        // Since these are falling waveforms
        // Kf1 == Kpd_on and Kf2 == Kpu_off
        Kf1[i] = (Ipu2*(Iout1+Igc1-Ipc1) + Ipu1*(Ipc2-Igc2-Iout2)) / den;
        Kf2[i] = (Ipd2*(Iout1+Igc1-Ipc1) + Ipd1*(Ipc2-Igc2-Iout2)) / den;

        //$strobe("\nTcom = %e\tKr1 = %e\tKr2 = %e\tKf1 = %e\tKf2 = %e", T_common[i], Kr1[i], Kr2[i], Kf1[i], Kf2[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vr2_common[1] = Vr2_common[2];
        Vf1_common[1] = Vf1_common[2];
        Vf2_common[1] = Vf2_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vr2_common[2] = Vr2_common[3];
        Vf1_common[2] = Vf1_common[3];
        Vf2_common[2] = Vf2_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kr2[i] = Kr2[Common_length];
        Kf1[i] = Kf1[Common_length];
        Kf2[i] = Kf2[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin   // Find initial logic state
        pu_on  = 1;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 1;
      end else begin
        pu_on  = 1;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 0;
      end
      //=======================================================================
    end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D, En_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5) or cross(V(En_D)-Vth, 0, 100p, 0.5)) begin
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin        // Find logic state
        pu_on  = 1;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 0;
        pd_off = 1;
        pd_on  = 0;
        pu_off = 1;
      end else begin
        pu_on  = 1;
        pd_off = 0;
        pd_on  = 1;
        pu_off = 0;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D or En_D changed
    //-------------------------------------------------------------------------
    @(cross(pu_on-0.5,  0)) begin
      Tpu_on_event  = $abstime;
      //$strobe("\nPU_ON  = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pu_on, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pu_off-0.5, 0)) begin
      Tpu_off_event = $abstime;
      //$strobe("\nPU_OFF = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pu_off, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_on-0.5,  0)) begin
      Tpd_on_event  = $abstime;
      //$strobe("\nPD_ON  = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pd_on, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_off-0.5, 0)) begin
      Tpd_off_event = $abstime;
      //$strobe("\nPD_OFF = %d\tTpu_on = %e\tTpu_off = %e\tTpd_on = %e\tTpd_off = %e", pd_off, Tpu_on_event, Tpu_off_event, Tpd_on_event, Tpd_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpu_on_event == 0.0 && Tpu_off_event == 0.0 &&
        Tpd_on_event == 0.0 && Tpd_off_event == 0.0) begin
                                                 // Initialization
      if (pu_on == 1) begin                      // Start with the end of the
        Kpu = Kr1[Common_length];                // Vt curves for those which
      end else if (pu_off == 1) begin            // are fully on initially
        Kpu = Kf2[Common_length];
      end else begin
        Kpu = 0.0;
      end

      if (pd_on == 1) begin
        Kpd = Kf1[Common_length];
      end else if (pd_off == 1) begin
        Kpd = Kr2[Common_length];
      end else begin
        Kpd = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pu_on  == 1) Kpu = $table_model($abstime - Tpu_on_event,  T_common, Kr1, "LL");
      else if (pu_off == 1) Kpu = $table_model($abstime - Tpu_off_event, T_common, Kf2, "LL");
      else                  Kpu = Kr1[1];

      if      (pd_on  == 1) Kpd = $table_model($abstime - Tpd_on_event,  T_common, Kf1, "LL");
      else if (pd_off == 1) Kpd = $table_model($abstime - Tpd_off_event, T_common, Kr2, "LL");
      else                  Kpd = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+ Kpu * -$table_model(V(pu), V_pu, I_pu, "LL") + k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+ Kpd *  $table_model(V(pd), V_pd, I_pd, "LL") + k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_3state
//=============================================================================
//=============================================================================
module IBIS_OpenSink (PU_ref, PD_ref, IO, In_D, En_D, PC_ref, GC_ref);
  input      In_D, En_D;
  electrical In_D, En_D;
  inout      IO, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical IO, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, IO) pc;    // branches driven in model
  branch (PU_ref, IO) pu;
  branch (IO, PD_ref) pd;
  branch (IO, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pd_ref    = 0.0,       // Pulldown reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpd_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pd[1:IVpd_length] = {-0.10,  0.00,  0.10,  0.20},
    V_pd[1:IVpd_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTf1_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {2.50,  2.50,     5.00,     5.00},
    Tr1[1:VTr1_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {5.00,  5.00,     2.50,     2.50},
    Tf1[1:VTf1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 5.0,                 // V_fixture values
    Vfx_f1 = 5.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_f1 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kf1[1:Max_length];

  integer pd_on;                  // Logic state variables
  integer pd_off;

  real Tpd_on_event;              // Variables to store event time
  real Tpd_off_event;

  real Kpd;                       // Output current scaling coefficients
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_f1;
  integer  new_index_r1;
  integer  new_index_f1;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vf1_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_f1;

  real     Iout1;
  real     Ipd1;
  real     Igc1;
  real     Ipc1;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_f1      = 1;
      new_index_r1  = VTr1_length;
      new_index_f1  = VTf1_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tf1[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tf1[VTf1_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tF1  = %e", T_common[i], Vr1_common[2], Vf1_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_f1 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdF1  = %e", dVwfm_r1, dVwfm_f1);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipd1  =        $table_model(Vr1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a rising waveform Kr1 == Kpd_off
        if (Ipd1 == 0) Kr1[i] = 0.0;
        else           Kr1[i] = (Ipc1 - Igc1 - Iout1) / Ipd1;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipd1  =        $table_model(Vf1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a falling waveform Kf1 == Kpd_on
        if (Ipd1 == 0) Kf1[i] = 0.0;
        else           Kf1[i] = (Ipc1 - Igc1 - Iout1) / Ipd1;

        //$strobe("\nTcom = %e\tKr1 = %e\tKf1 = %e", T_common[i], Kr1[i], Kf1[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vf1_common[1] = Vf1_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vf1_common[2] = Vf1_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kf1[i] = Kf1[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin   // Find initial logic state
        pd_off = 1;
        pd_on  = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pd_off = 0;
        pd_on  = 1;
      end else if ((V(En_D) < Vth)) begin
        pd_off = 1;
        pd_on  = 0;
      end else begin
        pd_off = 0;
        pd_on  = 1;
      end
      //=======================================================================
  end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D, En_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5) or cross(V(En_D)-Vth, 0, 100p, 0.5)) begin
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin        // Find logic state
        pd_off = 1;
        pd_on  = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pd_off = 0;
        pd_on  = 1;
      end else if ((V(En_D) < Vth)) begin
        pd_off = 1;
        pd_on  = 0;
      end else begin
        pd_off = 0;
        pd_on  = 1;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D or En_D changed
    //-------------------------------------------------------------------------
    @(cross(pd_on-0.5,  0)) begin
      Tpd_on_event  = $abstime;
      //$strobe("\nPD_ON  = %d\tTpd_on = %e\tTpd_off = %e", pd_on, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_off-0.5, 0)) begin
      Tpd_off_event = $abstime;
      //$strobe("\nPD_OFF = %d\tTpd_on = %e\tTpd_off = %e", pd_off, Tpd_on_event, Tpd_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpd_on_event == 0.0 && Tpd_off_event == 0.0) begin
                                                 // Initialization
      if (pd_on == 1) begin                      // Start with the end of the
        Kpd = Kf1[Common_length];                // Vt curves for those which
      end else if (pd_off == 1) begin            // are fully on initially
        Kpd = Kr1[Common_length];
      end else begin
        Kpd = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pd_on  == 1) Kpd = $table_model($abstime - Tpd_on_event,  T_common, Kf1, "LL");
      else if (pd_off == 1) Kpd = $table_model($abstime - Tpd_off_event, T_common, Kr1, "LL");
      else                  Kpd = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+                                                k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+ Kpd *  $table_model(V(pd), V_pd, I_pd, "LL") + k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_OpenSink
//=============================================================================
//=============================================================================
module IBIS_IO_OpenSink (PU_ref, PD_ref, IO, In_D, En_D, Rcv_D, PC_ref, GC_ref);
  input      In_D, En_D;
  electrical In_D, En_D;
  output     Rcv_D;
  electrical Rcv_D;
  inout      IO, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical IO, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, IO) pc;    // branches driven in model
  branch (PU_ref, IO) pu;
  branch (IO, PD_ref) pd;
  branch (IO, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    Vinh        = 2.0,       // Default Vinh value
    Vinl        = 0.8,       // Default Vinl value

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pd_ref    = 0.0,       // Pulldown reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpd_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pd[1:IVpd_length] = {-0.10,  0.00,  0.10,  0.20},
    V_pd[1:IVpd_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTf1_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {2.50,  2.50,     5.00,     5.00},
    Tr1[1:VTr1_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {5.00,  5.00,     2.50,     2.50},
    Tf1[1:VTf1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 5.0,                 // V_fixture values
    Vfx_f1 = 5.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_f1 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kf1[1:Max_length];

  integer pd_on;                  // Logic state variables
  integer pd_off;

  real Tpd_on_event;              // Variables to store event time
  real Tpd_off_event;

  real Kpd;                       // Output current scaling coefficients
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_f1;
  integer  new_index_r1;
  integer  new_index_f1;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vf1_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_f1;

  real     Iout1;
  real     Ipd1;
  real     Igc1;
  real     Ipc1;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_f1      = 1;
      new_index_r1  = VTr1_length;
      new_index_f1  = VTf1_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tf1[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tf1[VTf1_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tF1  = %e", T_common[i], Vr1_common[2], Vf1_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_f1 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdF1  = %e", dVwfm_r1, dVwfm_f1);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipd1  =        $table_model(Vr1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a rising waveform Kr1 == Kpd_off
        if (Ipd1 == 0) Kr1[i] = 0.0;
        else           Kr1[i] = (Ipc1 - Igc1 - Iout1) / Ipd1;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipd1  =        $table_model(Vf1_common[2] - V_pd_ref, V_pd, I_pd, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a falling waveform Kf1 == Kpd_on
        if (Ipd1 == 0) Kf1[i] = 0.0;
        else           Kf1[i] = (Ipc1 - Igc1 - Iout1) / Ipd1;

        //$strobe("\nTcom = %e\tKr1 = %e\tKf1 = %e", T_common[i], Kr1[i], Kf1[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vf1_common[1] = Vf1_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vf1_common[2] = Vf1_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kf1[i] = Kf1[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin   // Find initial logic state
        pd_off = 1;
        pd_on  = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pd_off = 0;
        pd_on  = 1;
      end else if ((V(En_D) < Vth)) begin
        pd_off = 1;
        pd_on  = 0;
      end else begin
        pd_off = 0;
        pd_on  = 1;
      end
      //=======================================================================
  end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D, En_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5) or cross(V(En_D)-Vth, 0, 100p, 0.5)) begin
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin        // Find logic state
        pd_off = 1;
        pd_on  = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pd_off = 0;
        pd_on  = 1;
      end else if ((V(En_D) < Vth)) begin
        pd_off = 1;
        pd_on  = 0;
      end else begin
        pd_off = 0;
        pd_on  = 1;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D or En_D changed
    //-------------------------------------------------------------------------
    @(cross(pd_on-0.5,  0)) begin
      Tpd_on_event  = $abstime;
      //$strobe("\nPD_ON  = %d\tTpd_on = %e\tTpd_off = %e", pd_on, Tpd_on_event, Tpd_off_event);
    end
    @(cross(pd_off-0.5, 0)) begin
      Tpd_off_event = $abstime;
      //$strobe("\nPD_OFF = %d\tTpd_on = %e\tTpd_off = %e", pd_off, Tpd_on_event, Tpd_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpd_on_event == 0.0 && Tpd_off_event == 0.0) begin
                                                 // Initialization
      if (pd_on == 1) begin                      // Start with the end of the
        Kpd = Kf1[Common_length];                // Vt curves for those which
      end else if (pd_off == 1) begin            // are fully on initially
        Kpd = Kr1[Common_length];
      end else begin
        Kpd = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pd_on  == 1) Kpd = $table_model($abstime - Tpd_on_event,  T_common, Kf1, "LL");
      else if (pd_off == 1) Kpd = $table_model($abstime - Tpd_off_event, T_common, Kr1, "LL");
      else                  Kpd = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+                                                k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+ Kpd *  $table_model(V(pd), V_pd, I_pd, "LL") + k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
    // A simple receiver logic
    //-------------------------------------------------------------------------
    if      (V(pd) > Vinh)   V(Rcv_D) <+ 1.0;
    else if (V(pd) < Vinl)   V(Rcv_D) <+ 0.0;
    else                     V(Rcv_D) <+ 0.5;
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_IO_OpenSink
//=============================================================================
//=============================================================================
module IBIS_OpenSource (PU_ref, PD_ref, IO, In_D, En_D, PC_ref, GC_ref);
  input      In_D, En_D;
  electrical In_D, En_D;
  inout      IO, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical IO, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, IO) pc;    // branches driven in model
  branch (PU_ref, IO) pu;
  branch (IO, PD_ref) pd;
  branch (IO, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pu_ref    = 5.0,       // Pullup reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpu_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pu[1:IVpu_length] = { 0.10,  0.00, -0.10, -0.20},
    V_pu[1:IVpu_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTf1_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {0.00,  0.00,     2.50,     2.50},
    Tr1[1:VTr1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {2.50,  2.50,     0.00,     0.00},
    Tf1[1:VTf1_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 0.0,                 // V_fixture values
    Vfx_f1 = 0.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_f1 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kf1[1:Max_length];

  integer pu_on;                  // Logic state variables
  integer pu_off;

  real Tpu_on_event;              // Variables to store event time
  real Tpu_off_event;

  real Kpu;                       // Output current scaling coefficients
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_f1;
  integer  new_index_r1;
  integer  new_index_f1;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vf1_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_f1;

  real     Iout1;
  real     Ipu1;
  real     Igc1;
  real     Ipc1;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_f1      = 1;
      new_index_r1  = VTr1_length;
      new_index_f1  = VTf1_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tf1[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tf1[VTf1_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tF1  = %e", T_common[i], Vr1_common[2], Vf1_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_f1 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdF1  = %e", dVwfm_r1, dVwfm_f1);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipu1  =        $table_model(V_pu_ref - Vr1_common[2], V_pu, I_pu, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a rising waveform Kr1 == Kpu_on
        if (Ipu1 == 0) Kr1[i] = 0.0;
        else           Kr1[i] = (Ipc1 - Igc1 - Iout1) / Ipu1;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipu1  =        $table_model(V_pu_ref - Vf1_common[2], V_pu, I_pu, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a falling waveform Kf1 == Kpu_off
        if (Ipu1 == 0) Kf1[i] = 0.0;
        else           Kf1[i] = (Ipc1 - Igc1 - Iout1) / Ipu1;

        //$strobe("\nTcom = %e\tKr1 = %e\tKf1 = %e", T_common[i], Kr1[i], Kf1[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vf1_common[1] = Vf1_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vf1_common[2] = Vf1_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kf1[i] = Kf1[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin   // Find initial logic state
        pu_on  = 1;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 1;
        pu_off = 0;
      end else begin
        pu_on  = 0;
        pu_off = 1;
      end
      //=======================================================================
  end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D, En_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5) or cross(V(En_D)-Vth, 0, 100p, 0.5)) begin
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin        // Find logic state
        pu_on  = 1;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 1;
        pu_off = 0;
      end else begin
        pu_on  = 0;
        pu_off = 1;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D or En_D changed
    //-------------------------------------------------------------------------
    @(cross(pu_on-0.5,  0)) begin
      Tpu_on_event  = $abstime;
      //$strobe("\nPU_ON  = %d\tTpu_on = %e\tTpu_off = %e", pu_on, Tpu_on_event, Tpu_off_event);
    end
    @(cross(pu_off-0.5, 0)) begin
      Tpu_off_event = $abstime;
      //$strobe("\nPU_OFF = %d\tTpu_on = %e\tTpu_off = %e", pu_off, Tpu_on_event, Tpu_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpu_on_event == 0.0 && Tpu_off_event == 0.0) begin
                                                 // Initialization
      if (pu_on == 1) begin                      // Start with the end of the
        Kpu = Kr1[Common_length];                // Vt curves for those which
      end else if (pu_off == 1) begin            // are fully on initially
        Kpu = Kf1[Common_length];
      end else begin
        Kpu = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pu_on  == 1) Kpu = $table_model($abstime - Tpu_on_event,  T_common, Kr1, "LL");
      else if (pu_off == 1) Kpu = $table_model($abstime - Tpu_off_event, T_common, Kf1, "LL");
      else                  Kpu = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+ Kpu * -$table_model(V(pu), V_pu, I_pu, "LL") + k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+                                                k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_OpenSource
//=============================================================================
//=============================================================================
module IBIS_IO_OpenSource (PU_ref, PD_ref, IO, In_D, En_D, Rcv_D, PC_ref, GC_ref);
  input      In_D, En_D;
  electrical In_D, En_D;
  output     Rcv_D;
  electrical Rcv_D;
  inout      IO, PC_ref, PU_ref, PD_ref, GC_ref;
  electrical IO, PC_ref, PU_ref, PD_ref, GC_ref;

  branch (PC_ref, IO) pc;    // branches driven in model
  branch (PU_ref, IO) pu;
  branch (IO, PD_ref) pd;
  branch (IO, GC_ref) gc;

  parameter real
  //---------------------------------------------------------------------------
  // IBIS parameters
  //---------------------------------------------------------------------------
    C_comp      = 5.0p,      // Default C_comp value
    k_C_comp_pc = 0.25,      // C_comp splitting coefficients
    k_C_comp_pu = 0.25,
    k_C_comp_pd = 0.25,
    k_C_comp_gc = 0.25,

    Vinh        = 2.0,       // Default Vinh value
    Vinl        = 0.8,       // Default Vinl value

    V_pc_ref    = 5.0,       // Power clamp reference voltage
    V_pu_ref    = 5.0,       // Pullup reference voltage
    V_gc_ref    = 0.0;       // Ground clamp reference voltage
  //---------------------------------------------------------------------------
  // Vectors of the IV curve tables
  //---------------------------------------------------------------------------
  parameter integer
    IVpc_length = 4,
    IVpu_length = 4,
    IVgc_length = 4;

  parameter real
    I_pc[1:IVpc_length] = { 0.08,  0.00,  0.00,  0.00},
    V_pc[1:IVpc_length] = {-5.00, -1.00,  5.00, 10.00},
    I_pu[1:IVpu_length] = { 0.10,  0.00, -0.10, -0.20},
    V_pu[1:IVpu_length] = {-5.00,  0.00,  5.00, 10.00},
    I_gc[1:IVgc_length] = {-0.08,  0.00,  0.00,  0.00},
    V_gc[1:IVgc_length] = {-5.00, -1.00,  5.00, 10.00};
  //---------------------------------------------------------------------------
  // Vectors of the Vt curve tables
  //---------------------------------------------------------------------------
  parameter integer
    VTr1_length = 4,
    VTf1_length = 4;

  parameter real
    Vr1[1:VTr1_length] = {0.00,  0.00,     2.50,     2.50},
    Tr1[1:VTr1_length] = {0.00,  1.00e-9,  2.00e-9,  3.00e-9},
    Vf1[1:VTf1_length] = {2.50,  2.50,     0.00,     0.00},
    Tf1[1:VTf1_length] = {0.00,  0.50e-9,  0.80e-9,  3.00e-9};
  //---------------------------------------------------------------------------
  parameter real
    Vfx_r1 = 0.0,                 // V_fixture values
    Vfx_f1 = 0.0,

    Rfx_r1 = 50.0,                // R_fixture values
    Rfx_f1 = 50.0,
  //===========================================================================
  // Non-IBIS parameters
  //---------------------------------------------------------------------------
    Max_dt = 1.0e-12,             // Maximum dt between time points
    Vth    = 0.5;                 // Logic input threshold voltage
//-----------------------------------------------------------------------------
// Internal variables
//-----------------------------------------------------------------------------
// Calculate maximum array length for common time and K_*** tables
// and declare the arrays
//-----------------------------------------------------------------------------
//  parameter real    Max_t      = max(max(max(Tr1[VTr1_length], Tr2[VTr2_length]), Tf1[VTf1_length]), Tf2[VTf2_length]);
//  parameter real    Min_t      = min(min(min( Tr1[1], Tr2[1]), Tf1[1]), Tf2[1]);
//  parameter integer Max_length = ceil((Max_t-Min_t)/Max_dt) + VTr1_length + VTr2_length + VTf1_length + VTf2_length;

    parameter integer Max_length = 3016;   // Delete this when the above 3 lines work

  real T_common[1:Max_length];    // Common time axis array
  real Kr1[1:Max_length];         // Scaling coefficient arrays
  real Kf1[1:Max_length];

  integer pu_on;                  // Logic state variables
  integer pu_off;

  real Tpu_on_event;              // Variables to store event time
  real Tpu_off_event;

  real Kpu;                       // Output current scaling coefficients
//-----------------------------------------------------------------------------
// variables used in Common_time / Common_wfm / Kcoef "function"
//-----------------------------------------------------------------------------
  integer  i;
  integer  Common_length;
  integer  index_r1;
  integer  index_f1;
  integer  new_index_r1;
  integer  new_index_f1;
  real     old_t;
  real     new_t;
  real     last_t;
  real     additional_pts;
  real     new_dt;

  real     Vr1_common[1:3];       // Need 3 elements in order to do derivatives
  real     Vf1_common[1:3];

  real     dVwfm_r1;
  real     dVwfm_f1;

  real     Iout1;
  real     Ipu1;
  real     Igc1;
  real     Ipc1;
//=============================================================================
  analog begin

    @(initial_step) begin
      //=======================================================================
      // This "function" makes a common time axis for all Vt curves and
      // calculates the corresponding voltage points with linear interpolation
      // and/or linear extrapolation if necessary.  After that it converts the
      // Vt curves to K coefficients which will be used in the analog equations
      // to scale the IV curves.  The K coefficients are stored in arrays.
      //-----------------------------------------------------------------------
      // Initialize variables
      //-----------------------------------------------------------------------
      Common_length = 1;
      index_r1      = 1;
      index_f1      = 1;
      new_index_r1  = VTr1_length;
      new_index_f1  = VTf1_length;
      //-----------------------------------------------------------------------
      // Put the earliest time value of all given time vectors into "old_t"
      //-----------------------------------------------------------------------
      old_t = min(Tr1[1], Tf1[1]);
      //-----------------------------------------------------------------------
      // Put the latest time value of all given time vectors into "new_t"
      //-----------------------------------------------------------------------
      last_t = max(Tr1[VTr1_length], Tf1[VTf1_length]);
      new_t  = last_t;
      //-----------------------------------------------------------------------
      // Loop until latest time value is reached in each given time vector
      //-----------------------------------------------------------------------
      while (old_t < new_t) begin
        //---------------------------------------------------------------------
        // Save the time value from old_t in T_common
        // and increment Common_length counter
        //---------------------------------------------------------------------
        T_common[Common_length] = old_t;
        Common_length = Common_length + 1;
        //---------------------------------------------------------------------
        // Find which given time vector(s) have the lowest time value and
        // advance their temporary indexes
        //---------------------------------------------------------------------
        if (Tr1[index_r1] <= old_t) new_index_r1 = index_r1 + 1;
        if (Tf1[index_f1] <= old_t) new_index_f1 = index_f1 + 1;
        //---------------------------------------------------------------------
        // Find the lowest value at the new indexes in the given vector(s)
        // and update indexes for next iteration
        //---------------------------------------------------------------------
        if (new_index_r1 <= VTr1_length) begin
          if (Tr1[new_index_r1] < new_t) new_t = Tr1[new_index_r1];
          index_r1 = new_index_r1;
        end
        if (new_index_f1 <= VTf1_length) begin
          if (Tf1[new_index_f1] < new_t) new_t = Tf1[new_index_f1];
          index_f1 = new_index_f1;
        end
        //---------------------------------------------------------------------
        // Add points if dt is larger than Max_dt
        //---------------------------------------------------------------------
        if ((new_t-old_t) > Max_dt) begin
          additional_pts = ceil((new_t-old_t)/Max_dt);
          new_dt = (new_t-old_t) / additional_pts;
          for (i = 1; i < additional_pts; i = i + 1) begin
            //$strobe("\ndif = %e\tAdd = %d\tCL = %d\t T = %e\tnew_dt = ", (new_t-old_t), additional_pts, Common_length, T_common[Common_length-1]+new_dt, new_dt);
            T_common[Common_length] = T_common[Common_length-1] + new_dt;
            Common_length = Common_length + 1;
          end
        end
        //---------------------------------------------------------------------
        // Update variables for next iteration
        //---------------------------------------------------------------------
        old_t = new_t;
        new_t = last_t;
        //---------------------------------------------------------------------
      end  // while loop
      T_common[Common_length] = last_t;
      //=======================================================================

      //=======================================================================
      // Calculate the scaling coefficients for each time point along T_common
      //-----------------------------------------------------------------------
      // Store the first voltage point in the common waveform variables
      //-----------------------------------------------------------------------
      Vr1_common[2] = $table_model(T_common[1], Tr1, Vr1, "LL");
      Vf1_common[2] = $table_model(T_common[1], Tf1, Vf1, "LL");

      for (i = 1; i <= Common_length; i = i + 1) begin
      //-----------------------------------------------------------------------
      // Store the next point (for the derivative calculations)
      //-----------------------------------------------------------------------
        if (i < Common_length) begin
          Vr1_common[3] = $table_model(T_common[i+1], Tr1, Vr1, "LL");
          Vf1_common[3] = $table_model(T_common[i+1], Tf1, Vf1, "LL");
        end

        //$strobe("\nTcom = %e\tR1  = %e\tF1  = %e", T_common[i], Vr1_common[2], Vf1_common[2]);
        //---------------------------------------------------------------------
        // Calculate the derivative of each waveform for C_comp compensation
        //---------------------------------------------------------------------
        if ((i <= 1) || (i >= Common_length)) begin
          dVwfm_r1 = 0.0;         // Force the end point derivatives to zero
          dVwfm_f1 = 0.0;
        end else begin
          dVwfm_r1 = ((Vr1_common[2] - Vr1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vr1_common[3] - Vr1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
          dVwfm_f1 = ((Vf1_common[2] - Vf1_common[1]) / (T_common[i]   - T_common[i-1]) +
                      (Vf1_common[3] - Vf1_common[2]) / (T_common[i+1] - T_common[i]))  / 2.0;
        end

        //$strobe("dR1 = %e\tdF1  = %e", dVwfm_r1, dVwfm_f1);
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kr1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vr1_common[2] - Vfx_r1) / Rfx_r1) + C_comp * dVwfm_r1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vr1_common[2], V_pc, I_pc, "LL");
        Ipu1  =        $table_model(V_pu_ref - Vr1_common[2], V_pu, I_pu, "LL");
        Igc1  =        $table_model(Vr1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a rising waveform Kr1 == Kpu_on
        if (Ipu1 == 0) Kr1[i] = 0.0;
        else           Kr1[i] = (Ipc1 - Igc1 - Iout1) / Ipu1;
        //---------------------------------------------------------------------
        // Calculate intermediate variables and the Kf1 coefficients
        //---------------------------------------------------------------------
        Iout1 = ((Vf1_common[2] - Vfx_f1) / Rfx_f1) + C_comp * dVwfm_f1;

        Ipc1  = -1.0 * $table_model(V_pc_ref - Vf1_common[2], V_pc, I_pc, "LL");
        Ipu1  =        $table_model(V_pu_ref - Vf1_common[2], V_pu, I_pu, "LL");
        Igc1  =        $table_model(Vf1_common[2] - V_gc_ref, V_gc, I_gc, "LL");

        // Since these is a falling waveform Kf1 == Kpu_off
        if (Ipu1 == 0) Kf1[i] = 0.0;
        else           Kf1[i] = (Ipc1 - Igc1 - Iout1) / Ipu1;

        //$strobe("\nTcom = %e\tKr1 = %e\tKf1 = %e", T_common[i], Kr1[i], Kf1[i]);
        //---------------------------------------------------------------------
        // Shift points by one index for next iteration (for C_comp compensation)
        //---------------------------------------------------------------------
        Vr1_common[1] = Vr1_common[2];
        Vf1_common[1] = Vf1_common[2];

        Vr1_common[2] = Vr1_common[3];
        Vf1_common[2] = Vf1_common[3];

      end  // for loop
      //-----------------------------------------------------------------------
      // Fill the unused portion of the arrays with reasonable values
      // to prevent any misbehaving later with $table_model
      //=======================================================================
      for (i = Common_length + 1; i <= Max_length; i = i + 1) begin
        T_common[i] = T_common[i-1] + Max_dt;
        Kr1[i] = Kr1[Common_length];
        Kf1[i] = Kf1[Common_length];
      end  // for loop
      //-----------------------------------------------------------------------
      // End of Common_time / Common_wfm / Kcoef "function"
      //=======================================================================

      //=======================================================================
      // Initialize logic state variables
      //-----------------------------------------------------------------------
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin   // Find initial logic state
        pu_on  = 1;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 1;
        pu_off = 0;
      end else begin
        pu_on  = 0;
        pu_off = 1;
      end
      //=======================================================================
  end // @(initial_step)
   
    //=========================================================================
    // Catch: process to catch threshold crossings in In_D, En_D
    //-------------------------------------------------------------------------
    @(cross(V(In_D)-Vth, 0, 100p, 0.5) or cross(V(En_D)-Vth, 0, 100p, 0.5)) begin
      if ((V(En_D) > Vth) && (V(In_D) > Vth)) begin        // Find logic state
        pu_on  = 1;
        pu_off = 0;
      end else if ((V(En_D) > Vth) && (V(In_D) < Vth)) begin
        pu_on  = 0;
        pu_off = 1;
      end else if ((V(En_D) < Vth)) begin
        pu_on  = 1;
        pu_off = 0;
      end else begin
        pu_on  = 0;
        pu_off = 1;
      end
    end
    //-------------------------------------------------------------------------
    // Save time of event when In_D or En_D changed
    //-------------------------------------------------------------------------
    @(cross(pu_on-0.5,  0)) begin
      Tpu_on_event  = $abstime;
      //$strobe("\nPU_ON  = %d\tTpu_on = %e\tTpu_off = %e", pu_on, Tpu_on_event, Tpu_off_event);
    end
    @(cross(pu_off-0.5, 0)) begin
      Tpu_off_event = $abstime;
      //$strobe("\nPU_OFF = %d\tTpu_on = %e\tTpu_off = %e", pu_off, Tpu_on_event, Tpu_off_event);
    end
    //-------------------------------------------------------------------------
    // End of process Catch;
    //=========================================================================

    //=========================================================================
    // This section contains the simultaneous analog equations to find the
    // appropriate scaling coefficients according to the state the buffer.
    //-------------------------------------------------------------------------
    if (Tpu_on_event == 0.0 && Tpu_off_event == 0.0) begin
                                                 // Initialization
      if (pu_on == 1) begin                      // Start with the end of the
        Kpu = Kr1[Common_length];                // Vt curves for those which
      end else if (pu_off == 1) begin            // are fully on initially
        Kpu = Kf1[Common_length];
      end else begin
        Kpu = 0.0;
      end

    end else begin                // Look up coefficients in normal operation

      if      (pu_on  == 1) Kpu = $table_model($abstime - Tpu_on_event,  T_common, Kr1, "LL");
      else if (pu_off == 1) Kpu = $table_model($abstime - Tpu_off_event, T_common, Kf1, "LL");
      else                  Kpu = Kf1[1];

    end
    //=========================================================================
    // This section contains the analog expressions which calculate
    // the output currents of the IV curves and C_comp capacitors.
    //-------------------------------------------------------------------------
    I(pc) <+       -$table_model(V(pc), V_pc, I_pc, "LL") + k_C_comp_pc*C_comp*ddt(V(pc));
    I(pu) <+ Kpu * -$table_model(V(pu), V_pu, I_pu, "LL") + k_C_comp_pu*C_comp*ddt(V(pu));
    I(pd) <+                                                k_C_comp_pd*C_comp*ddt(V(pd));
    I(gc) <+        $table_model(V(gc), V_gc, I_gc, "LL") + k_C_comp_gc*C_comp*ddt(V(gc));
    //-------------------------------------------------------------------------
    // A simple receiver logic
    //-------------------------------------------------------------------------
    if      (V(pd) > Vinh)   V(Rcv_D) <+ 1.0;
    else if (V(pd) < Vinl)   V(Rcv_D) <+ 0.0;
    else                     V(Rcv_D) <+ 0.5;
    //-------------------------------------------------------------------------
  end
endmodule // IBIS_IO_OpenSource
//=============================================================================
//=============================================================================

////********************************************************************
//// Separate module for each N-coupled tline ...
////********************************************************************
//module ibis_buffer_in (pad, pc, gc, dout) ;
//  inout pad, pc, gc ;
//  output dout ;
//  electrical pad, pc, gc, dout ;
//  parameter string ibis_file ;
//  parameter string ibis_model ;
//  // TBD: implementation
//endmodule
////********************************************************************
//module ibis_buffer_out (pad, pu, pd, pc, gc, din, enable) ;
//  inout pad, pu, pd, pc, gc ;
//  input din, enable ;
//  electrical pad, pu, pd, pc, gc, din, enable ;
//  parameter string ibis_file ;
//  parameter string ibis_model ;
//  parameter real pu_scale ;
//  parameter real pd_scale ;
//  // TBD: implementation
//endmodule
////********************************************************************
//module ibis_buffer_io (pad, pu, pd, pc, gc, dout, din, enable) ;
//  inout pad, pu, pd, pc, gc ;
//  input din, enable ;
//  output dout ;
//  electrical pad, pu, pd, pc, gc, dout, din, enable ;
//  parameter string ibis_file ;
//  parameter string ibis_model ;
//  parameter real pu_scale ;
//  parameter real pd_scale ;
//  // TBD: implementation
//endmodule
////********************************************************************
////********************************************************************
//module ibis_buffer_series
//  inout pad1, pad2 ;
//  electrical pad1, pad2 ;
//  parameter string ibis_file ;
//  parameter string ibis_model ;
//  // TBD: implementation
//endmodule
////********************************************************************
//module ibis_sparam2 (a1, b1, ref1, a2, b2, ref2) ;
//  inout a1, b1, ref1, a2, b2, ref2 ;
//  electrical a1, b1, ref1, a2, b2, ref2 ;
//  parameter string tstone_file ;
//  // TBD: implementation
//endmodule
////********************************************************************