////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Copyright @ 2015 The Regents of the University of California
 
//The terms under which the software and associated documentation (the Software) is provided are as the following:
 
//The Software is provided "as is", without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose and noninfringement. In no event shall the authors or copyright holders be liable for any claim, damages or other liability, whether in an action of contract, tort or otherwise, arising from, out of or in connection with the Software or the use or other dealings in the Software.
 
//MIT grants, free of charge, to any users the right to modify, copy, and redistribute the Software, both within the user's organization and externally, subject to the following restrictions:
 
//1. The users agree not to charge for the MIT code itself but may charge for additions, extensions, or support.
 
//2. In any product based on the Software, the users agree to acknowledge the MIT VS Model Research Group that developed the software. This acknowledgment shall appear in the product documentation.
 
//3. The users agree to obey all U.S. Government restrictions governing redistribution or export of the software.
 
//4. The users agree to reproduce any copyright notice which appears on the software on any copy or modification of such made available to others.
 
// Agreed to by 
// S. Adair Gerke
// May 8, 2015
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



// UC Berkeley VCSEL (Vertical Cavity Surface Emitting Laser) Compact Model Version 1.00
// A compact model for VCSEL laser direct modulation dynamics, including self-heating, for use in optical interconnects
// Implemented by S. Adair Gerke and C.J. Chang-Hasnain, UC Berkeley, May 8 2015


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

