%*****************Runge Kutta 4 steps algorithm incorporating noise.*****************************
function [ret_s, ret_n] = runge4_n_sum(current, s, n, tincrement, data,time)
% runge kutta 4 step algorithm

ngroup=4;                   %refractive index
echarge	= 1.6021e-19;       %electron charge [Coulomb]
echargec=echarge;		
vlightcm	= 2.997925e10;  %speed of light [cm s-1]
kappa=data(1);              %s-1, total loss rate
volume=data(2);             %cm3

gamma_cons=0.25;            %mode confinement factor
gslope=2.5E-16 ;            %[cm2 s-1] gain slope coefficient
Anr=2e+8;                   %non-radiative recombination constant
Ccons=1e-29;                %cm6 non linear recombination constant
Bcons=1E-10;                %cm3 radiative recombination constant
beta=1E-4;
n0=(1E+018)*volume;         %carrier number at transparency
epsi=5e-18;                 %cm3
time_s=10e-9;               %time after which noise is introduced
Fne=0;
Fse=0;

% Runge-Kutta 4 step integrator.
	temps	= s;
	tempn	= n;
    
% photon and carrier # in the 1st step of Runge-Kutta method.
	gn			= gamma_cons * gslope *(( n - n0)/volume)* (1.0 - ( epsi * s/volume) );
    gain_varnce = gamma_cons * gslope * (( n - n0)/volume)* (1.0 - ( epsi * s/volume) );    
    gn=gn * vlightcm / ngroup;
    gain_varnce=gain_varnce * vlightcm / ngroup;	
    netgain	= gn-kappa;   
  	Rsp  = Bcons * n*n/volume;       
    temp = ( Anr *n ) + Rsp + ( Ccons * n^3 )/volume^2;    
    if time>time_s   
	Vni =sqrt(abs( ( current /  echargec ) + temp +( gain_varnce * s )));
    count_i=0;
    tmp1n = -1000*n;
    tmptx=tincrement;
while count_i < 999 & tmp1n+ n < 0      
    Fne = randn*Vni/sqrt(tmptx);
    tmp1n = Fne*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1n+ n < 0
    Fne=0;
end
    end
%     Fne=sqrt(abs((2*temp+2*gn*s)/tincrement))   
   ret_val = current / echarge  - temp - gn * s + Fne;
   
   nk1	= tincrement*ret_val;
 if time>time_s  
   Vsi = sqrt(abs(gain_varnce*s+kappa*s + beta*Rsp));
 count_i=0;
    tmp1s = -1000*s;
    tmptx=tincrement;
while count_i < 999 & tmp1s+ s < 0      	 
    Fse = randn*Vsi/sqrt(tmptx);
    tmp1s = Fse*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1s+ s < 0
    Fse=0;
end
end 
%    Fse=sqrt(abs((2*gn*s+2*beta*Rsp)/tincrement))  
    ret_val = netgain * s + beta * Rsp + Fse;
    sk1	= tincrement*ret_val;
	n = nk1*0.5+tempn;
	s = sk1*0.5+temps;
   
% photon and carrier # in the 2nd step of Runge-Kutta method.
	gn			= gamma_cons * gslope *(( n - n0)/volume)* (1.0 - ( epsi * s/volume) ); 
    gain_varnce = gamma_cons * gslope * (( n - n0)/volume)* (1.0 - ( epsi * s/volume) );	
    gn=gn * vlightcm / ngroup;
    gain_varnce=gain_varnce * vlightcm / ngroup;	
    netgain	= gn-kappa;  
  	Rsp  = Bcons * n*n/volume;   
    temp = ( Anr *n ) + Rsp + ( Ccons * n^3 )/volume^2;   
    if time>time_s    
	Vni =sqrt( abs(( current /  echargec ) + temp +( gain_varnce * s )));
    count_i=0;
    tmp1n = -1000*n;
    tmptx=tincrement;
while count_i < 999 & tmp1n+ n < 0       	 
    Fne = randn*Vni/sqrt(tmptx);
    tmp1n = Fne*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1n+ n < 0
    Fne=0;
end
    end
%     Fne=sqrt(abs((2*temp+2*gn*s)/tincrement))   
   ret_val = current / echarge  - temp - gn * s + Fne;   
   nk2	= tincrement*ret_val;
 if time>time_s  
   Vsi = sqrt(abs(gain_varnce*s+kappa*s + beta*Rsp));
 count_i=0;
    tmp1s = -1000*s;
    tmptx=tincrement;
