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

`define E0 8.85e-12  // Permittivity of Vacuum [F/m]
`define MPRoo(nam,def,uni,lwr,upr,des) (*units=uni,                   desc=des*) parameter real    nam=def from(lwr:upr);
`define MPIcc(nam,def,uni,lwr,upr,des) (*units=uni,                   desc=des*) parameter integer nam=def from[lwr:upr];


module ferroelectric (a,b,c,d); // Extra terminals c and are added to obtain varible "Q" and "V(aux)" in SPICE simulator
  inout a, b,c,d;
  electrical a, b, delay, aux, c, d;
  
  `MPRoo( EFE0             ,40             ,"F/m"            ,1          ,inf         ,"Background Permittivity of the Ferroelectric layer" )
  `MPRoo( TFE              ,10e-9          ,"m"              ,1e-9       ,inf         ,"Ferroelectric Thickness" )
  `MPRoo( EC               ,1.2e8          ,"V/m"            ,1e7        ,inf         ,"Coercive Field" )
  `MPRoo( PR               ,0.185           ,"C/m^2"          ,1e-6       ,inf         ,"Remnant Polarization" )
  `MPRoo( PS               ,0.19            ,"C/m^2"          ,PR+1e-6       ,inf         ,"Saturation Polarization" )
  `MPRoo( VMAX             ,3.0            ,"V"              ,0          ,inf         ,"Maximum Applied Voltage" )
  `MPRoo( RDELAY           ,1              ,"ohm"            ,1e-3       ,inf         ,"Resistor value to produce delayed version" )
  `MPRoo( CDELAY           ,1e-15          ,"F"              ,1e-18      ,inf         ,"Capacitor value to produce delayed version" )
  `MPRoo( TAUV             ,0.1e-6         ,"s"              ,1e-9       ,inf         ,"Relaxation Time for auxiliary Voltage" )
  `MPRoo( TAUP             ,1e-10          ,"s"              ,0      ,inf         ,"Dielectric Relaxation Time" )
  `MPRoo( KN               ,1e10           ,"C.s/V"          ,1          ,inf         ,"Coupling Coefficient" )
  `MPIcc( INITYPE          ,1              ,""               ,0          ,1           ,"Initialization 1:Forward loop and 0: Reverse loop" )

  real type,m_up,P_up,P_down,m_down,slope,path,efe,alpha;
  real Vfe,Paux,V_turn_up,V_turn_down,P_turn_down,P_turn_up,Q, P;

analog begin

    Vfe=V(a,b); // Voltage across the Ferroelectric Layer
	efe=EFE0*`E0;
	slope=(1.0/(2.0*EC))*ln((PS+PR)/(PS-PR)); // Equation (4)
	
// Initialization	
   @(initial_step)
   begin
	  P_down = 0; 
	  P_up = 0;
	  alpha=1.0;
      type=INITYPE;
	  V_turn_up=VMAX;
	  V_turn_down=-VMAX;
	  P_turn_up=PS*tanh(slope*(VMAX/TFE-EC));
	  P_turn_down=PS*tanh(slope*(-VMAX/TFE+EC));
	  path=1;
	  m_down=1;
	  m_up=1;
   end
	
///Tracing of turning points in auxiliary voltage using R-C circuit (Figure 4(a))
  I(delay)<+ V(delay)/RDELAY + CDELAY*ddt(V(delay));
  I(delay)<+ -alpha*V(aux);
	
// Reverse Loop
  if (type==0)
				begin				
					Paux=m_down*PS*tanh(slope*(V(aux)/TFE+EC))+P_down ; // Equation (5)
					P=Paux+(efe*V(aux))/TFE;
					
//Condition for forward loop            				
					if (V(aux) > V(delay))
							begin
								type=1;
								V_turn_down=V(aux);
								P_turn_down=Paux;					
					end 
					else if(V(aux)<V_turn_down)
							begin
								m_down=(P_turn_down-PS*tanh(slope*(-VMAX/TFE+EC)))/(PS*(tanh(slope*(V_turn_down/TFE+EC))-tanh(slope*(-VMAX/TFE+EC)))+1e-6); //Equation (14) // Added small value (1e-6) in the denominator to avoid divide by zero error.
								P_down=P_turn_down-m_down*PS*(tanh(slope*(V_turn_down/TFE+EC))); //Equation (15)
								path=0;
					end
					else
							begin
								m_down=(P_turn_up-P_turn_down)/(PS*(tanh(slope*(V_turn_up/TFE+EC))-tanh(slope*(V_turn_down/TFE+EC)))); // Equation (12)
								P_down=P_turn_up-m_down*PS*(tanh(slope*(V_turn_up/TFE+EC))); // Equation (13)
								path=1;
					end
					
					
				end
// Forward Loop	
			else
				begin
				if(path==0)
							begin
								m_up=(P_turn_down-PS*tanh(slope*(VMAX/TFE-EC)))/(PS*(tanh(slope*(V_turn_down/TFE-EC))-tanh(slope*(VMAX/TFE-EC)))+1e-6); // Equation (16) // Added small value (1e-6) in the denominator to avoid divide by zero error.
								P_up=P_turn_down-m_up*PS*(tanh(slope*(V_turn_down/TFE-EC))); // Equation (17)
								Paux=m_up*PS*tanh(slope*((V(aux))/TFE-EC))+P_up; 
				                P=Paux+(efe*V(aux))/TFE;
					end 
					else
							begin
								m_up=(P_turn_down-P_turn_up)/(PS*(tanh(slope*(V_turn_down/TFE-EC))-tanh(slope*(V_turn_up/TFE-EC)))+1e-6); // Added small value (1e-6) in the denominator to avoid divide by zero error.
								P_up= P_turn_down-m_up*PS*(tanh(slope*(V_turn_down/TFE-EC)));
								Paux=m_up*PS*tanh(slope*((V(aux))/TFE-EC))+P_up; // Equation (3)
				                P=Paux+(efe*V(aux))/TFE;
								if(V(aux)>V_turn_up)
								begin
								path=0;
								P_turn_down=Paux;
								V_turn_down=V_turn_up;
								end
								
					end		 

				
// Condition for reverse loop               				
					if (V(aux) < V(delay))
							begin
								type=0;
								V_turn_up=V(aux);
								P_turn_up=Paux;
								m_down=(P_turn_up-P_turn_down)/(PS*(tanh(slope*(V_turn_up/TFE+EC))-tanh(slope*(V_turn_down/TFE+EC)))+1e-6); // Equation (10) // Added small value (1e-6) in the denominator to avoid divide by zero error.
								P_down=P_turn_up-m_down*PS*(tanh(slope*(V_turn_up/TFE+EC))); // // Equation (11)
					end
						
					
			end
V(aux)<+Vfe-TAUV*ddt(V(aux));	// Auxiliary voltage (Equation (1))
Q= P - TAUP*((ddt(P))/(1+(KN/TFE)*(ddt(V(aux)))));	// Ferroelectric charge density (Equation (6))			
//Q=P;	
I(a,b)<+ ddt(Q);
I(c)<+ Q;
I(d)<+ V(aux);		
end

endmodule