Laser Rate Equation Modelling

67 views (last 30 days)
Donal Mckay on 23 Jan 2024
Answered: Alan Stevens on 24 Jan 2024
Hi,
I am modelling a laser diode using the laser rate equations. The coupled rate equations are solved in the laserRateEquations function. These equations describe the carrier density (N) and photon density (S). In the for loop in the test script I am also solving the Phase shift and the Output Power.
I keep getting the following error when I try to solve the Phase equation:
Array indices must be positive integers or logical values.
Error in test (line 62)
Phase(i) = (2.981/2)*(parameters.optical_conFactor*parameters.g_po(y(i,1)- ...
I have no issue when calculating the Power so I am confused as I am using the same logic.
I am unsure why this is happening and would appreciate any help on this. Thanks
Here is my code:
%Script
% Define model parameters
temperature = 25;
wavelength = 1.312e-6;
parameters.frequency = 3e8/wavelength;
parameters.I = 18e-3;
parameters.e = 1.602e-19;
parameters.V = 9e-11;
parameters.tau_n = 3e-9;
parameters.tau_p = 1e-12;
parameters.beta_sp = 4e-4;
parameters.g_po = 3e-6;
parameters.g_comp = 3.4e-17;
parameters.N_t = 1.2e18;
parameters.total_quantum_efficiency = 0.2;
parameters.plancksconstant = 6.624e-34;
parameters.optical_conFactor = 0.2;
% Set initial conditions and time span
y0 = [1e15; 0;];
tspan = [0 1e-9];
% Vary injection current values
injection_current_values = linspace(0, 20e-3, 3585); % Example: Varying from 0 to 1 mA
Power = zeros(size(injection_current_values));
Phase = zeros(size(injection_current_values));
% Solve the rate equations for each injection current value
for i = 1:length(injection_current_values)
% Set the current injection
parameters.input_power = injection_current_values(i);
% Solve the differential equations
[t, y] = ode45(@(t, y) laserRateEquations(t, y, parameters), tspan, y0);
Phase(i) = (2.981/2)*(parameters.optical_conFactor*parameters.g_po(y(i,1)- ...
parameters.N_t)-(1/parameters.tau_p));
Power(i) = (y(i,2)*parameters.V*parameters.total_quantum_efficiency* ...
parameters.plancksconstant*parameters.frequency)/ ...
(2*parameters.optical_conFactor*parameters.tau_p);%Output Power
end
%Laser Rate Equation Function
function dydt = laserRateEquations(t, y, parameters)
% Lang-Kobayashi model for carrier density (N) and photon density (S)
I = parameters.I;
e = parameters.e;
V = parameters.V;
tau_n = parameters.tau_n;
tau_p = parameters.tau_p;
beta_sp = parameters.beta_sp;
g_po = parameters.g_po;
g_comp = parameters.g_comp;
N_t = parameters.N_t;
% Rate equations
dydt = zeros(2, 1);
dydt(1) = I/(e*V) - (y(1)/tau_n) - g_po*((y(1) - N_t)/(1+y(2)))*y(2); %Carrier Density (N)
dydt(2) = g_po*((y(1) - N_t)/(1+g_comp*y(2)))*y(2) ...
- (y(2)/tau_p) + (beta_sp*(y(1)/tau_n)); %Photon Density (S)
end

Alan Stevens on 24 Jan 2024
Like this? (Look at lines 27, 28 and 38)
%Script
% Define model parameters
temperature = 25;
wavelength = 1.312e-6;
parameters.frequency = 3e8/wavelength;
parameters.I = 18e-3;
parameters.e = 1.602e-19;
parameters.V = 9e-11;
parameters.tau_n = 3e-9;
parameters.tau_p = 1e-12;
parameters.beta_sp = 4e-4;
parameters.g_po = 3e-6;
parameters.g_comp = 3.4e-17;
parameters.N_t = 1.2e18;
parameters.total_quantum_efficiency = 0.2;
parameters.plancksconstant = 6.624e-34;
parameters.optical_conFactor = 0.2;
% Set initial conditions and time span
y0 = [1e15; 0;];
tspan = [0 1e-9];
% Vary injection current values
injection_current_values = linspace(0, 20e-3, 3585);% Example: Varying from 0 to 1 mA
Power = zeros(size(injection_current_values))'; %%%%%%%%%%%%%
Phase = zeros(size(injection_current_values))'; %%%%%%%%%%%%%
% Solve the rate equations for each injection current value
for i = 1:length(injection_current_values)
% Set the current injection
parameters.input_power = injection_current_values(i);
% Solve the differential equations
[t, y] = ode45(@(t, y) laserRateEquations(t, y, parameters), tspan, y0);
Phase(i) = (2.981/2)*(parameters.optical_conFactor*parameters.g_po*(y(i,1)- ... %%%%%%%%%%%%%
parameters.N_t)-(1/parameters.tau_p));
Power(i) = (y(i,2)*parameters.V*parameters.total_quantum_efficiency* ...
parameters.plancksconstant*parameters.frequency)/ ...
(2*parameters.optical_conFactor*parameters.tau_p);%Output Power
end
figure
subplot(2,1,1)
plot(t,Phase),grid
xlabel('t'),ylabel('Phase')
subplot(2,1,2)
plot(t,Power),grid
xlabel('t'),ylabel('Power')
%Laser Rate Equation Function
function dydt = laserRateEquations(~, y, parameters)
% Lang-Kobayashi model for carrier density (N) and photon density (S)
I = parameters.I;
e = parameters.e;
V = parameters.V;
tau_n = parameters.tau_n;
tau_p = parameters.tau_p;
beta_sp = parameters.beta_sp;
g_po = parameters.g_po;
g_comp = parameters.g_comp;
N_t = parameters.N_t;
% Rate equations
dydt = zeros(2, 1);
dydt(1) = I/(e*V) - (y(1)/tau_n) - g_po*((y(1) - N_t)/(1+y(2)))*y(2); %Carrier Density (N)
dydt(2) = g_po*((y(1) - N_t)/(1+g_comp*y(2)))*y(2) ...
- (y(2)/tau_p) + (beta_sp*(y(1)/tau_n)); %Photon Density (S)
end