`include "disciplines.vams"
`include "constants.vams"

module myfet(d, g, s);
  inout electrical d, g, s;
  parameter real kp = 1m;
  parameter real vt = 1;
  real vgst;
  real cur;
  analog begin
   vgst = V(g,s)-vt;
   cur = 0.0;
   if (vgst > 0)
    if (vgst > V(d,s))
     cur = (vgst-0.5*V(d,s))*V(d,s);
    else
     cur = 0.5*pow(vgst, 2);
     I(d,s) <+ kp*cur;
   end
endmodule


module lccap(a,b);
  inout a,b; //Interface ports
  electrical a,b;
  electrical force_node; //Internal node

  real vext=0;
  real vv=0;
  real vforce=0;

  parameter real rforce=1;            // Rd
  parameter real cforce=0.001;        // Cd 
  parameter real tau = rforce*cforce; // tau=Rd*Cd

  real vcontrol= 0;  // RMS Voltage to control the angle 

  // capacitance 
  real cap = 0; 
  real qlc = 0;
  parameter real cmin=2.5e-15; // C  
  parameter real cmax=8.0e-15; // Cparallel

  real alpha = 0;
  parameter real delta_sq = 0.1;
  parameter real vtc = 2;
  parameter real vmc = 0.1;
 
  // transmittance / reflectance 
  real beta = 0; 
  real trans=0;
  
  parameter real tmin = 0.01;
  parameter real eta_sq = 0.1;
  parameter real vmo = 0.9;
  parameter real vto = 2;

  analog
  begin
    @(initial_step)
    begin
      //Initial voltage for internal node
      vforce = 0;
    end

    begin
      //Probing terminal voltage
      vext = V(a,b);
      vv = vext * vext;

      // Calculation of Internal node voltage
      I(force_node) <+ ddt(cforce * V(force_node)); // current into Cd 
      I(force_node) <+ (V(force_node)-vv) / rforce; // current from Rd     
      vforce = V(force_node);  
      vcontrol = sqrt(vforce); 			    // calculate RMS Voltage

      // C-V calculation
      alpha = (vcontrol - vtc) / vmc;
      cap = cmin + (2/`M_PI) * (cmax-cmin) * atan( (alpha + sqrt(alpha * alpha + delta_sq ))/2);
      qlc = cap * vext;   // delta Q
      I(a,b) <+ ddt(qlc); // dQ / dt

      // T-V calculation
      beta = (vcontrol - vto)/ vmo;
      trans = 1 - (1-tmin)*tanh( (beta+sqrt(beta*beta + eta_sq)) / 2 );

    end
  end
endmodule


module demo;
  electrical src_out, gate, out, gnd;
  ground gnd;

  analog begin
    V(out) <+ idt(6.0 * 100, 0, 0); // ramp 0-6V in 10ms
    $bound_step(10u/100);           // bessere Aufloesung
  end
  
  lccap lcc(out, gnd);
endmodule