Frequency Comb generation in Matlab

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)