#### Sur76

##### Newbie level 1

Hello,

I am designing a optical comb to be used in DWDM systems, can anyone please help me in giving some conceptual inputs regarding how the laser can be gain switched to achieve the comb. I have a code that solves the rate equations to give photon density, carrier density and phase of the light emmitted. So can you please throw some light on this topic as I'm a bit new to this domain.

Thanking you in advance.

Sur76

P.S: The following is the code I have.

clear all

ind2 = 1;

steps = 0:0.1:20;

% Bind = 1;

% Bias_step = 16:4:24;

%

% for Biasing = Bias_step

for ind = steps

%% %%%%%%%%%%%%%%%%%%%%% Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nbits = 1024*10;

freq = ind*1e9;

%freq = 0.1*1e9;

w = 2*pi*freq;

% bitrate = 9e10;

bitrate = 300e9;

dt = 1/bitrate;

df = 1/dt;

time = (-Nbits/2:Nbits/2-1)*dt;

FREQ=(-(Nbits/2)Nbits/2-1))*1/(dt*Nbits);

FREQGHz = FREQ.*1e-9;

%% %%%%%%%%%%%%%%%% Drive Current %%%%%%%%%%%%%%%%%%%%%%%%%%

% amp = 10e-3; %mA

%

% Bias = 50e-3; %mA

%%%%default%%%%

amp = 4e-3; %mA

Bias = 50e-3; %mA

%%%%%%%%%%%%%%%

I = (amp/2)*sin(2*pi*freq*time);

%Bias = Biasing*1e-3;

I = I + Bias;

%I = Bias.*ones(size(time));

% I = I - mean(I);

%% %%%%%%%%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%

g0 = 1e-12; % Differential gain coeff. !!!!1.6e-12

Nom = 1.4e23; % Transparency density (m^-3) 1.4e23

V = 11e-17; % Volume of active layer 11e-17

tp = 2e-11; % Photon lifetime (s) !!!!2.4e-11

tn = .3e-9; % Carrier recombination lifetime !!!!.03e-8

OCF = 0.35; % Mode confinement factor 0.35

B = 0.0000; % Beta, spontaneous emission factor 0.0

q = 1.6e-19; % Charge of electron (C) 1.6e-19

alpha = 6.8; % Linewidth Enhancement Factor

%Ac = 8e-3; % Amplitude of carriers

Ac = amp;

%I_bias = 70e-3; % Bias Current

I_bias = Bias;

deltaf = -01e9; % Frequency Detuning -11

deltaw = 2*pi*deltaf; %Convert frequency to radians

Si = 0e20; %Injection Level (Injected photon density) - change to 0.7e20

Kc = 2.5e11; %Injected light coulping coefficient

if Si == 0;

deltaw = 0;%it there's no injection then ignore the detuning in the phase equation

end

h = 6.625e-34; % Plancks constant

R = 0.32; % Reflectivity in cavity

c = 3e8; % Speed of light

n = 3.63; % Refractive index

Ar = .03e-12; % Area of the active region

lamda = 1550e-9; % Wavelength of output light

f_laser = c/lamda;% Frequency of the output light

responsivity = 0.6; % Responsivity of the detector

fm1 = 1e7; %Data Modulation Rate

bitperiod = 1/fm1; %Data bit period Period = 1/f

numcarriers = 5; %number of subcarriers

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% THIS SECTION GETS THE STEADY STATE VALUES FOR PHASE PHOTON DENSITY AND CARRIER % % % % DENSITY IT THEN USES THEM TO WORK OUT THE SMALL SIGNAL RESPONSE OF THE LASER WHICH % % IT PLOTS THIS IS THE MODULATION RESPONSE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the coefficients for a quartic equation in N0 to determine the steady

% state carrier density

% a_var*x^4 + b_var*x^3 + c_var*x^2 + d_var*x + e_var = 0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a_var = OCF^2*V^2*q^2*g0^2*tp^2*(alpha^2 + B^2 - 2*B + 1);

b_var = OCF*V*q*g0*tp*(I_bias*OCF*g0*tn*tp*(alpha^2 - B + 1) ...

+ V*q*(Nom*OCF*g0*tp*(alpha^2 + (B-1)^2) + alpha^2 ...

+ alpha*2*tp*deltaw -B + 1));

