hello, am unable to debug the following code, I need somebody to help me run it successfully.

function Barow() % SOA_vsPower - SOA simulations versus input power % Determines signal gain and approximate noise figure versus bias current % distributions using paper 'Wideband Semiconductor Optical Amplifier % Steady-State Numerical Model' by Michael J. Connelly, IEEE Journal % Quantum Electron, vol. 37, no. 3, pp.439-447, 2001. % Author: M. Connelly, University of Limerick, Ireland. % $Revision: 1.0. $Date: 23/5/2007

clear all clc

%*********************************************************

% universal constants

global c e k h hbar kT

c = 3e8; % speed light h = 6.63e-34; % Planck constant (J-s) k = 1.38e-23; % Boltzmann constant (J/K) m0 =9.11e-31; % electron rest mass (kg) e = 1.6e-19; % electron charge (C) meV = 1e-3*e; T = 300; % Temperature (K)

hbar = h/(2*pi); % planck constant/2pi kT = k*T;

% SOA parameters

global yAs Lc Lt d W confine Kg n1 n2 neq0 eta_in eta_out R1 R2 K0 K1 Arad Brad Anrad Bnrad Caug Dleak... me mhh mlh

yAs = 0.892; % molar fraction of Arsenide in the active region Lc = 600; % Central active region length (um) Lt = 100; % Tapered active region length (um) d = 0.4e-6; %Active region thickness (m) W = 0.4e-6; % Central active region width (m) confine = 0.45; % Optical confinement factor Kg = 0.9e-10; % Bandgap shrinkage coefficient (eVm) n1 = 3.22; % InGaAsP active region refractive index n2 = 3.167; % InP region refractive index dn1dn = -1.8e-26; % Differential of active region refractive index with respect to carrier density (m^-3) neq0 = 3.22; % Equivalent refractive index at zero carrier density (from (29)) dneqdn = n1*confine/(neq0*(2-confine))*dn1dn; % Differential of equivalent refractive index with respect to carrier density (equation (30)) (m^-3) eta_in = 3; % Input coupling loss (dB) eta_out = 3; % Output coupling loss (dB) R1 = 5e-5; % Input facet reflectivity R2 = 5e-5; % Output facet reflectivity K0 = 6200; % Carrier independent absorption loss coefficient (m^-1) K1 = 7500; % Carrier dependent absorption loss coefficient (m^2) Arad = 1e7; % Linear radiative recombination coefficient (s^-1) Brad = 5.6e-16; % Bimolecular radiative recombination coefficient (m^3 s^-1) Anrad = 3.5e8; % Linear nonradiative recombination coefficient due to traps (s^-1) Bnrad = 0e-16; % Bimolecular nonradiative recombination coefficient (m^3 s^-1) Caug = 3e-41; % Auger recombination coefficient (m^6s^-1) Dleak = 0e48; % Leakage recombination coefficient (m^13.5s^-1) Nz =10; %number of sections me = m0*0.045; % Effective mass of electron in the CB mhh = m0*0.46; % Effective mass of a heavy hole in the VB mlh = m0*0.056; % Effective mass of a light hole in the VB Nm = 10;

% SOA parameters derived from above parameters

L = (Lc + Lt)*1e-6; % average length

r1 = sqrt(R1); % amplitude reflectivities r2 = sqrt(R2);

eta_in = 10^(-eta_in/10); % converting from dB to linear quantities eta_out = 10^(-eta_out/10);

neq = ((n1^2 - n2^2)*confine/(2+confine) + n2^2)^0.5; delta_Em = h*c/(2*neq*L); % longitudinal mode energy spacing (assumed carrier independent) (J)

% parameters used in nilsson function

global mdh nc nv

mdh = (mhh^1.5 + mlh^1.5)^(2/3); nc = 2*(me*kT/(2*pi*hbar^2))^1.5; nv = 2*(mdh*kT/(2*pi*hbar^2))^1.5;

% parameters used in gain_coeff function

global k0

k1 = c^2*(h^1.5)/(4*sqrt(2)*(pi^1.5)*n1^2); k2 = (2*me*mhh/(hbar*(me + mhh)))^1.5; k0 = k1*k2;

%****************************************************************** disp('SOA versus input power simulation') %constants; % load in SOA and universal constants %simul_params; % load in numerical model parameters and vectors used in spectrum calculations

bias = 100e-3; % bias current (A)

% input signal

