%***************Runge Kutta 4-step algorithm.*****************************
function [ret_s,ret_n]=runge4(current,sdensity,ndensity,tincrement,data)
% runge kutta 4 step algorithm

echargec	= 1.6021e-19;	  		%electron charge [Coulomb]
echarge  = echargec*1.0e9*1.0e3  ; 	%electron charge		
vlightcm	= 2.997925e10;  		%velocity of light [cm s-1]
vlight   = vlightcm/1.0e9/1.0e-7 ;	%velocity of light [nm s-1]	

%start reading parameters from the main program
gamma_cons=data(1);
gslope=data(2);
n0density=data(3);
epsi=data(4);
kappa=data(5);
Bcons=data(6);
Anr=data(7);
Ccons=data(8);
volume=data(9);
beta=data(10);
%finish reading parameters from the main program.

%Runge-Kutta 4 step integrator.
	temps	= sdensity;
	tempn	= ndensity;

%photon and carrier number in the 1st step of Runge-Kutta method.
	gn		= gamma_cons*gslope*(ndensity-n0density)*(1.0-epsi*sdensity);
	netgain	= gn-kappa;
   
  	Rsp  = Bcons * ndensity * ndensity;
	temp = Anr * ndensity + Rsp + Ccons * ndensity * ndensity * ndensity ;
   ret_val = current / echarge / volume - temp - gn * sdensity ;
   nk1	= tincrement*ret_val;
   
   ret_val = netgain * sdensity + beta * Rsp ;
	sk1	= tincrement*ret_val;

	ndensity = nk1*0.5+tempn;
	sdensity = sk1*0.5+temps;

%photon and carrier number in the 2nd step of Runge-Kutta method.
	gn		= gamma_cons*gslope*(ndensity-n0density)*(1.0-epsi*sdensity);
	netgain	= gn-kappa;
   
  	Rsp  = Bcons * ndensity * ndensity;
	temp = Anr * ndensity + Rsp + Ccons * ndensity * ndensity * ndensity ;
   ret_val = current / echarge / volume - temp - gn * sdensity ;
   nk2	= tincrement*ret_val;
   
   ret_val = netgain * sdensity + beta * Rsp ;
	sk2	= tincrement*ret_val;

	ndensity = nk2*0.5+tempn;
	sdensity = sk2*0.5+temps;

%photon and carrier number in the 3rd step of Runge-Kutta method.
	gn		= gamma_cons*gslope*(ndensity-n0density)*(1.0-epsi*sdensity);
	netgain	= gn-kappa;
   
  	Rsp  = Bcons * ndensity * ndensity;
	temp = Anr * ndensity + Rsp + Ccons * ndensity * ndensity * ndensity ;
   ret_val = current / echarge / volume - temp - gn * sdensity ;
   nk3	= tincrement*ret_val;
   
   ret_val = netgain * sdensity + beta * Rsp ;
	sk3	= tincrement*ret_val;

	ndensity = nk3+tempn;
	sdensity = sk3+temps;

%photon and carrier number in the 4th step of Runge-Kutta method.
%double gain(double freq, double sdensity, double ndensity, double Ef, 
%double En, double beta)
	gn		= gamma_cons*gslope*(ndensity-n0density)*(1.0-epsi*sdensity);
	netgain	= gn-kappa;
   
  	Rsp  = Bcons * ndensity * ndensity;
	temp = Anr * ndensity + Rsp + Ccons * ndensity * ndensity * ndensity ;
   ret_val = current / echarge / volume - temp - gn * sdensity ;
   nk4	= tincrement*ret_val;
   
   ret_val = netgain * sdensity + beta * Rsp ;
	sk4	= tincrement*ret_val;

	ret_s = temps+(sk1+2.0*sk2+2.0*sk3+sk4)/6.0;
	ret_n = tempn+(nk1+2.0*nk2+2.0*nk3+nk4)/6.0;

	sdensity=ret_s;
	ndensity=ret_n;

%Runge-Kutta function ends here.