c_var = (I_bias^2*OCF^2*g0^2*tn^2*tp^2*(alpha^2 + 1) ...

+ 2*I_bias*OCF*V*q*g0*tn*tp*(2*Nom*OCF*g0*tp*(alpha^2 - B + 1) + 2*alpha^2 ...

+ 4*alpha*tp*deltaw - B + 2) + V^2*q^2*(4*Kc^2*Si*g0*tn*tp^2 ...

+ Nom^2*OCF^2*g0^2*tp^2*(alpha^2 + B^2 - 2*B + 1) ...

+ 2*Nom*OCF*g0*tp*(alpha^2 + 2*alpha*tp*deltaw - B + 1) + alpha^2 ...

+ 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1));

d_var = tn*(I_bias^2*OCF*g0*tn*tp*(Nom*OCF*g0*tp*(alpha^2 + 1) + alpha^2 ...

+ 2*alpha*tp*deltaw + 1) + I_bias*V*q*(2*Kc^2*Si*g0*tn*tp^2 ...

+ Nom^2*OCF^2*g0^2*tp^2*(alpha^2 - B + 1) ...

+ Nom*OCF*g0*tp*(2*alpha^2 + 4*alpha*tp*deltaw - B + 2) ...

+ alpha^2 + 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1) ...

+ 2*Kc^2*Nom*Si*V^2*q^2*g0*tp^2);

e_var = I_bias^2*tn^2*(Nom^2*OCF^2*g0^2*tp^2*(alpha^2 + 1) ...

+ 2*Nom*OCF*g0*tp*(alpha^2 + 2*alpha*tp*deltaw +1) + alpha^2 ...

+ 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1) ...

+ 4*I_bias*Kc^2*Nom*Si*V*q*g0*tn^2*tp^2;

%Get the roots of the above equation

N0_eq = [a_var -2*b_var c_var -2*d_var e_var];

N0_Roots = roots(N0_eq);

N0 = N0_Roots(4); %N0, the steady state carrier density 4!!

% S0 the steady state photon density is worked out by rearranging the

% carrier density rate equation

S0 = ((I_bias*tn)- (N0*q*V))/(tn*g0*(N0-Nom)*q*V);

% Injection Ratio

ratio = Si/S0;

% Output power is worked out using Photon Density

Pout = S0.*c./(2*OCF*n)*h*f_laser*Ar*(1-R);

%To work out steady state phase, Phi0, both equations give same answer but different %sign

if Si == 0;

Phi0 = 0;

else

% Phi0 = acos(((S0/tp) - (OCF*g0*(N0-Nom)*S0))/(2*Kc*sqrt(S0*Si)));

Phi0 = asin((-alpha/(2*tp) - (deltaw) + ...

((alpha/2)*OCF*g0*(N0-Nom)))/(Kc*sqrt(Si/S0)));

end

% Set up frequency points to plot the frequency response

% freqM = [10e6:10e6:20e9];

% wm = 2*pi*freqM;

%

% % %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% % % Using the small signal derivation of the rate equations to get the modulation %response Modulation response is given as change in photon number with respect to % %change in input current Input current is Ac and change in photon density is s1

% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% X = 2*Kc*sqrt(Si*S0)*cos(Phi0);

% Y = Kc*sqrt(Si/S0)*sin(Phi0);

% O = (j.*wm + (X/(2*S0)));

% P = (j.*wm + (1/tn) + g0*S0);

% if Si == 0

% s1 = -((Ac/(q*V))*g0*S0*OCF)./(wm.^2 - (wm).*(i*g0*S0) + (wm).*(i/(tn)) - (g0*S0)/tp);

% end

% if Si > 0

% s1 = (((OCF*Ac*g0*S0)./(P.*q*V)).*(1 - Y*alpha./O)) ... %top line

% ./((j.*wm) + (X/(2*S0)) + (Y^2./O) + (X*Y*alpha*g0./(P.*O)) - ...

% (S0*Y*alpha*g0./(P.*O*tp)) - (X*g0./P) + (g0*S0./(P.*tp)));

% end

%% %%%%%%% Plot Resonance Freq. %%%%%%%%%

% figure(2);

% s = abs(s1); %Get absolute values of the change in photon density

% logs = 20*log10(s/Ac); %Get the log of change in photon number with respect to change

