lmaamp_back_wave.m

clear
close('all')
%%Rate media rate equations, quasi 2-level system

%for ASE see Pedersen, B., The design of erbium-doped fiber amplifiers. 1991, Journal of Lightwave Tech. 

%-----plotting
% numPlots = 11;
% 
% m = floor(sqrt(numPlots));      %rows of plots
% n = ceil(numPlots/m);           %columns of plots
% p = 1;                          %current plotting position
% 
% scrsz = get(0,'ScreenSize');
% f = figure('Position',[1 scrsz(4) 2*scrsz(3)/3 2*scrsz(4)/3], 'Visible', 'off');

%color parameters
nCol = 12;
pcol = [1 0 0];
sigcol = [1 0 1];
asefcol = [0 0 1];
asebcol = [0 1 1];
%-----


%constants
h = 6.62606957E-34;     %J*s
c = 299792458;          %m/s

%parameters
%crystal parameters
L = 1.7;                  %m fiber length
%N = 6.1815E25; %30.1*10^24;  %number density of doopant atoms, #/m^3
abs_p = 6.0/(10/log(10));
alpha_p = 0.1;          %additional pump loss
alpha_s = 1/(10/log(10));

%cross sections, from Nufern file for Gen VIII fiber.
%s_ap = 1.7806E-24;     %absorption pump, m^2
%s_ep = 1.7623E-24;     %emission pump, m^2
%s_as = 0.0562E-24;     %abs signal, m^2
%s_es = 0.7634E-24;     %emi signal, m^2

%create wavelength array for ase
%ase_lim = [1010,1040]*1E-9;
%lambda_ase = linspace(ase_lim(1), ase_lim(2),300);

%wavelengths
l_p = 0.976E-6;         %pump wavelength, m
l_s = 1.030E-6;         %signal wavelength, m
v_p = c/l_p;            %pump freq, Hz
v_s = c/l_s;            %signal freq, Hz

%import cross sections from file
[folder, name, ext] = fileparts(which('lmaamp_back_wave.m'));
M = csvread(fullfile(folder,'Nufern VIII.csv'),12,0);
lambda_ase = M(:,1)'*1E-9;      %column 1, wavelength, nm
s_a_ase = M(:,2)'*1E-24;        %column 2, abs cs, pm^2
s_e_ase = M(:,3)'*1E-24;        %column 3, ems cs, pm^2

dl_ase = gradient(lambda_ase);
nu_ase = c./lambda_ase;
dnu_ase = gradient(nu_ase);

s_a = interp1(lambda_ase,s_a_ase, [l_s,l_p]);
s_e = interp1(lambda_ase,s_e_ase, [l_s,l_p]);
s_as = s_a(1);
s_ap = s_a(2);
s_es = s_e(1);
s_ep = s_e(2);

ref_index = 1.447;
NA_core = 0.06;

%calculate spantaneous emission lifetime from cross section
tau_se = 1/((8*pi*ref_index^2/c^2)*trapz(-nu_ase,nu_ase.^2.*s_e_ase));

%calculate crosssection based on fits from nufern
% ase_em_coefs = [-0.0009083478128*1E18,1.86299326*1E9,-954.4494928];
% ase_ab_coefs = [-0.00002194150379*1E18,0.04204067705*1E9,-19.96593069];
% s_e_ase = polyval(ase_em_coefs,lambda_ase)*1E-24;
% s_a_ase = polyval(ase_ab_coefs,lambda_ase)*1E-24;
%tau_se = 770E-6;        %spontaneous emission lifetime, s, see Barnard 1994 j quant elec.

dCore = 30E-6;           %core diameter (m)
dClad = 250E-6;          %clad diameter
MFD = 21E-6;
MFA = pi*(MFD/2)^2;    %dopant mode field area, ~estimate from NA
Ap = pi*(dClad/2)^2;    %pump (clad) area
As = pi*(dCore/2)^2;    %signal (core) area

Gp = MFA/Ap;    %Pump overlap
Gs = MFA/As;    %signal overlap

N = (abs_p/s_ap)/(As/Ap);

