j_f_peng.m

孤子传输孤子传输孤子传输孤子传输孤子传输孤子传输

clear
clc
% h 迭代时选择的步长
% dT 选择的时间步长
% T 时间坐标
% z 位置坐标
% n 将T分成N个区域dT=T/N,n=0,1,2,3,...,N-1
% w0 w0=2*pi/N
% w 文章中采用负频故w=-w0*n,n,w均是个向量
% gn 关于的函数
% A 单根光纤中模场的慢变振幅
% At1 双芯耦合器第一个纤芯中模场的归一化慢变振幅
% At2 双芯耦合器第二个纤芯中模场的归一化慢变振幅
% Ath1 三芯耦合器第一个纤芯中模场的归一化慢变振幅
% Ath2 三芯耦合器第二个纤芯中模场的归一化慢变振幅
% Ath3 三芯耦合器第三个纤芯中模场的归一化慢变振幅
h=0.0005;
dT=0.01;
T=-15:dT:15;                    
z=0:h:10;                        
n=0:1:(length(T)-2);                       
w0=2*pi/(length(T)-1);           
w=-w0*n;                         
gn=-i*w.^2/2;                    
A(1,:)=sech(T(1):dT:T(length(T)-1));        
At1(1,:)=sech(T(1):dT:T(length(T)-1));      
At2(1,:)=zeros(length(T)-1,1);           
Ath1(1,:)=sech(T(1):dT:T(length(T)-1));     
Ath2(1,:)=zeros(length(T)-1,1);          
Ath3(1,:)=zeros(length(T)-1,1);          
% 对称分步傅里叶变换法求解光脉冲在N芯耦合器中的耦合模方程组的数值解的迭代过程
% 由于这里只考虑二阶色散，故可以不考虑fft()和ifft()的交换问题
% 用fftshift(),ifftshift()将零频矢量移到频率中心
for k=1:1:length(z)-1
% 孤子（光脉冲）在单根光纤中的传输演化过程
    An=fftshift(fft(exp(gn*h/2).*fftshift(ifft(A(k,:)))));
    Ap=An.*exp(i*h*abs(A(k,:)).^2);
    A(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ap))));
    for l=1:2
        Ap=An.*exp(i*h/2*(abs(A(k,:)).^2+abs(A(k+1,:)).^2));
        A(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ap))));
    end
% 孤子（光脉冲）在双芯耦合器中的传输演化过程
    At1n=fftshift(fft(exp(gn*h/2).*fftshift(ifft(At1(k,:))+i*h/2*ifft(At2(k,:)))));
    At1p=At1n.*exp(i*h*abs(At1(k,:)).^2);
    At2n=fftshift(fft(exp(gn*h/2).*fftshift(ifft(At2(k,:))+i*h/2*ifft(At1(k,:)))));
    At2p=At2n.*exp(i*h*abs(At2(k,:)).^2);
    At1(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(At1p)+i*h/2*ifft(At2p))));
    At2(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(At2p)+i*h/2*ifft(At1p))));
    for l=1:2
        At1p=At1n.*exp(i*h/2*(abs(At1(k,:)).^2+abs(At1(k+1,:)).^2));
        At2p=At2n.*exp(i*h/2*(abs(At2(k,:)).^2+abs(At2(k+1,:)).^2));
        At1(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(At1p)+i*h/2*ifft(At2p))));
        At2(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(At2p)+i*h/2*ifft(At1p))));
    end