% %in input current (The modulation response)

% norms = logs - logs(1); %Normalise it to the first value

% plot(freq.*1e-9,norms); % and plot

% axis([0 20 -30 40]);

% title('Modulation Response');

% xlabel('Frequency (GHz)'); ylabel('Response (dB)');

% hold on; grid;

%% %%%%%%%%%%%%%%%%%%% Initial Parameters %%%%%%%%%%%%%%%%%%%%%

S0 = abs(S0);

N0 = abs(N0);

Phi0 = abs(Phi0);

P1 = S0;

N1 = N0;

Phi1 = Phi0;

% P1 = 1e23;

% N1 = Nom;

% Phi1 = -3;

%% %%%%%%%%%%%%%%%%%%%%%% Calculate Rate Equations %%%%%%%%%%%%%%%%%

dt=time(2)-time(1);

N=N1*ones(size(time));

P=P1*ones(size(time));

Phi=Phi1*ones(size(time));

for k=2:length(time),

dPdt= ((OCF*g0*(N(k-1)-Nom) - (1/tp))*P(k-1)) + ((OCF*B*N(k-1))/tn)...

+ ((2*Kc)*sqrt(Si*P(k-1))*cos(Phi(k-1)));

dNdt= (I(k-1)/(q*V)) - (g0*(N(k-1)-Nom)*P(k-1)) - (N(k-1)/tn);

P(k)=dPdt*dt+P(k-1); %Photon Density

N(k)=dNdt*dt+N(k-1); %Carrier Density

if P(k) <= 0

P(k) = 100;

end

if Si == 0

Phi(k) = 0;

else

dPhidt = (alpha/2)*(OCF*g0*(N(k-1)-Nom) - 1/tp) - (deltaw) - (Kc*sqrt(Si/P(k-1))*sin(Phi(k-1)));

Phi(k)=dPhidt*dt+Phi(k-1);

end

end

%% %%%%%%%%%%%%%%%%%%%% Get Power %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Power_leftside = P.*c/(2*OCF*n)*h*f_laser*Ar*(1-R);

Power_leftside = Power_leftside*1e3; %mW

%Power_leftside = Power_leftside - mean(Power_leftside);

Power_leftsidef=abs(fftshift(fft(fftshift(Power_leftside))));

%Note: Taking FFT includes transcients at start => high power

%Normal for frequency response graph to start at 0dBm for small

%signal analysis

% Sig_Pwr = max(Power_leftside(1,(length(Power_leftside))/2:length(Power_leftside)));

% response(ind2) = Sig_Pwr; %mW

Sig_Pwr = max(Power_leftside(1,(length(Power_leftside)-5000):length(Power_leftside)))...

- min(Power_leftside(1,(length(Power_leftside)-5000):length(Power_leftside))); % Power signal amplitude after transcients mW

response(ind2) = 10*log10(Sig_Pwr); %dBm

%response(ind2) = 10*log10(Sig_Pwr./(amp.^2))

Sig_Count(ind2) = Sig_Pwr;

ind2 = ind2 + 1;

end

% Mult_Res(Bind, = response;

% Bind = Bind + 1;

% ind2 = 1;

% end

%% %%%%%%%%%%%%%%%%%%%%%% Graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure(3)

% plot(FREQGHz, 10*log10(Power_leftsidef/max(Power_leftsidef)))

% title('Power v Frequency')

% xlabel('Frequency (GHz))')

% ylabel('Power(dB)')

% axis([min(FREQGHz) max(FREQGHz) -60 0])

response_graph = steps;

figure, plot(response_graph, response)

title('Laser Frequency response')

xlabel('Frequency (GHz)')

ylabel('Power (dBm)')

hold on; grid;

% figure(6)

% plot(10*log10(abs(fftshift(fft(fftshift(P(length(P)/2:length(P))))))))

%

%

% ax=[0 max(FREQGHz) -60 0];

% If=abs(fftshift(fft(fftshift(I)))).^2;

% figure(1)

% plot (time,I.*1e3)

% title('Tx time signal')

% xlabel('time')

% ylabel('mA')

%

%

% figure(2)

% plot(FREQGHz,10*log10(If/max(If)))

% title('Tx spectrum')

% axis(ax)

I am designing a optical comb to be used in DWDM systems, can anyone please help me in giving some conceptual inputs regarding how the laser can be gain switched to achieve the comb. I have a code that solves the rate equations to give photon density, carrier density and phase of the light emmitted. So can you please throw some light on this topic as I'm a bit new to this domain.

Thanking you in advance.

Sur76

P.S: The following is the code I have.

clear all

ind2 = 1;

steps = 0:0.1:20;

% Bind = 1;

% Bias_step = 16:4:24;

%

% for Biasing = Bias_step

for ind = steps

%% %%%%%%%%%%%%%%%%%%%%% Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nbits = 1024*10;

freq = ind*1e9;

%freq = 0.1*1e9;

w = 2*pi*freq;

% bitrate = 9e10;

bitrate = 300e9;

dt = 1/bitrate;

df = 1/dt;

time = (-Nbits/2:Nbits/2-1)*dt;

FREQ=(-(Nbits/2)Nbits/2-1))*1/(dt*Nbits);