Pin_dBm = [-40.0:5:10]; % input power vector in dBm Pin = 1e-3*10.^(Pin_dBm/10); % input signal power in W wavelength_s = 1550*1e-9; % signal wavelength in m Es = h*c./(wavelength_s); % signal energy beta_s = 2*pi*neq*Es/(h*c); % signal propagation coefficient delta_E=h*c./1450*1e-9-Es; % model initial conditions dz=1.22; for Ipin = 1:length(Pin_dBm)

    disp(['Input power = ' num2str(Pin_dBm(Ipin)) ' dBm'])
    Ain = sqrt(Pin(Ipin))./sqrt(Es);  % input signal amplitude (square root of input photon rate)
    weight(1:Nz) = 0.1;  % weighting factor
    n(1:Nz) = 1.2e24; % initial guess for carrier density

    Asp = zeros(Nz+1); % forward travelling signal amplitude
    Asm = zeros(Nz+1); % backward travelling signal amplitude

    Nsp = zeros(Nz+1,Nm+1);  % forward travelling ASE  (single polarisation)
    Nsm = Nsp; % backward travelling ASE (single polarisation)

    oldsignQ(1:Nz) = 1;  % assuming a positive sign for initial value of Q (the algorithm will still converge regardless of reasonable initial conditions)

    tol = 999;  %
   tolerence= 1e-6;
    % Start main iteration 

    while(tol > tolerence) % calculate coefficients of the travelling-wave equations    

        Nsp_old = Nsp;  % storing old values of Nsp and Nsm for use in tolerance calculation
        Nsm_old = Nsm;

        % attenuation coefficient - energy independent    

        alpha_s =  K0 + confine*K1*n/1e24;  % signal attenuation coefficient

        for J = 1:Nm+1;
            alpha(:,J) = K0 + confine*K1*n/1e24;  % attenuation coefficient
        end

        % material gain coefficients and spontaneous emission term in equations (37,38)

        for I = 1:Nz
            dummy = max(n(I),Es);
            gm_s(I) = dummy(1,:);  % material gain coefficient for signal
            dummy = max(n(I),Es);
            gm(I,:) = dummy(1,:);  % material gain coefficient for ASE
            Rsp(I,:) = dummy(1,:);   % Spontaneous emission term in equations (37,38)
        end

        %    boundary conditions - signal

        Asp(1) = (1-r1)*sqrt(eta_in)*Ain + r1*Asm(1);
        Asm(Nz+1) = r2*Asp(Nz+1);

        %    boundary conditions - ASE

        Nsp(1) = R1*Nsm(1);  % input facet
        Nsm(Nz+1) = R2*Nsp(Nz+1);  % output facet

        % solve travelling-wave equations for signal

        for I = 2:(Nz+1)
            Asp(I) = Asp(I-1).*exp((-j*beta_s + 0.5*(confine*gm_s(I-1)-alpha_s(I-1)))*dz);
        end

        for I = Nz:-1:1
            Asm(I) = Asm(I+1).*exp((-j*beta_s + 0.5*(confine*gm_s(I)-alpha_s(I)))*dz);
        end

        % solve travelling-wave equations for ASE

        for I = 2:(Nz+1)
            Nsp(I,:) = Nsp(I-1,:).*exp((confine*gm(I-1,:)-alpha(I-1))*dz) +...
                + Rsp(I-1,:).*(exp((confine*gm(I-1,:)-alpha(I-1))*dz)-1)./(confine*gm(I-1,:)-alpha(I-1));
        end

        for I = Nz:-1:1
            Nsm(I,:) = Nsm(I+1,:).*exp((confine*gm(I,:)-alpha(I))*dz) +...
                Rsp(I,:).*(exp((confine*gm(I,:)-alpha(I))*dz)-1)./(confine*gm(I,:)-alpha(I));
        end

        % normalise noise photon rates

        Gs = exp(sum((confine*gm(I,:)-alpha(I))*dz));    % single pass gain     
        gamma = 4*Gs*sqrt(R1*R2)./(1-sqrt(R1*R2)*Gs).^2;    
        K = 1./sqrt(1+gamma.^2);

        % percentage tolerance using ASE, at energies above maximum bandgap only  

        Eg = max((n));    
        index = max(find(Es <= Eg)) + 1;  

        % tolerence calculation uses ASE throughout the SOA

        tol = max(50*max(abs(Nsp(2:(Nz+1),index:Nm+1) - Nsp_old(2:(Nz+1),index:Nm+1))./Nsp(2:(Nz+1),index:Nm+1) +...
            abs(Nsm(1:Nz,index:Nm+1) - Nsm_old(1:Nz,index:Nm+1))./Nsm(1:Nz,index:Nm+1)));

      disp(['error = ' num2str(tol) ' %'])  % display current error

        % update weights and carrier density 

        for I = 1:Nz        
            Q(I) = bias/(e*d*W*L) - ( (Anrad+Arad)*n(I) + (Bnrad+Brad)*n(I)^2 + Caug*n(I)^3 + Dleak*n(I)^5.5) -...
                0.5*confine/(d*W)*(gm_s(I).*(abs(Asp(I)).^2 + abs(Asp(I+1)).^2 + abs(Asm(I)).^2 + abs(Asm(I+1)).^2)) -...
                confine/(d*W)*sum(gm(I,:).*K.*(Nsp(I,:)+Nsp(I+1,:)+Nsm(I,:)+Nsm(I+1,:)));   % RHS of carrier density rate equation          

            if sign(Q(I)) ~= oldsignQ(I)   % weight adjustment
                weight(I) = weight(I)/2;
            end
            oldsignQ(I) = sign(Q(I));

            if Q(I) > 0    % new carrier density
                n(I) = n(I)*(1+weight(I));
            else
                n(I) = n(I)/(1+weight(I));
            end
        end

    end % of tolerence iteration

    disp('Numerical algorithm converged')

    % Output signal powers

    Aout = (1-r2)*sqrt(eta_out)*Asp(Nz+1);  % output signal amplitude from output facet
    Pout(Ipin) = Es.*(abs(Aout).^2); % output signal power from output facet

    % total ASE output powers

    Nout(Ipin) = 2*eta_out*(1-R2)*sum(K.*Nsp(Nz+1,:).*Es);  % output total ASE power

    % Calculate and displace SOA output power spectrum on an OSA

    sigmaN_spec    = 2*eta_out*(1-R2)*sum(K.*Nsp(Nz+1,:).*Es*h/delta_E);  % output ASE power spectral density (Watts/Hz) for energy vector E

    sigmaN(Ipin) =interp1(Es,sigmaN_spec,Es,'spline'); % interpolating to find output ASE power spectral density at signal wavelength