%calculated constants (see notebook 1, page , or equation below)
b_p = (s_ap + s_ep)/(h*v_p);
b_s = (s_as + s_es)/(h*v_s);
b_ase = (s_a_ase + s_e_ase)./(h*nu_ase);
a_p = s_ap/(h*v_p);
a_s = s_as/(h*v_s);
a_ase = s_a_ase./(h*nu_ase);


%With these defined contsants, rate equation looks like
%dn2/dt = a_p*I_p + a_s*I_s - n2*(b_p*I_p + b_s*I_s + 1/tau_se)



%%
%Position-dependent pump intensity
dz = 0.002;                %spacial slice, 2mm
z = 0:dz:L;

P_pi = 5;        %lowest pump power (W)
P_pf = 25;        %highest pump power (W)

%s_ap = s_ap*(dCore/dClad)^2;
%s_ep = s_ep*(dCore/dClad)^2;

FSeed = 500E3; %in rep rate s^-1
ESeed = 80E-9; % in J

I_0p = (P_pi:P_pi:P_pf)/Ap;          %Initial pump intensity
I_0s = (ESeed*FSeed)/As;                %total seed signal intensity

%dv_ase = 53E-9*(v_s/l_s);

%Initialize parameters
num_I = length(I_0p);
num_z = length(z);
num_l = length(lambda_ase);

I_p = zeros(num_I,num_z);        %pump intensity
I_sig = zeros(num_I,num_z);      %singal (of interest) intensity
I_asef = zeros(num_I,num_l,num_z);     %forward ase intensity, also initial cond.
I_aseb = zeros(num_I,num_l,num_z);     %backward ase intensity, also initial cond.
n_2 = zeros(num_I,num_z);        %excited state population
dg = zeros(num_I,num_z);         %gain contribution
g = zeros(num_I,num_z);          %gain coefficient
G = zeros(num_I,num_z);          %gain  

%arrays for cross sections
s_e_ase_mat = repmat(s_e_ase, [num_I,1]);
s_a_ase_mat = repmat(s_a_ase, [num_I,1]);
lambda_ase_mat = repmat(lambda_ase, [num_I,1]);
dl_ase_mat = repmat(dl_ase, [num_I,1]);

%initial Pump/Signal conditions
%backpumped fiber
I_p(:,end) = I_0p;
I_sig(:,1) = I_0s;

%DE's for light intensity
dI_p = @(z,I_p,n2) (-Gp*(s_ap*N*(1-n2) - s_ep*N*n2) - alpha_p).*I_p;
dI_sig = @(z,I_sig,n2) (-Gs*(s_as*N*(1-n2) - s_es*N*n2) - alpha_s).*I_sig;
dI_ase = @(z,I_ase, n2) (-Gs.*(s_a_ase_mat.*N.*(1-n2) - s_e_ase_mat.*N.*n2)).*I_ase + 2.*Gs.*n2.*h.*c^2.*N.*s_e_ase_mat.*dl_ase_mat./(As*lambda_ase_mat.^3);
%note ASE has same DE, but z directions are different, this is implemented
%below


%while-loop determinants
j = 0;
cGain = zeros(num_I,1);
errG = 1;
err = [];
loop_max = 100;
gainerr_max = 1E-3;

%Main loop
wb = waitbar(0, 'Performing iterations...');