FREQGHz = FREQ.*1e-9;

%% %%%%%%%%%%%%%%%% Drive Current %%%%%%%%%%%%%%%%%%%%%%%%%%

% amp = 10e-3; %mA

%

% Bias = 50e-3; %mA

%%%%default%%%%

amp = 4e-3; %mA

Bias = 50e-3; %mA

%%%%%%%%%%%%%%%

I = (amp/2)*sin(2*pi*freq*time);

%Bias = Biasing*1e-3;

I = I + Bias;

%I = Bias.*ones(size(time));

% I = I - mean(I);

%% %%%%%%%%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%

g0 = 1e-12; % Differential gain coeff. !!!!1.6e-12

Nom = 1.4e23; % Transparency density (m^-3) 1.4e23

V = 11e-17; % Volume of active layer 11e-17

tp = 2e-11; % Photon lifetime (s) !!!!2.4e-11

tn = .3e-9; % Carrier recombination lifetime !!!!.03e-8

OCF = 0.35; % Mode confinement factor 0.35

B = 0.0000; % Beta, spontaneous emission factor 0.0

q = 1.6e-19; % Charge of electron (C) 1.6e-19

alpha = 6.8; % Linewidth Enhancement Factor

%Ac = 8e-3; % Amplitude of carriers

Ac = amp;

%I_bias = 70e-3; % Bias Current

I_bias = Bias;

deltaf = -01e9; % Frequency Detuning -11

deltaw = 2*pi*deltaf; %Convert frequency to radians

Si = 0e20; %Injection Level (Injected photon density) - change to 0.7e20

Kc = 2.5e11; %Injected light coulping coefficient

if Si == 0;

deltaw = 0;%it there's no injection then ignore the detuning in the phase equation

end

h = 6.625e-34; % Plancks constant

R = 0.32; % Reflectivity in cavity

c = 3e8; % Speed of light

n = 3.63; % Refractive index

Ar = .03e-12; % Area of the active region

lamda = 1550e-9; % Wavelength of output light

f_laser = c/lamda;% Frequency of the output light

responsivity = 0.6; % Responsivity of the detector

fm1 = 1e7; %Data Modulation Rate

bitperiod = 1/fm1; %Data bit period Period = 1/f

numcarriers = 5; %number of subcarriers

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% THIS SECTION GETS THE STEADY STATE VALUES FOR PHASE PHOTON DENSITY AND CARRIER % % % % DENSITY IT THEN USES THEM TO WORK OUT THE SMALL SIGNAL RESPONSE OF THE LASER WHICH % % IT PLOTS THIS IS THE MODULATION RESPONSE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the coefficients for a quartic equation in N0 to determine the steady

% state carrier density

% a_var*x^4 + b_var*x^3 + c_var*x^2 + d_var*x + e_var = 0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a_var = OCF^2*V^2*q^2*g0^2*tp^2*(alpha^2 + B^2 - 2*B + 1);

b_var = OCF*V*q*g0*tp*(I_bias*OCF*g0*tn*tp*(alpha^2 - B + 1) ...

+ V*q*(Nom*OCF*g0*tp*(alpha^2 + (B-1)^2) + alpha^2 ...

+ alpha*2*tp*deltaw -B + 1));

