clear all
close all
clc

global P_pump f gamma_B L

L=10e-3;                 % effective interaction length, unit: km
gamma_B=50e6;            % Brillouin linewidth, unit: Hz
P_pump=100;              % pump power, unit: mW, modulation instability threshold 20 dBm
f_offset=250e6;                   % scanning frequency range, unit: Hz
f_step=1e3;                       % frequency step, unit: Hz
f=-f_offset:f_step:f_offset;      % frequenct detuning span
alpha_log=0.2;                    % fiber logarithm attenuation, unit: dB/km
alpha_lin=alpha_log/4.343;        % fiber linear attenuation, unit: /km

L_fiber=60;                       % fiber length, unit: km
delta_L=50e-3;                    % fiber section length, unit: km
l=0:delta_L:L_fiber;

m=0.1:0.01:1;                     % ratio of maximum gain and maximum loss
d=0.4:0.02:3;                     % frequency seperation ratio to the FWHM of gain profile

n=500;                            % times of iteration

noise_level=0.012;              % noise level

for j=1:length(l)
    
    sigma_g_g=0;                              % initial value of standard deviation of frequency error for gain only case
    sigma_g_gl=zeros(size(d,2),size(m,2));    % initial value of standard deviation of frequency error for gain with losses case
    
    for s=1:n                                 % iteration
                
        [SBS_g_log SBS_g_phase]=SBS_g(P_pump*exp(-alpha_lin*l(j)));                 % generation of SBS gain profile    
        SBS_g_log_n=SBS_g_log+noise_level*randn(1,size(SBS_g_log,2));             % additional noise on gain spectrum of gain-only case
        
        % least square method fit for gain-only case
        % fit function: Lorentzian

        % initial values of fitting          
        P=P_pump*exp(-alpha_lin*l(j));      % P_pump
        gamma=gamma_B;                      % width             
        f0=0;                               % peak center
        y0=0;                               % baseline
           
        vStart_g_g=[P,gamma,f0,y0];
        warning('off','all')
        vEnd_g_g=nlinfit(f,SBS_g_log_n,@Lorentz_g_gain_fun,vStart_g_g);
    
        delta_g_g=vEnd_g_g(3);              % peak frequency error   
        
        % standard deviation calculation
        sigma_g_g=sigma_g_g+delta_g_g^2;    % standard deviation of frequency error for gain spectrum in gain only case
       
        for i=1:length(m)
            for k=1:length(d)

                [SBS_gl_log SBS_gl_phase]=SBS_gl(P_pump*exp(-alpha_lin*l(j)),m(i),d(k));     % generation of SBS gain with losses prifle
                SBS_gl_log_n=SBS_gl_log+noise_level*randn(1,size(SBS_gl_log,2));           % additional noise on gain spectrum of gain with losses case              

                % least square method fit for gain with losses case
                % fit function: superimposed Lorentzian

                % additional initial values of fitting
                mm=m(i);                        % m value
                dd=d(k);                        % d value

                vStart_g_gl=[P,gamma,f0,y0,mm,dd];
                vEnd_g_gl=nlinfit(f,SBS_gl_log_n,@Lorentz_gl_gain_fun,vStart_g_gl);

                delta_g_gl=vEnd_g_gl(3);          % frequency offset

                % standard deviation calculation           
                sigma_g_gl(k,i)=sigma_g_gl(k,i)+delta_g_gl^2;        % standard deviation of frequency error for gain with losses case            
            end
        end
        
        f_g_g(j)=sqrt(sigma_g_g/n);                % standard deviation of frequency error for gain only case, unit: Hz
        
        for i=1:length(m)
            for k=1:length(d)
                f_g_gl(j,k,i)=sqrt(sigma_g_gl(k,i)/n);     % standard deviation of frequency error gain with losses case, unit: Hz
            end
        end
    end
end