while j < loop_max && abs(errG) > gainerr_max
    pGain = cGain;
    
    %RK4 method
    %single-pass propagation
    for i = 1:num_z
        k = num_z - i + 1;
        
        %I_s(:,i) = I_sig(:,i) + sum(I_asef(:,i,:),3) + sum(I_aseb(:,i,:),3);
        %n_2(:,i) = (a_p*I_p(:,i) + a_s*I_s(:,i))./(b_p*I_p(:,i) + b_s*I_s(:,i) + 1/tau_se);
        
        n_2(:,i) = ((Gp*Ap)*a_p*I_p(:,i) + (Gs*As)*a_s*I_sig(:,i) + (Gs*As)*sum(repmat(a_ase,[num_I,1]).*(I_asef(:,:,i)+I_aseb(:,:,i)),2))./...
           ((Gp*Ap)*b_p*I_p(:,i) + (Gs*As)*b_s*I_sig(:,i) + (Gs*As)*sum(repmat(b_ase,[num_I,1]).*(I_asef(:,:,i)+I_aseb(:,:,i)),2) + As/tau_se);
        
        n_2(:,k) = ((Gp*Ap)*a_p*I_p(:,k) + (Gs*As)*a_s*I_sig(:,k) + (Gs*As)*sum(repmat(a_ase,[num_I,1]).*(I_asef(:,:,k)+I_aseb(:,:,k)),2))./...
            ((Gp*Ap)*b_p*I_p(:,k) + (Gs*As)*b_s*I_sig(:,k) + (Gs*As)*sum(repmat(b_ase,[num_I,1]).*(I_asef(:,:,k)+I_aseb(:,:,k)),2) + As/tau_se);
        
        if i<num_z
            
            %update I_p
            k1 = dI_p(z(k), I_p(:,k), n_2(:,k));
            k2 = dI_p(z(k) + dz/2, I_p(:,k) + k1*dz/2, n_2(:,k));
            k3 = dI_p(z(k) + dz/2, I_p(:,k) + k2*dz/2, n_2(:,k));
            k4 = dI_p(z(k) + dz, I_p(:,k) + k3*dz, n_2(:,k));
            I_p(:,k-1) = I_p(:,k) + (k1+2*k2+2*k3+k4)*dz/6;
            
            %update I_sig
            k1 = dI_sig(z(i), I_sig(:,i), n_2(:,i));
            k2 = dI_sig(z(i) + dz/2, I_sig(:,i) + k1*dz/2, n_2(:,i));
            k3 = dI_sig(z(i) + dz/2, I_sig(:,i) + k2*dz/2, n_2(:,i));
            k4 = dI_sig(z(i) + dz, I_sig(:,i) + k3*dz, n_2(:,i));
            I_sig(:,i+1) = I_sig(:,i) + (k1+2*k2+2*k3+k4)*dz/6;
            
            %update I_asef
            k1 = dI_ase(z(i), (I_asef(:,:,i)), repmat(n_2(:,i),[1,num_l]));
            k2 = dI_ase(z(i) + dz/2, (I_asef(:,:,i)) + k1*dz/2, repmat(n_2(:,i),[1,num_l]));
            k3 = dI_ase(z(i) + dz/2, (I_asef(:,:,i)) + k2*dz/2, repmat(n_2(:,i),[1,num_l]));
            k4 = dI_ase(z(i) + dz, (I_asef(:,:,i)) + k3*dz, repmat(n_2(:,i),[1,num_l]));
            I_asef(:,:,i+1) = (I_asef(:,:,i)) + (k1+2*k2+2*k3+k4)*dz/6;
            
            %update I_aseb
            k1 = dI_ase(z(k), (I_aseb(:,:,k)), repmat(n_2(:,k),[1,num_l]));
            k2 = dI_ase(z(k) + dz/2, (I_aseb(:,:,k)) + k1*dz/2, repmat(n_2(:,k),[1,num_l]));
            k3 = dI_ase(z(k) + dz/2, (I_aseb(:,:,k)) + k2*dz/2, repmat(n_2(:,k),[1,num_l]));
            k4 = dI_ase(z(k) + dz, (I_aseb(:,:,k)) + k3*dz,repmat(n_2(:,k),[1,num_l]));
            I_aseb(:,:,k-1) = (I_aseb(:,:,k)) + (k1+2*k2+2*k3+k4)*dz/6;
            
            
        end
        
        %single-pass gain
%         g(:,k) = (s_es*N*n_2(:,k) - s_as*N*(1-n_2(:,k)));
%         dg(:,k) = g(:,k)*dz;
%         G(:,k) = exp(sum(dg(:,i:k),2));
%         cGain = G(:,end);

        
    end
    
    g = Gs*(s_es*N*n_2 - s_as*N*(1-n_2));
    dg = g*dz;
    G = exp(sum(dg,2));
    cGain = G;
    
    errG = max(abs((cGain - pGain)./cGain));
    err = [err, errG];
    
    j = j+1;
    
    waitbar(j/loop_max, wb,sprintf('Performing iterations... %d; Error = %4.3e', [j,errG]))
    
    if abs(errG) > 1E100
        display('Error exceeds maximum; simulation terminated');
        break
    end
    