c_var = (I_bias^2*OCF^2*g0^2*tn^2*tp^2*(alpha^2 + 1) ...

+ 2*I_bias*OCF*V*q*g0*tn*tp*(2*Nom*OCF*g0*tp*(alpha^2 - B + 1) + 2*alpha^2 ...

+ 4*alpha*tp*deltaw - B + 2) + V^2*q^2*(4*Kc^2*Si*g0*tn*tp^2 ...

+ Nom^2*OCF^2*g0^2*tp^2*(alpha^2 + B^2 - 2*B + 1) ...

+ 2*Nom*OCF*g0*tp*(alpha^2 + 2*alpha*tp*deltaw - B + 1) + alpha^2 ...

+ 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1));

d_var = tn*(I_bias^2*OCF*g0*tn*tp*(Nom*OCF*g0*tp*(alpha^2 + 1) + alpha^2 ...

+ 2*alpha*tp*deltaw + 1) + I_bias*V*q*(2*Kc^2*Si*g0*tn*tp^2 ...

+ Nom^2*OCF^2*g0^2*tp^2*(alpha^2 - B + 1) ...

+ Nom*OCF*g0*tp*(2*alpha^2 + 4*alpha*tp*deltaw - B + 2) ...

+ alpha^2 + 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1) ...

+ 2*Kc^2*Nom*Si*V^2*q^2*g0*tp^2);

e_var = I_bias^2*tn^2*(Nom^2*OCF^2*g0^2*tp^2*(alpha^2 + 1) ...

+ 2*Nom*OCF*g0*tp*(alpha^2 + 2*alpha*tp*deltaw +1) + alpha^2 ...

+ 4*alpha*tp*deltaw + 4*deltaw^2*tp^2 + 1) ...

+ 4*I_bias*Kc^2*Nom*Si*V*q*g0*tn^2*tp^2;

%Get the roots of the above equation

N0_eq = [a_var -2*b_var c_var -2*d_var e_var];

N0_Roots = roots(N0_eq);

N0 = N0_Roots(4); %N0, the steady state carrier density 4!!

% S0 the steady state photon density is worked out by rearranging the

% carrier density rate equation

S0 = ((I_bias*tn)- (N0*q*V))/(tn*g0*(N0-Nom)*q*V);

% Injection Ratio

ratio = Si/S0;

% Output power is worked out using Photon Density

Pout = S0.*c./(2*OCF*n)*h*f_laser*Ar*(1-R);

%To work out steady state phase, Phi0, both equations give same answer but different %sign

if Si == 0;

Phi0 = 0;

else

% Phi0 = acos(((S0/tp) - (OCF*g0*(N0-Nom)*S0))/(2*Kc*sqrt(S0*Si)));

Phi0 = asin((-alpha/(2*tp) - (deltaw) + ...

((alpha/2)*OCF*g0*(N0-Nom)))/(Kc*sqrt(Si/S0)));

end

% Set up frequency points to plot the frequency response

% freqM = [10e6:10e6:20e9];

% wm = 2*pi*freqM;

%

% % %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% % % Using the small signal derivation of the rate equations to get the modulation %response Modulation response is given as change in photon number with respect to % %change in input current Input current is Ac and change in photon density is s1

% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% X = 2*Kc*sqrt(Si*S0)*cos(Phi0);

% Y = Kc*sqrt(Si/S0)*sin(Phi0);

% O = (j.*wm + (X/(2*S0)));

% P = (j.*wm + (1/tn) + g0*S0);

% if Si == 0

% s1 = -((Ac/(q*V))*g0*S0*OCF)./(wm.^2 - (wm).*(i*g0*S0) + (wm).*(i/(tn)) - (g0*S0)/tp);

% end

% if Si > 0

% s1 = (((OCF*Ac*g0*S0)./(P.*q*V)).*(1 - Y*alpha./O)) ... %top line

% ./((j.*wm) + (X/(2*S0)) + (Y^2./O) + (X*Y*alpha*g0./(P.*O)) - ...

% (S0*Y*alpha*g0./(P.*O*tp)) - (X*g0./P) + (g0*S0./(P.*tp)));

% end

%% %%%%%%% Plot Resonance Freq. %%%%%%%%%

% figure(2);

% s = abs(s1); %Get absolute values of the change in photon density

% logs = 20*log10(s/Ac); %Get the log of change in photon number with respect to change

