Почему я получаю противоречивые графики (не линейно уменьшающиеся по мере необходимости) с изменением количества рассеивателей? - PullRequest
Новая
       56

Почему я получаю противоречивые графики (не линейно уменьшающиеся по мере необходимости) с изменением количества рассеивателей?

0 голосов
/ 28 апреля 2020

Я пытаюсь смоделировать мощность и затухание обратного рассеянного сигнала для фазово-оптической рефлектометрии во временной области (phi-OTDR). Я использовал формулу, приведенную в научных работах:

formula

и общую мощность p = p1 + p2

Я должен получить график зависимости длины от общей мощности в линейно убывающим образом для любого значения N (количество рассеивателей), как показано ниже. Здесь N = 100.

graph for N=100

Но когда я изменяю N на любое другое значение, я не получаю тот же результат. Вместо этого я получаю что-то вроде этого: Здесь, N = 200:

graph for N=200

Мой код, как показано ниже: 


NA = 0.13; %numerical aperture
ncore = 1.47 ; %fiber core refractive index
ndiff = 3e-3 ; %difference in index of refraction
lamda = 1.55 ;%Lambda in micrometer
c = 3e5; %speed of light in km/s
f = c/lamda*1e9; %pulse frequency
m = 4.55 ; %for single mode fiber
S = (NA/ncore)^2 / m ; %Capture coefficient 
asdB= (0.76+0.51*ndiff)/lamda^4; % Rayleigh attenuationin dB/km
as=asdB/4.34; % Rayleigh attenuation in km-1
adB= 0.22 ; % attenuationin dB/km
a=adB/4.34; % attenuation in km-1
ng=1.4681; % group refractive index
vg=c/ng; %group velocity in km/s
W = 1e-7; %Pulse width 100 ns
Z = 1; %distance of the fiber in km
N = 100; %number of scatterers
z = 0:Z/N:Z; %scatterers across the fiber length
fa=1.1e8; %Acousto Optic Modulator frequency in Hz
tau = (2*ng*z)/c;%delay time of scattering centers
t1=2*ng*Z/c; %roundtrip time
t=0:t1/100:t1; % at 100m resolution
L=vg*t/2; %length of fiber with 100m resolution- so steps of 0.01
Pin=0.8; % optical power of input lightwave

%% Sum of the optical power of the backscattered lightwave from each independent scatterer without distrbtion
P1=0;
for i=1:N+1
    P1 = P1 + Pin*S*as*exp((-2*a*c*tau(i))/ng)*(heaviside(t-tau(i)+W/2)-heaviside(t-tau(i)-W/2)); 
end
P2=0;
for i=1:N+1
    for j=i+1:N+1
        phi=2*pi*(f+fa)*(tau(i)-tau(j));
        P2 = P2 + Pin*cos(phi)*S*as*exp(-a*c*(tau(i)+tau(j))/ng)*(heaviside(t-tau(i)+W/2)-heaviside(t-tau(i)-W/2))...
            .*(heaviside(t-tau(j)+W/2)-heaviside(t-tau(j)-W/2));
    end
end
P2 = 2*P2;
Ptotal = P1+P2;
figure(1);
plot(L,Ptotal);
grid on;
xlabel('Distance in km ');
ylabel('Output Power in Watts ');
...