%Chapt9Exercise5.m
% solves single mode laser rate equations
% using 4th order Runge-Kutta method for fixed time increment and plots
% light output, L, and carrier density, n, 
% as a function of injection current, I.

%	gain=g=gamma*gslope*(n-n0)*(1-epsi*s)
%	carrier density, n, current, I, and photon density, s, are related via 
%	fcn=n'=(I/ev)-(n/tau_n)-(g*s)
%	fcs=s'=(g-K)*s+beta*Rsp

%  start the program and clear all the previous definitions.
clear;
clf;
FSplot=14;      %font size for plots 
LWplot=2;       %line width for plots
TLplot=0.03;    %tick length for plots
% declare constants
hconstjs=6.6262e-34;   		%Plancks constant [J s]
echargec=1.6021e-19;	  	%electron charge [Coulomb]
vlightcm=2.997925e10;  		%velocity of light [cm s-1]
pi2=2*pi;			  	  	%value of 2pi

ngplot=100;					%number of points plotted
nsteps=2000;				%number of steps in time integration
curentI =zeros(1,ngplot);
carrierN=zeros(1,ngplot);
photonL =zeros(1,ngplot);
%********************  start the main program ********************************/
% start reading parameters from the data file datainLI.txt
fp = fopen('datainLI.txt','r');
for i=1:1:20
	x(i) = fscanf(fp, '%e', 1);
end   
fclose(fp);

ngroup     = x(1);			%refractive index

clength    = x(2);			%cavity length (cm)
thick      = x(3);			%thickness of active region (cm)
width      = x(4);			%width of active region (cm)

tincrement = x(5);			%time increment (s)
initialn   = x(6);			%initial value of carrier density (cm-3)

Anr        = x(7);			%non-radiative recombination rate (s-1)
Bcons      = x(8);			%radiative recombination coefficent (cm3 s-1) 
Ccons      = x(9);			%non-linear recombination coefficient (cm6 s-1)

n0density  = x(10);			%transparency carrier density (cm-3)
gslope     = x(11);			%optical gain-slope coefficient (cm2 s-1)
epsi       = x(12);			%gain compression (cm3)
beta       = x(13);			%spontaneous emission coefficient
gamma_cons = x(14);			%optical confinement factor
wavelength = x(15);			%optical emission wavelength

mirrone    = x(16);			%optical reflectivity from first mirror
mirrtwo    = x(17);			%optical reflectivity from second mirror
alfa_i     = x(18);			%internal optical loss (cm-1)

Imin	     = x(19);		%start value of current
Imax   	  = x(20);			%end value of current

% finish reading parameters from the data file
% initialize sdensity, ndensity and calculate total optical loss [cm-1]
	sdensity = 0.0;
	ndensity = initialn * 1.0e-21;

	alfa_m   = log ( 1.0/(mirrone * mirrtwo) ) / ( 2.0 * clength ) ;
	alfasum  = alfa_m + alfa_i ;

% rescale integrator: time[ns], current[mA], length[nm] and 1nm3=1e-21cm3
	echarge   	= echargec*1.0e9*1.0e3  ; 	% electron charge		
   vlight    	= vlightcm/1.0e9/1.0e-7 ;	% velocity of light [m/s] to [nm s-1]	

	tincrement 	= tincrement * 1.0e9 ;
	gslope    	= gslope * 1.0e14 ;
   n0density 	= n0density * 1.0e-21 ;
   
% length scale in nm.
   clength   	= clength / 1.0e-7 ;
	width     	= width / 1.0e-7 ;
   thick			= thick / 1.0e-7 ;
   
%  recombination constants   
	Anr       	= Anr / 1.0e9 ;
	Bcons   		= Bcons / (1.0e9*1.0e-21) ;
   Ccons   		= Ccons / (1.0e9*1.0e-42) ;
   
	epsi    		= epsi / 1.0e-21 ;
	gslope  		= gslope * vlight / ngroup ;
	alfasum 		= alfasum * 1.0e-7 ;
	kappa   	 	= alfasum * vlight / ngroup ;
	volume  	 	= clength * width * thick ;