end
close(wb)
%display final loop count and err
display(['Number of loops: ',num2str(j)]);
display(['Gain error: ', num2str(errG)]);


%%
%Plotting
%-----plotting
numPlots = 6;

m = floor(sqrt(numPlots));      %rows of plots
n = ceil(numPlots/m);           %columns of plots
p = 1;                          %current plotting position

scrsz = get(0,'ScreenSize');
f = figure('Position',[1 scrsz(4) 2*scrsz(3)/3 2*scrsz(4)/3], 'Visible', 'off');

figure(f)
colors = hsv(num_I);
set(gcf, 'Colormap', colors);


%Pump plot
subplot(m,n,p)
p=p+1;
set(gca, 'ColorOrder', colors);
hold on
plot(z, I_p*Ap)
title('Pump intensity')
xlabel('Position along fiber (m)')
ylabel('Pump Power (W)')
%annotation('textbox', [0.6,0.6,0.1,0.1],...
%           'String', ['I_0 = ' num2str(I_p(1)/(1E7)) ' kW/cm^2']);
hold off


%Signal plot
subplot(m,n,p)
p=p+1;
set(gca, 'ColorOrder', colors);
hold on
plot(z, I_sig*As)
title('Signal intensity')
xlabel('Position along fiber (m)')
ylabel('Signal Power (W)')
%annotation('textbox', [0.6,0.6,0.1,0.1],...
%           'String', ['I_0 = ' num2str(I_p(1)/(1E7)) ' kW/cm^2']);
hold off


%ASE plot
subplot(m,n,p)
p=p+1;
set(gca, 'ColorOrder', colors);
hold on
plot(z, As*squeeze(sum(I_asef,2)), '-')
plot(z, As*squeeze(sum(I_aseb,2)), '-.')
title('ASE Signal intensity')
xlabel('Position along fiber (m)')
ylabel('ASE Signal Power (W)')
%annotation('textbox', [0.6,0.6,0.1,0.1],...
%           'String', ['I_0 = ' num2str(I_p(1)/(1E7)) ' kW/cm^2']);
hold off

%Inversion plot
subplot(m,n,p)
p=p+1;
set(gca, 'ColorOrder', colors);
hold on
plot(z, n_2)
title('Inversion')
xlabel('Position along fiber (m)')
ylabel('Upper state population')
%annotation('textbox', [0.6,0.6,0.1,0.1],...
%           'String', ['I_0 = ' num2str(I_p(1)/(1E7)) ' kW/cm^2']);
hold off

%Gain plot
subplot(m,n,p)
p=p+1;
set(gca, 'ColorOrder', colors);
hold on
plot(z, g)
title('Gain')
xlabel('Position along fiber (m)')
ylabel('Gain Coef')
%annotation('textbox', [0.6,0.6,0.1,0.1],...
%           'String', ['I_0 = ' num2str(I_p(1)/(1E7)) ' kW/cm^2']);
hold off

%Output plot
subplot(m,n,p)
p=p+1;
set(gca, 'ColorOrder', colors);
hold on
gaindB = 10*log10(I_sig(:,end)/I_0s);
plotyy(I_0p*Ap, I_sig(:,end)*As,I_0p*Ap, gaindB)
title('Pump Efficiency')
xlabel('Pump power (W)')
ylabel('Output power (W)')


%annotation('textbox', [0.6,0.6,0.1,0.1],...
%           'String', ['I_0 = ' num2str(I_p(1)/(1E7)) ' kW/cm^2']);
hold off
%
% %%
% %Analysis
% 
% d = 6E-6;
% %P_t = (I_p(:,end) + I_sig(:,end) + I_asef(:,end) + I_aseb(:,end))*pi*(d/2)^2;
% P_0p = I_0p*pi*(d/2)^2;
% 
% effPlot = figure('Name', 'Total Output Power');
% plot(P_0p,P_t)