while count_i < 999 & tmp1s+ s < 0    	 
    Fse = randn*Vsi/sqrt(tmptx);
    tmp1s = Fse*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1s+ s < 0
    Fse=0;
end
end 
%    Fse=sqrt(abs((2*gn*s+2*beta*Rsp)/tincrement))  
    ret_val = netgain * s + beta * Rsp + Fse;
    sk2	= tincrement*ret_val;
	n = nk2*0.5+tempn;
	s = sk2*0.5+temps;

% photon and carrier # in the 3rd step of Runge-Kutta method.
	gn			= gamma_cons * gslope * (( n - n0)/volume)* (1.0 - ( epsi * s/volume) ); 
    gain_varnce = gamma_cons * gslope * (( n - n0)/volume)* (1.0 - ( epsi * s/volume) );	
    gn=gn * vlightcm / ngroup;
    gain_varnce=gain_varnce * vlightcm / ngroup;	
    netgain	= gn-kappa;   
  	Rsp  = Bcons * n*n/volume;    	    
    temp = ( Anr *n ) + Rsp + ( Ccons * n^3 )/volume^2;
    
    if time>time_s    
	Vni =sqrt( abs(( current /  echargec ) + temp +( gain_varnce * s )));
    count_i=0;
    tmp1n = -1000*n;
    tmptx=tincrement;
while count_i < 999 & tmp1n+ n < 0       	 
    Fne = randn*Vni/sqrt(tmptx);
    tmp1n = Fne*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1n+ n < 0
    Fne=0;
end
    end
%     Fne=sqrt(abs((2*temp+2*gn*s)/tincrement))   
   ret_val = current / echarge  - temp - gn * s + Fne;   
   nk3	= tincrement*ret_val;
 if time>time_s  
   Vsi = sqrt(abs(gain_varnce*s+kappa*s + beta*Rsp));
 count_i=0;
    tmp1s = -1000*s;
    tmptx=tincrement;
while count_i < 999 & tmp1s+ s < 0      	 
    Fse = randn*Vsi/sqrt(tmptx);
    tmp1s = Fse*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1s+ s < 0
    Fse=0;
end
end
%    Fse=sqrt(abs((2*gn*s+2*beta*Rsp)/tincrement)) 
   ret_val = netgain * s + beta * Rsp + Fse;
	sk3	= tincrement*ret_val;
	n = nk3+tempn;
	s = sk3+temps;
  
% photon and carrier # 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 *(( n - n0)/volume)* (1.0 - ( epsi * s/volume) ); 
    gain_varnce = gamma_cons * gslope *(( n - n0)/volume)* (1.0 - ( epsi * s/volume) );    
gn=gn * vlightcm / ngroup;
    gain_varnce=gain_varnce * vlightcm / ngroup;		
    netgain	= gn-kappa;   
  	Rsp  = Bcons * n*n/volume;   	    
    temp = ( Anr *n ) + Rsp + ( Ccons * n^3 )/volume^2;    
    if time>time_s   
	Vni =sqrt(abs( ( current /  echargec ) + temp +( gain_varnce * s )));
    count_i=0;
    tmp1n = -1000*n;
    tmptx=tincrement;
while count_i < 999 & tmp1n+ n < 0       	 
    Fne = randn*Vni/sqrt(tmptx);
    tmp1n = Fne*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1n+ n < 0
    Fne=0;
end
    end
%     Fne=sqrt(abs((2*temp+2*gn*s)/tincrement))   
   ret_val = current / echarge  - temp - gn * s + Fne;  
   nk4	= tincrement*ret_val;
 if time>time_s  
   Vsi = sqrt(abs(gain_varnce*s+kappa*s + beta*Rsp));
 count_i=0;
    tmp1s = -1000*s;
    tmptx=tincrement;
while count_i < 999 & tmp1s+ s < 0       	 
    Fse = randn*Vsi/sqrt(tmptx);
    tmp1s = Fse*tincrement;
    count_i=count_i+1;
end
if count_i>=999 & tmp1s+ s < 0
    Fse=0;
end
end 
%    Fse=sqrt(abs((2*gn*s+2*beta*Rsp)/tincrement))
  
ret_val = netgain * s + beta * Rsp + Fse;
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;

% Runge-Kutta function ends here.

