Laser Rate Equation Modelling

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

Answer 1

[Alan Stevens]

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