% USE: Simulation of Optical Soliton Transmission in Optical Fiber Coupler
% AUTHOR: J.F.Peng (jingyujiafu@163.com)
% Date：******2008.03.31******
% Copyright: wspace.m(Thomas E. Murphy),J_F_Peng.m(J.F.Peng)
% 理論指導(dǎo): 周青春. 江蘇科技大學(xué). 數(shù)理學(xué)院. 應(yīng)用物理專業(yè)
% 理論來源：施娟,侯韶華. 光纖耦合器中光孤子傳輸仿真研究[J]--詳見參考文獻(xiàn)
clear
clc
% h 迭代時選擇的步長
% z 位置坐標(biāo)
% T 時間窗口(-10:10)
% nt 時間區(qū)域個數(shù)（2^M個）
% dt 選擇的時間步長
% t 時間坐標(biāo)
% w 傅里葉域頻率
% gn 關(guān)于w的函數(shù)
% 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芯耦合器中的耦合模方程組的數(shù)值解的迭代過程
% 由于這里只考慮二階色散，故可以不考慮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)
%參考文獻(xiàn):
%[1] 施娟,侯韶華. 光纖耦合器中光孤子傳輸?shù)姆抡嫜芯縖J]. 電子元器件應(yīng)用. 2008,Vol.10,No.4:65-67.
%[2] 張曉光等. 非線性薛定諤方程數(shù)值解法中傅立葉正逆變換選取的討論[J]. 計算物理. 2003,Vol.20,No.3:267-272. 
%[3] 王志斌,李志全. 光孤子傳輸演化的分步傅里葉法研究[J]. 應(yīng)用光學(xué). 2007,Vol.28,No.1:82-85.
%[4] 項鵬,鄭亞彬. 光纖波導(dǎo)中光脈沖演化方程的數(shù)值解及計算機(jī)仿真研究[J]. 光子技術(shù). 2006(2):117-120.
%[5] 施娟. 基于對稱分步傅立葉算法的光孤子仿真[J]. 電子元器件應(yīng)用. 2008,Vol.10,No.1:73-75.
%[6] 呂理想,張曉萍. 不同形式非線性薛定諤方程及其分步傅里葉法求解[J]. 計算物理. 2007,Vol.24,No.3:373-377.
%[7] C.查布拉,P.卡納爾. 工程數(shù)值方法(第5版)[M]. 北京: 清華大學(xué)出版社,2007.12.
%[8] P.阿戈沃. 非線性光纖光學(xué)原理及應(yīng)用[M]. 北京: 電子工業(yè)出版社,2002.12.