clear
clc
tic;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%Parameters%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pp = 4.2e-3 ; %W
Ps = 48.6e-6; %W
N = 1.4677; % refractive index
NA=0.2; % numerical aperture
r = 9.2E-6/2; % core radius (m)
c0 = 2.9979E8; % light velocity
lambda = 1550E-9; % laser wavelength (both pump & probe)
va = 5960; % acoustic velocity
rho = 2.201E3; % density
p12 = 0.2678; % longitudinal elacto-optical coefficient
gamma = 0.5; % polarisation factor
Delta_fB = 36E6; % Brillouin gain linewidth
Att_km_fiber = 0.217; % fiber attenuation (dB/km)
for Ps = 48.6e-6 ; %W
%%%%%%%%%%%%%%%%
%Effective area%
%%%%%%%%%%%%%%%%
A_eff = 86E-12; % Effective area value to validate the model
L_fiber_begin=1e3; % length of the fiber
L_fiber_step=10e3; % step for the loop 1e3
L_fiber_end=100e3; % length of the fiber
L_fiber = L_fiber_begin:L_fiber_step:L_fiber_end;
%L_fiber = L_fiber_begin:L_fiber_end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Brillouin frequency shift%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
f_B=2*N*va/lambda; % Brillouin frequency shift%
L_i = size(L_fiber);
a=0; % used for the loop
s=0; % used for the loop
alpha_m = Att_km_fiber/(10*log10(exp(1)))/1000; % fiber attenuation in 1/m
gBmax =1.2E-11; % Brillouin frequency shift to validate the model
Ip = Pp/A_eff; % Pump intensity in W/m^A
Is_in = Ps/A_eff; % Probe intensity in W/m^A
for L_fiber = L_fiber_begin:L_fiber_step:L_fiber_end; % in meter
a=a+1;
%%%%%%%%%%%%%%%%%%%%%%%%%
%Attenuation coefficients%
Att_fiber = (L_fiber./1000).*Att_km_fiber;
%%%%%%%%%%%%%%%%%
%Effective length
L_eff_fiber = (1 - exp((-1).*alpha_m.*(L_fiber)))./alpha_m;
%%%%%%%%%%%%%%%
%Brillouin gain
%gBmax = (2*pi*N^7.*p12.^2*gamma)./(c0*lambda^2*rho*va*Delta_fB) % m/W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Resolution of the non-linear coupled waves system%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s=0;
while 1
s=s+ 1;
A=Is_in/(s)*exp(gBmax*Pp*L_eff_fiber/A_eff -alpha_m.*L_fiber);
Is=A;
sys = @(z,I) [
-gBmax*I(1)*I(2) - alpha_m*I(1);
-gBmax*I(1)*I(2) + alpha_m*I(2);
];
[z,I] = ode45(sys,[0,L_fiber],[Ip Is]);
clear A
if (abs((Ps - I(end,2)*A_eff))/Ps) > 0.6;
clear I
clear z
else
break
end
end
P_Stokes(a)=I(1,2)*A_eff*1e3;
end
ind= find(P_Stokes==max(P_Stokes));
L_opt = (L_fiber_begin)+(ind-1).*L_fiber_step;
L_i;
a=1:length(P_Stokes);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%Plots%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
hold on
plot(((L_fiber_begin)+(a-1).*L_fiber_step)/1000,P_Stokes,'-')
hold on
figure(2)
 plot(z/1000,I(:,2)*A_eff*1e3)
 hold on
 plot(z/1000,I(:,1)*A_eff*1e3,'r-.')
 %plot((L_fiber_begin)+(a-1),P_Stokes)
%axis([0:10:100 1:0.2:2])
%legend('Stokes wave at the fiber input')
xlabel('Longueur de fibre en km')
ylabel('Puissance en mW')
end
%figure(2)
%plot(z/1000,I(:,2)*A_eff*1e3)
%hold on
toc