module vcsel(cathode,anode,refnode,P,L,tsink,tqw);
	// Nodes:
	// Real electrical nodes (cathode, anode), real internal electical nodes (anodeint, cathodeint, midp) 
	// Reference node (refnode) for all non-electrical quantities below, represented as voltages vs. this node
	// Rate equation unknowns: (M,N,Y) as voltages, scaled by Mo, No
	// Light output: (L) voltage represents wavelength (um), (P) voltage represents power (mW)
	// Thermal nodes: (tsink, tqww) temperatures at heat sink and quantum well (K)
	electrical cathode,anode,anodeint,cathodeint,midp;
	electrical refnode,P,L,tsink,tqw,M,N,Y;

	// Branches:
	branch (anodeint,cathodeint) active;
	branch (anodeint,cathodeint) jncap; 
	branch (anode,anodeint) dbr1;
	branch (cathodeint,cathode) dbr2;
	
	// Parameters:
	parameter real version 	= 1.00 from [0:inf);

	// Electrical  parameters
	parameter real is		= 1.0e-15 from (0:1.0);  	// [A] Saturation current 
	parameter real rp 		= 90 from [0:inf); 		// [Ohm] pad resistance
	parameter real cp 		= 50e-15 from [0:inf); 	// [F] pad capacitance
	parameter real cjt 		= 500e-15 from [0:inf); 	// [F] junction capacitance
	parameter real rdbr1		= 30 from [0:inf); 		// [Ohm] DBR resistance
	parameter real rdbr2 		= 30 from [0:inf); 		// [Ohm] DBR resistance
	parameter real ideality 	= 2 from (0:inf); 			// [#] ideality factor of diode

	// Geometrical parameters
	parameter real rapt		= 4.5e-4 from (0:1.0); 	// [cm] radius of aperture
	parameter real dqw		= 8e-7 from (0:inf); 		// [cm] thickness of quantum wells
	parameter real nqw 		= 3 from (0:inf); 			// [#] number of quantum wells
	parameter real Gammaq		= 0.1 from (0:inf); 		// [#] qw/sch confinement factor (volume ratio)

	// Optical parameters
	parameter real Gamma		= 0.025 from (0:1.0); 		// [#] Mode overlap factor from quantum well, fit or use optical simulations to calculate
	parameter real Beta		= 1.69e-4 from (0:inf); 	// [s] spontaneous emission coefficient 
	parameter real vg		= 9.0E9 from (0:3.0E10); 		// [cm/s] Mode group velocity, must be < speed of light
	parameter real alphai0		= 12 from [0:inf); 		// [cm^-1] Internal modal loss
	parameter real alphatm		= 18 from (0:inf); 		// [cm^-1] Top mirror modal loss
	parameter real alphabm		= 1.6 from [0:inf); 		// [cm^-1] Bottom mirror modal loss
	parameter real etai0		= 0.849 from [0:1.0); 		// [#] injection efficiency at 0C (273.15 K)
	parameter real Epsicm		= 3.8E-17 from [0:inf); 	// [#] Photon dependence of gain compression (see Kern)
	parameter real amfracK		= 1e-3 from (-inf:inf); 		// [K^-1] Mirror reflectivity changes with temp (see Westbergh)
	parameter real EpsicmK		= 1e-19 from (-inf:inf); 	// [cm^3/K] Gain compression change with temp
	
	// Carrier dynamics parameters
	parameter real tcap		= 1e-12 from (0:inf);  	// [s] capture time to quantum well from clad. 
	parameter real tesc		= 60e-12 from (0:inf); 	// [s] escape time from quantum well to clad. 

	// Logarighmic Gain vs N parameters
	parameter real Ntrlog		= 1.5e18 from (0:inf); 	// [cm^-3] transparency density 
	parameter real Nslog		= 2.0e17 from (0:inf); 		// [cm^-3] quantum well gain parameter: coldren and corzine
	parameter real gain0 		= 4000 from (0:inf); 		// [cm^-1] gain at TrefK
	parameter real gnperK		= -3.4e-3 from (-inf:inf); 	// [1/K] Gain rolloff vs temperature
		
	// Thermo-optical parameters: Chuang, Michalzik (fig 2.17), Lascola. 
	parameter real TrefK		= 300 from (0:500);		 // [K] Temp at which thermally-dependent parameters are measured
	parameter real lambdaK		= 0.062 from (-inf:inf); 		// [nm/K] cavity wavelength change vs tempreature
	parameter real lambda0 		= 0.850 from (0:inf); 		// [um] wavelength at Tfit
	parameter real Arec0		= 1.0e8 from [0:inf); 		// [1/s] NR recombination coefficient
	parameter real Brec0		= 1.5e-10 from (0:inf); 	// [1/s*cm^3] Radiative recombination coefficient
	parameter real Crec0		= 4e-30 from [0:inf);		// [1/s*cm^6] SRH recombination coefficient
	parameter real Rtherm0		= 1.9e3 from [0:inf); 		// K/W, thermal resistance to heat sink
	parameter real Ttherm 		= 1e-6 from [0:inf); 		//sec, thermal time constant
	
	// Thermal dependence of VCSEL parameters (Baveja)
	parameter real alphaiK		= 0.035 from (-inf:inf); 		// [1/cm/K]
	parameter real etaiK1		= 0.0022 from (-inf:inf);	// [K^-1] Quantum well capture vs. temperature, 1st deriv
	parameter real etaiK2		= -3.3e-5 from (-inf:inf); 	// [K^-2] Quantum well capture vs. temperature, 1st deriv: Sign very important
	parameter real RthermK		= 0 from (-inf:inf); 			// [(K/W)/K] Change in Rtherm vs temp (= 5 in Baveja)
	parameter real BrecK1 		= -7.5e-13 from (-inf:inf); 	// [1/s*cm^3/K] thermal fit of Brec
	parameter real BrecK2 		= 3.125e-15 from (-inf:inf); // [1/s*cm^3/K^2] thermal fit of Brec
	
	// Scaling to make internal node voltages, currents reasonable magnitude. Improves convergence, but does not change model equations
	parameter real No 		= 1e18 from (0:inf); 
	parameter real Mo 		= 1e16 from (0:inf); 
	parameter real Po 		= 1e15 from (0:inf); 
 	parameter real delta		= 5e-8 from (0:inf); 	

	// Variable declarations
	real Volcm,Brec,TempQW, EpsicmF,lambdaum,gcf,gain,theta, alphai,alpham,Tdc,Tc,S,tn,frec,heatgen,etai,Tfit,Rtherm,tp; 
	real IY1,IY2,IY3,IN1,IN2,IN3,IN4,IM1,IM2,IM3; //Rate terms

	analog begin
		//*********VARIABLE CALCULATIONS
		Volcm		= `M_PI *rapt*rapt*dqw*nqw; 		//[cm^3] volume of ACTIVE 

		// Thermally Parameter Adjustments
		TempQW	= V(tqw,refnode);
		Tfit 			= TempQW - TrefK;
		Tc 			= TempQW - `P_CELSIUS0;		 //Celsius quantum well temp	
		alphai 		= alphai0 + alphaiK*Tfit;
		etai			= etai0 + Tc*etaiK1 + Tc*Tc*etaiK2; 
		Rtherm		= Rtherm0 + RthermK*Tfit;
		Brec 		= Brec0  + BrecK1*Tfit + BrecK2*Tfit*Tfit;		
		alpham		= (alphatm + alphabm)*(1.0 + amfracK*Tfit); 		//[cm^-1] total modal mirror loss	
		EpsicmF 	= Epsicm + Tfit*EpsicmK;

		tp	= 1.0/(alpham + alphai)/vg; 				// [sec] photon lifetime
		tn 	= 1.0/(Arec0 + Brec*V(N,refnode)*No + Crec0*V(N,refnode)*V(N,refnode)*No*No) ; // [sec] carrier lifetime
		frec	= Brec*V(N,refnode)*No*tn + 0.0001; 		//Fraction of recombination raditive; small number to avoid going negative
		lambdaum	= lambda0 + lambdaK*Tfit;		// [um] wavelength due to thermal tuning
		theta	= (lambdaum*1e-6*Gamma)/(`P_H*`P_C*Volcm*vg*alphatm); 	// [cm^-3/W]  normalizing factor, see Javro
			
		// Gain calculations
		S		= theta*pow(V(Y,refnode) + delta,2); 	// [cm^-3] photon concentration
		gcf 		= (1.0/(1.0+EpsicmF*S));				// gain compression factor term
		gain 	= gain0*gcf*(ln((V(N,refnode)*No + Nslog)/(Ntrlog + Nslog)) + gnperK*Tfit); //[cm^-1] optical gain in quantum well

		//********* RATE EQUATIONS:  term calculations. 
		// M equation: Electrons in SCH
		IM1	= I(active)*(Gammaq*etai/`P_Q/Volcm)/Mo;	// [cm^-3/s] injection into sch. 
		IM2	= -1.0*V(M,refnode)/tcap;					// [cm^-3/s]capture into qw
		IM3	= Gammaq*V(N,refnode)*No/tesc/Mo;		// [cm^-3/s]re-escape from qw
		
		// N equation: Electrons in quantum well
		IN1 	= V(M,refnode)*Mo/(Gammaq*tcap)/No; 	// [cm^-3/s] Capture into qw from sch
		IN2	= -V(N,refnode)*(1.0/tn + 1.0/tesc); 			// [cm^-3/s] recombination (ALL KINDS) and escape
		IN3 	= -(gain*vg*S)/No; 						// [cm^-3/s] stim reombination
		IN4 = limexp(-V(N,refnode)*10000) - 1.0; 		// [cm^-3/s] Serves to keep N positive. 
		
		// Y equation: square root of photon population; Y has dimensions of electric field amplitude per volume
		IY1 	= Gamma*gain*vg*(V(Y,refnode) + delta);	//Stim emission into lasing mode
		IY2	= Gamma*Beta*V(N,refnode)*No/(theta*tn*(V(Y,refnode) + delta))*frec; //Spont emission into lasing mode
		IY3 	= -(V(Y,refnode) + delta)/tp;				//Photon escape	
		
		// THERMAL equations
		heatgen = I(active)*V(active) + I(dbr1)*V(dbr1) + I(dbr2)*V(dbr2) - V(P,refnode); //Heat generated: electrial minus optical power
		Tdc	= V(tsink) + heatgen*Rtherm*0.5*(1.0+tanh(1e4*heatgen)); 	//Avoids negative heat flow
		
		//********* CONTRIBUTIONS

		// optical contributions: core rate equations
		I(M,refnode)		<+ IM1 + IM2 + IM3;
		I(M,refnode)		<+ -ddt(V(M,refnode));
		I(N,refnode) 		<+ IN1 + IN2 + IN3 + IN4;
		I(N,refnode) 		<+ -ddt(V(N,refnode));
		I(Y,refnode) 		<+ IY1 + IY2 + IY3;
		I(Y,refnode) 		<+ -ddt(V(Y,refnode))*2.0; 			//factor of 2 from y=(s+d)^2 chain rule		
		V(P,refnode) 	<+ pow((V(Y,refnode) + delta),2); 
		V(L,refnode) 		<+ lambdaum;
		
		// electrical node contributions
		I(dbr1) 			<+ V(anode,anodeint)/rdbr1;
		I(active) 			<+ is*(limexp(V(anodeint,cathodeint)/(V(tsink)*(`P_K/`P_Q)*ideality)) - 1.0); // Diode equation
		I(jncap) 			<+ ddt(cjt*V(anodeint,cathodeint));
		I(dbr2) 			<+ V(cathodeint,cathode)/rdbr2; 
		I(anode,midp) 	<+ V(anode,midp)/rp; 				// this is in SERIES with the pad leakage path
		I(midp,cathode) 	<+ ddt(cp*V(midp,cathode)); 		// pad capacitance

		// thermal contributions
		V(tqw,refnode)	<+ Tdc - ddt(Ttherm*V(tqw,refnode)); 

	end //of analog begin
endmodule