% 孤子（光脉冲）在三芯耦合器中的传输演化过程
    Ath1n=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath1(k,:))+i*h/2*(ifft(Ath2(k,:))))));
    Ath1p=Ath1n.*exp(i*h*abs(Ath1(k,:)).^2);
    Ath2n=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath2(k,:))+i*h/2*(ifft(Ath1(k,:))+ifft(Ath3(k,:))))));
    Ath2p=Ath2n.*exp(i*h*abs(Ath2(k,:)).^2);  
    Ath3n=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath3(k,:))+i*h/2*(ifft(Ath2(k,:))))));
    Ath3p=Ath3n.*exp(i*h*abs(Ath3(k,:)).^2);    
    Ath1(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath1p)+i*h/2*ifft(Ath2p))));
    Ath2(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath2p)+i*h/2*(ifft(Ath1p)+ifft(Ath3p)))));
    Ath3(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath3p)+i*h/2*ifft(Ath2p)))); 
    for l=1:2
        Ath1p=Ath1n.*exp(i*h/2*(abs(Ath1(k,:)).^2+abs(Ath1(k+1,:)).^2));
        Ath2p=Ath2n.*exp(i*h/2*(abs(Ath2(k,:)).^2+abs(Ath2(k+1,:)).^2));
        Ath3p=Ath3n.*exp(i*h/2*(abs(Ath3(k,:)).^2+abs(Ath3(k+1,:)).^2));
        Ath1(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath1p)+i*h/2*ifft(Ath2p))));
        Ath2(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath2p)+i*h/2*(ifft(Ath1p)+ifft(Ath3p)))));
        Ath3(k+1,:)=fftshift(fft(exp(gn*h/2).*fftshift(ifft(Ath3p)+i*h/2*ifft(Ath2p)))); 
    end
end
plot(T(1:1:length(T)-1),A(1,:))
xlabel('t')
ylabel('A')
title('光脉冲A=sech(t)')
figure
waterfall(T(1:1:length(T)-1),z(1:(length(z)-1)/50:length(z)-1),abs(A(1:(length(z)-1)/50:length(z)-1,:)))
%mesh(T(1:1:length(T)-1),z(1:1:length(z)-1),abs(A(1:1:length(z)-1,:)))
xlabel('t')
ylabel('z')
zlabel('A')
axis([T(1) T(length(T)) 0 z(length(z)) 0 1])
title('光脉冲在单根光纤中的传输')
view(129,44)
figure
waterfall(T(1:1:length(T)-1),z(1:(length(z)-1)/50:length(z)-1),abs(At1(1:(length(z)-1)/50:length(z)-1,:)))
%mesh(T(1:1:length(T)-1),z(1:1:length(z)-1),abs(At1(1:1:length(z)-1,:)))
xlabel('t')
ylabel('z')
zlabel('At1')
axis([T(1) T(length(T)) 0 z(length(z)) 0 1])
title('双芯耦合器中光脉冲在第一根光纤中传输')
view(129,44)
figure
waterfall(T(1:1:length(T)-1),z(1:(length(z)-1)/50:length(z)-1),abs(At2(1:(length(z)-1)/50:length(z)-1,:)))
%mesh(T(1:1:length(T)-1),z(1:1:length(z)-1),abs(At2(1:1:length(z)-1,:)))
xlabel('t')
ylabel('z')
zlabel('At2')
axis([T(1) T(length(T)) 0 z(length(z)) 0 1])
title('双芯耦合器中光脉冲在第二根光纤中传输')
view(129,44)
figure
waterfall(T(1:1:length(T)-1),z(1:(length(z)-1)/50:length(z)-1),abs(Ath1(1:(length(z)-1)/50:length(z)-1,:)))
%mesh(T(1:1:length(T)-1),z(1:1:length(z)-1),abs(Ath1(1:1:length(z)-1,:)))
xlabel('t')
ylabel('z')
zlabel('Ath1')
axis([T(1) T(length(T)) 0 z(length(z)) 0 1])
title('三芯耦合器中光脉冲在第一根光纤中传输')
view(129,44)
figure
waterfall(T(1:1:length(T)-1),z(1:(length(z)-1)/50:length(z)-1),abs(Ath2(1:(length(z)-1)/50:length(z)-1,:)))
%mesh(T(1:1:length(T)-1),z(1:1:length(z)-1),abs(Ath2(1:1:length(z)-1,:)))
xlabel('t')
ylabel('z')
zlabel('Ath2')
axis([T(1) T(length(T)) 0 z(length(z)) 0 1])
title('三芯耦合器中光脉冲在第二根光纤中传输')
view(129,44)
figure
waterfall(T(1:1:length(T)-1),z(1:(length(z)-1)/50:length(z)-1),abs(Ath3(1:(length(z)-1)/50:length(z)-1,:)))
%mesh(T(1:1:length(T)-1),z(1:1:length(z)-1),abs(Ath3(1:1:length(z)-1,:)))
xlabel('t')
ylabel('z')
zlabel('Ath3')
axis([T(1) T(length(T)) 0 z(length(z)) 0 1])
title('三芯耦合器中光脉冲在第三根光纤中传输')
view(129,44)