//////////////////////////////////////////////////
// VerilogA model for    
//
// A Computationally Efficient Verilog-A ReRAM Model
//
// Nano Research Group, Electronics and Computer Science Department, 
// University of Southampton
// Department of Physics, Aristotle University of Thessaloniki
//
// April 2017
//
//////////////////////////////////////////////////


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

module ReRAM_model(p, n);
	inout p, n;

electrical p, n;

//'Switching sensitivity' parameters for v>0
parameter real Ap = 0.1675;

//'Switching sensitivity' parameters for v<0
parameter real An = -0.5633;

//Function 'r' parameters for v>0
parameter real a0 = 2738.9570;
parameter real a1 = 1969.2233; 

//Function 'r' parameters for v<0
parameter real b0 = 6462.8692;
parameter real b1 = 658.7840;  

//Initial resistive state 
parameter real Rinit = 5400; 

//I-V relationship parameters
parameter real ap=0.237;
parameter real an=0.247;
parameter real bp=3;
parameter real bn=3;

//Stp function parameters
parameter real c1=1e-3;     

//set-reset direction
parameter real s=1;    

real Rmp; 		// Function 'r' for v>0
real Rmn; 		// Function 'r' for v<0
real stpPOSv; 	// step function for v>0
real stpNEGv; 	// step function for v<0 
real stpWFpos; 	// step function for R<Rmax - the discontinuous step function is used
real stpWFneg; 	// step function for R>Rmin - the discontinuous step function is used
real swSENpos; 	// 'Switching sensitivity' function for v>0
real swSENneg; 	// 'Switching sensitivity' function for v>0
real vin; 		// Input voltage parameter 
real Rpos;		// Positive analytical time response expression
real Rneg;		// Negative analytical time response expression
real Rresp;		// Full time response expression 
real IVpos;		// Positive branch of the IV relationship 
real IVneg;		// Negative branch of the IV relationship  
real IV;		// Full IV relationship


//Local variables
real first_iteration; 
real R_last;		 
real dt;
real it;

//Implementation of the 'switching sensitivity' function
analog function real sense_fun;
input A,in;
real A,in;
begin
sense_fun=A*(-1+exp(abs(in)));
end
endfunction

//Implementation of the discontinuous step function
analog function integer stp; //Stp function
	real arg; input arg;
	stp = (arg >= 0 ? 1 : 0 );
endfunction


analog begin

//For the first iteration use 'Rinit' for the initial device resistance.
//'it' is a parameter that helps track each time-step duration
if (first_iteration==0) begin
it=0;
R_last=Rinit;
end

////Tracks the timestep duratio
dt=$realtime-it;

//Assigns the voltage applied at the terminals of the device on 'vin'
vin=V(p,n);

//Implementations of function 'r' for v>0 and v<0
Rmp=a0+a1*vin;
Rmn=b0+b1*vin;

//Continuous approximations of the step function which apply
//on the input voltage 'vin' 
//'limexp' is used instead of 'exp' to prevent numerical overflows 
//imposed by the stp functions
stpPOSv=1/(1+limexp(-vin/c1)); 
stpNEGv=1/(1+limexp(vin/c1)); 

//Step functions for the window function expression 'f'
stpWFpos=stp(s*(Rmp-R_last));
stpWFneg=stp((-s)*(Rmn-R_last));

//The 'switching sensitivity' functions
swSENpos=sense_fun(Ap,vin);
swSENneg=sense_fun(An,vin);

//Time-evolution expressions of device resistance for v>0 and v<0
Rpos=(R_last+swSENpos*Rmp*(Rmp-R_last)*stpWFpos*dt)/(1+swSENpos*(Rmp-R_last)*stpWFpos*dt);
Rneg=(R_last+swSENneg*Rmn*(Rmn-R_last)*stpWFneg*dt)/(1+swSENneg*(Rmn-R_last)*stpWFneg*dt);

//Time-evolution expression of device resistance
Rresp=Rpos*stpPOSv+Rneg*stpNEGv;     

//The two branches of the I-V expression
IVpos=ap*(1/Rresp)*sinh(bp*vin);
IVneg=an*(1/Rresp)*sinh(bn*vin);

//Device I-V expression
IV=IVpos*stpPOSv+IVneg*stpNEGv;

//Current flowing through the modules port
I(p, n)<+ IV; // Ohms law

//Updates the initial resistance for the next iteration
R_last=Rresp;

//Local parameters update
first_iteration=1;
it=$realtime;

end

endmodule