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)