% scale light output to mW from mirror one
	tmp  	 = hconstjs * (vlightcm * 1.0e-2)/(wavelength * 1.0e-6);	% use m [J]	
	tmp	 = tmp * alfa_m * vlightcm;						% use cm [s-1]	*/
	tmp	 = tmp * 1000.0*volume/ngroup;					% use mW and nm [nm3]*/
	lightout=tmp*(1.0-mirrone)/(2.0-mirrone-mirrtwo);% output mirror one*/

%********************  the main integration loop  ****************************

data(1)=gamma_cons;
data(2)=gslope;
data(3)=n0density;
data(4)=epsi;
data(5)=kappa;
data(6)=Bcons;
data(7)=Anr;
data(8)=Ccons;
data(9)=volume;
data(10)=beta;

photonold=0;
current=0;
dcurrent=(Imax-Imin)/ngplot;

for j=1:ngplot
   current=current+dcurrent;
   
   for i = 1:nsteps;   
   	[sdensity, ndensity]=runge4( current, sdensity, ndensity, tincrement, data);
   	carrier = ndensity * 1.0e21/1.0e18	;
      photon = (lightout*sdensity);
      photonold=photon;
   end;
currentI(j)=current;
photonL(j)=photon;
carrierN(j)=carrier;


end

close all;

figure(1);
hold on;
%grid on;
plot(currentI,photonL,'r-');
temp=['Light output as function of injection current'];
title(temp, 'fontsize', 14)
xlabel(['Current, I (mA)'],'fontsize', 14);
ylabel(['Light, L (mW/facet)'],'fontsize', 14);
hold off;

figure(2);
%grid on;
plot(currentI,log(photonL),'r-');
title(temp, 'fontsize', 14)
xlabel(['Current, I (mA)'],'fontsize', 14);
ylabel(['log(light), log(L) (mW/facet)'],'fontsize', 14);

figure (3);
hold on;
temp2=['Carrier density as function of injection current'];
title(temp2, 'fontsize', 14)
xlabel(['Current, I (mA)'],'fontsize', 14);
ylabel(['Carrier density, n (\times 10^{18} cm^{-3})'],'fontsize', 14);
plot(currentI,carrierN,'b');
hold off;


figure (4);
set(gca,'Fontsize',FSplot);                     %set current gca font size
hl1 = line(currentI,carrierN,'Linewidth',LWplot,'Color','k');
ax1 = gca;                                      %ax1 is first axis
set(ax1,'XColor','k','YColor','k','Linewidth',LWplot,...
    'TickLength',[TLplot TLplot]);
ylabel(['Carrier density, \it n \rm (\times 10^{18} cm^{-3})']);

ax2 = axes('Position',get(ax1,'Position'),...   %ax2 is second axis
'YAxisLocation','right','XAxisLocation','bottom',...
'Color','none',...
'XColor','k','YColor','r','Linewidth',LWplot,'TickLength',[TLplot TLplot]);
set(gca,'Fontsize',FSplot);                     %set current gca font size
hl2 = line(currentI,photonL,'Linewidth',LWplot,'Color','r','Parent',ax2);

temp2=['Carrier density and light output as function of injection current'];
title(temp2);
xlabel(['Current, \it I \rm (mA)']);
ylabel(['Light, \it L \rm (mW/facet)']);


% [AX,H1,H2]=plotyy(log10(currentI),carrierN,log10(currentI),photonL,'plot');
% set(get(AX(1),'Ylabel'),'String','Carrier density, n (\times 10^{18} cm^{-3})') ;
% set(get(AX(2),'Ylabel'),'String','Light, L (mW/facet)') ;
% temp2=['Carrier density and light output as function of injection current'];
% title(temp2, 'fontsize', 12);
% xlabel(['Log_{10} (current), log_{10} I (mA)'],'fontsize', 12);
% set(H1,'Color','k','LineStyle','-');
% set(H2,'Color','g','LineStyle','-');
% ax1 = gca;
% set(ax1,'XColor','k','YColor','k');



