j_f_peng.m

此源码是文档simulation.rar__关于光孤子在光纤耦合器中传输的仿真研究的matlab程序

% USE: Simulation of Optical Soliton Transmission in Optical Fiber Coupler
% AUTHOR: J.F.Peng ([e-mail])
% Date：******2008.03.31******
% Copyright: wspace.m(Thomas E. Murphy),J_F_Peng.m(J.F.Peng)
% 理论指导: 周青春. 江苏科技大学. 数理学院. 应用物理专业
% 理论来源：施娟,侯韶华. 光纤耦合器中光孤子传输仿真研究[J]--详见参考文献
clear
clc
% h 迭代时选择的步长
% z 位置坐标
% T 时间窗口(-10:10)
% nt 时间区域个数（2^M个）
% dt 选择的时间步长
% t 时间坐标
% w 傅里叶域频率
% gn 关于w的函数
% A 单根光纤中模场的慢变振幅
% At1 双芯耦合器第一个纤芯中模场的归一化慢变振幅
% At2 双芯耦合器第二个纤芯中模场的归一化慢变振幅
% Ath1 三芯耦合器第一个纤芯中模场的归一化慢变振幅
% Ath2 三芯耦合器第二个纤芯中模场的归一化慢变振幅
% Ath3 三芯耦合器第三个纤芯中模场的归一化慢变振幅 
h=0.001;                    
z=0:h:10;     
T = 20;                                 
nt = 2^10;                              
dt = T/nt;                             
t = ((1:nt)'-(nt+1)/2)*dt;              
w = wspace(T,nt); 
% 上式引用了Thomas E. Murphy的程序wspace.m计算圆频率
gn=-i*w'.^2/2;
% 输入初始脉冲
A(1,:)=sech(t);        
At1(1,:)=sech(t);      
At2(1,:)=zeros(length(t),1);           
Ath1(1,:)=sech(t);     
Ath2(1,:)=zeros(length(t),1);          
Ath3(1,:)=zeros(length(t),1);          
% 对称分步傅里叶变换法求解光脉冲在N芯耦合器中的耦合模方程组的数值解的迭代过程
% 由于这里只考虑二阶色散，故可以不考虑fft()和ifft()的交换问题（本程序有考虑到）
for k=1:1:length(z)-1
% 孤子（光脉冲）在单根光纤中的传输演化过程
    An=fft(exp(gn*h/2).*ifft(A(k,:)));
    Ap=An.*exp(i*h*abs(A(k,:)).^2);
    A(k+1,:)=fft(exp(gn*h/2).*ifft(Ap));
    for l=1:2
        Ap=An.*exp(i*h/2*(abs(A(k,:)).^2+abs(A(k+1,:)).^2));
        A(k+1,:)=fft(exp(gn*h/2).*ifft(Ap));
    end
% 孤子（光脉冲）在双芯耦合器中的传输演化过程
    At1n=fft(exp(gn*h/2).*(ifft(At1(k,:))+i*h/2*ifft(At2(k,:))));
    At1p=At1n.*exp(i*h*abs(At1(k,:)).^2);
    At2n=fft(exp(gn*h/2).*(ifft(At2(k,:))+i*h/2*ifft(At1(k,:))));
    At2p=At2n.*exp(i*h*abs(At2(k,:)).^2);
    At1(k+1,:)=fft(exp(gn*h/2).*(ifft(At1p)+i*h/2*ifft(At2p)));
    At2(k+1,:)=fft(exp(gn*h/2).*(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,:)=fft(exp(gn*h/2).*(ifft(At1p)+i*h/2*ifft(At2p)));
        At2(k+1,:)=fft(exp(gn*h/2).*(ifft(At2p)+i*h/2*ifft(At1p)));
    end
% 孤子（光脉冲）在三芯耦合器中的传输演化过程
    Ath1n=fft(exp(gn*h/2).*(ifft(Ath1(k,:))+i*h/2*(ifft(Ath2(k,:)))));
    Ath1p=Ath1n.*exp(i*h*abs(Ath1(k,:)).^2);
    Ath2n=fft(exp(gn*h/2).*(ifft(Ath2(k,:))+i*h/2*(ifft(Ath1(k,:))+ifft(Ath3(k,:)))));
    Ath2p=Ath2n.*exp(i*h*abs(Ath2(k,:)).^2);  
    Ath3n=fft(exp(gn*h/2).*(ifft(Ath3(k,:))+i*h/2*(ifft(Ath2(k,:)))));
    Ath3p=Ath3n.*exp(i*h*abs(Ath3(k,:)).^2);    
    Ath1(k+1,:)=fft(exp(gn*h/2).*(ifft(Ath1p)+i*h/2*ifft(Ath2p)));
    Ath2(k+1,:)=fft(exp(gn*h/2).*(ifft(Ath2p)+i*h/2*(ifft(Ath1p)+ifft(Ath3p))));
    Ath3(k+1,:)=fft(exp(gn*h/2).*(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,:)=fft(exp(gn*h/2).*(ifft(Ath1p)+i*h/2*ifft(Ath2p)));
        Ath2(k+1,:)=fft(exp(gn*h/2).*(ifft(Ath2p)+i*h/2*(ifft(Ath1p)+ifft(Ath3p))));
        Ath3(k+1,:)=fft(exp(gn*h/2).*(ifft(Ath3p)+i*h/2*ifft(Ath2p))); 
    end
end
save('J_F_Peng','t','z','A','At1','At2','Ath1','Ath2','Ath3')
% waterfall 用于画流水线图样，mesh 用于画彩色曲面，注意采样点少点好画
figure
plot(t(1:1:length(t)),A(1,:))
xlabel('t')
ylabel('A')
title('光脉冲A=sech(t)')
figure
waterfall(t(1:1:length(t)),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)),z(1:(length(z)-1)/200:length(z)-1)',abs(A(1:(length(z)-1)/200: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)),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)),z(1:(length(z)-1)/200:length(z)-1)',abs(At1(1:(length(z)-1)/200: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)),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)),z(1:(length(z)-1)/200:length(z)-1)',abs(At2(1:(length(z)-1)/200: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)),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)),z(1:(length(z)-1)/200:length(z)-1)',abs(Ath1(1:(length(z)-1)/200: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)),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)),z(1:(length(z)-1)/200:length(z)-1)',abs(Ath2(1:(length(z)-1)/200: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)),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)),z(1:(length(z)-1)/200:length(z)-1)',abs(Ath3(1:(length(z)-1)/200: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)
%参考文献:
%[1] 施娟,侯韶华. 光纤耦合器中光孤子传输的仿真研究[J]. 电子元器件应用. 2008,Vol.10,No.4:65-67.
%[2] 张晓光等. 非线性薛定谔方程数值解法中傅立叶正逆变换选取的讨论[J]. 计算物理. 2003,Vol.20,No.3:267-272. 
%[3] 王志斌,李志全. 光孤子传输演化的分步傅里叶法研究[J]. 应用光学. 2007,Vol.28,No.1:82-85.
%[4] 项鹏,郑亚彬. 光纤波导中光脉冲演化方程的数值解及计算机仿真研究[J]. 光子技术. 2006(2):117-120.
%[5] 施娟. 基于对称分步傅立叶算法的光孤子仿真[J]. 电子元器件应用. 2008,Vol.10,No.1:73-75.
%[6] 吕理想,张晓萍. 不同形式非线性薛定谔方程及其分步傅里叶法求解[J]. 计算物理. 2007,Vol.24,No.3:373-377.
%[7] C.查布拉,P.卡纳尔. 工程数值方法(第5版)[M]. 北京: 清华大学出版社,2007.12.
%[8] P.阿戈沃. 非线性光纤光学原理及应用[M]. 北京: 电子工业出版社,2002.12.