end

Pout_dBm = db(Pout/1e-3,'power'); Gain = Pout./Pin; % fibre-to-fibre gain NF = sigmaN./(Es.*Gain) + eta_out./Gain; % noise figure (equation 63)

disp('SOA versus input power simulation complete')

% plot signal gain and noise figure vs. bias

subplot(2,1,1) [AX,H1,H2] = plotyy(Pin_dBm,db(Gain,'power'),Pin_dBm,db(NF,'power'),'plot'); xlabel('Input power (dBm)','Fontsize',14); title(['Fibre-to-fibre gain and noise figure at ' num2str(wavelength_s/1e-9) ' nm. Bias current = ' num2str(bias/1e-3) ' mA'],'Fontsize',14)

set(get(AX(1),'Ylabel'),'String','Fibre-to-fibre gain (dB)','Fontsize',14) set(AX(1),'XLim',[min(Pin_dBm) max(Pin_dBm)]); set(AX(1),'YLim',[-10 30],'Ytick',[-10:5:30]); set(H1,'LineStyle','-','Marker','+')

set(get(AX(2),'Ylabel'),'String','Noise figure (dB)','Fontsize',14) set(H2,'LineStyle','-','Marker','o'); set(H2,'LineStyle','-.'); set(AX(2),'XLim',[min(Pin_dBm) max(Pin_dBm)]); set(AX(2),'YLim',[0 20],'Ytick',[0:5:20]);

subplot(2,1,2) [AX,H1,H2] = plotyy(Pout_dBm,db(Gain,'power'),Pout_dBm,db(NF,'power'),'plot'); xlabel('Output power (dBm)','Fontsize',14); title(['Fibre-to-fibre gain and noise figure at ' num2str(wavelength_s/1e-9) ' nm. Bias current = ' num2str(bias/1e-3) ' mA'],'Fontsize',14)

set(get(AX(1),'Ylabel'),'String','Fibre-to-fibre gain (dB)','Fontsize',14) set(AX(1),'XLim',[-10 10]); set(AX(1),'YLim',[-10 30],'Ytick',[-10:5:30]); set(H1,'LineStyle','-','Marker','+')

set(get(AX(2),'Ylabel'),'String','Noise figure (dB)','Fontsize',14) set(H2,'LineStyle','-','Marker','o'); set(H2,'LineStyle','-.'); set(AX(2),'XLim',[-10 10]); set(AX(2),'YLim',[0 20],'Ytick',[0:5:20]);

[Adam]

Please read this on formatting code:

https://uk.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup#answer_18099

You haven't really asked a question, you just threw a load of code on the screen. What do you expect people to do with it?

At a glance it has something of the order of 30 global variables thrown into it so no wonder you are having trouble with the code, it could basically be doing anything and I doubt anyone here is going to be able to help you understand such code - it is incomprehensible, even if it were formatted.

[Jan]

Sorry, Edwin, I'm not even able to read your code. Please apply a proper formatting. Thanks.

If you want somebody to help you, it would be great if you mention the occurring problems. Do you get an error message?

Answer 1

[Yazan Khlefat]

The Correct Code is in attachment

Good Luck

[Guillaume]

One would hope that Edwin is no longer waiting for help after more than a year!

Nonetheless, it would be useful if you explained what you modified as it's not obvious.