% %in input current (The modulation response)

% norms = logs - logs(1); %Normalise it to the first value

% plot(freq.*1e-9,norms); % and plot

% axis([0 20 -30 40]);

% title('Modulation Response');

% xlabel('Frequency (GHz)'); ylabel('Response (dB)');

% hold on; grid;

%% %%%%%%%%%%%%%%%%%%% Initial Parameters %%%%%%%%%%%%%%%%%%%%%

S0 = abs(S0);

N0 = abs(N0);

Phi0 = abs(Phi0);

P1 = S0;

N1 = N0;

Phi1 = Phi0;

% P1 = 1e23;

% N1 = Nom;

% Phi1 = -3;

%% %%%%%%%%%%%%%%%%%%%%%% Calculate Rate Equations %%%%%%%%%%%%%%%%%

dt=time(2)-time(1);

N=N1*ones(size(time));

P=P1*ones(size(time));

Phi=Phi1*ones(size(time));

for k=2:length(time),

dPdt= ((OCF*g0*(N(k-1)-Nom) - (1/tp))*P(k-1)) + ((OCF*B*N(k-1))/tn)...

+ ((2*Kc)*sqrt(Si*P(k-1))*cos(Phi(k-1)));

dNdt= (I(k-1)/(q*V)) - (g0*(N(k-1)-Nom)*P(k-1)) - (N(k-1)/tn);

P(k)=dPdt*dt+P(k-1); %Photon Density

N(k)=dNdt*dt+N(k-1); %Carrier Density

if P(k) <= 0

P(k) = 100;

end

if Si == 0

Phi(k) = 0;

else

dPhidt = (alpha/2)*(OCF*g0*(N(k-1)-Nom) - 1/tp) - (deltaw) - (Kc*sqrt(Si/P(k-1))*sin(Phi(k-1)));

Phi(k)=dPhidt*dt+Phi(k-1);

end

end

%% %%%%%%%%%%%%%%%%%%%% Get Power %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Power_leftside = P.*c/(2*OCF*n)*h*f_laser*Ar*(1-R);

Power_leftside = Power_leftside*1e3; %mW

%Power_leftside = Power_leftside - mean(Power_leftside);

Power_leftsidef=abs(fftshift(fft(fftshift(Power_leftside))));

%Note: Taking FFT includes transcients at start => high power

%Normal for frequency response graph to start at 0dBm for small

%signal analysis

% Sig_Pwr = max(Power_leftside(1,(length(Power_leftside))/2:length(Power_leftside)));

% response(ind2) = Sig_Pwr; %mW

Sig_Pwr = max(Power_leftside(1,(length(Power_leftside)-5000):length(Power_leftside)))...

- min(Power_leftside(1,(length(Power_leftside)-5000):length(Power_leftside))); % Power signal amplitude after transcients mW

response(ind2) = 10*log10(Sig_Pwr); %dBm

%response(ind2) = 10*log10(Sig_Pwr./(amp.^2))

Sig_Count(ind2) = Sig_Pwr;

ind2 = ind2 + 1;

end

% Mult_Res(Bind, = response;

% Bind = Bind + 1;

% ind2 = 1;

% end

%% %%%%%%%%%%%%%%%%%%%%%% Graphs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure(3)

% plot(FREQGHz, 10*log10(Power_leftsidef/max(Power_leftsidef)))

% title('Power v Frequency')

% xlabel('Frequency (GHz))')

% ylabel('Power(dB)')

% axis([min(FREQGHz) max(FREQGHz) -60 0])

response_graph = steps;

figure, plot(response_graph, response)

title('Laser Frequency response')

xlabel('Frequency (GHz)')

ylabel('Power (dBm)')

hold on; grid;

% figure(6)

% plot(10*log10(abs(fftshift(fft(fftshift(P(length(P)/2:length(P))))))))

%

%

% ax=[0 max(FREQGHz) -60 0];

% If=abs(fftshift(fft(fftshift(I)))).^2;

% figure(1)

% plot (time,I.*1e3)

% title('Tx time signal')

% xlabel('time')

% ylabel('mA')

%

%

% figure(2)

% plot(FREQGHz,10*log10(If/max(If)))

% title('Tx spectrum')

% axis(ax)

Last edited: