非线性光纤光学

%观察群速度色散和kell效应同时扰动作用下孤子对的传输特性


clear
clc

%?注意时间域与频域的选择问题
T = 40;                               
tn= 2^10;  % 时间区域个数（2^M个）                            
dt = T/tn;   %选择的时间步长  
t = ((1:tn)'-(tn+1)/2)*dt;  
w = wn(T,tn);%傅里叶域频率,调用计算频率函数 

m=ran(0,1)%n为0-1之间随机变化数
r=ran(-1,1)%l为-1和1之间随机变化数
%上面两个式子调用了随机数产生函数

b=0.03;
c=-i*(1+m*r)*w'.^2/2+i*b*w'.^3;
%步长的选择影响计算效率
dz=0.1;  %迭代时选择的步长
z=0:dz:100;  % 位置坐标  
%z=0:dz:50;  % 位置坐标  

n=ran(0,1)%n为0-1之间随机变化数
l=ran(-1,1)%l为-1和1之间随机变化数
%上面两个式子调用了随机数产生函数

% 输入初始脉冲
q=3.5;r=1;%q为初始间距，r为相对振幅
u(1,:)=sech(t-q)+r*sech(r*(t+q));  

%?注意复数的表示，傅里叶变换函数的使用
%由于输入负频，fft对应逆变换，ifft对应傅里叶变换
for k=1:1:length(z)
% 分布傅里叶变换解方程
    u1=fft(exp(c*dz/2).*ifft(u(k,:)));
    u2=u1.*exp(i*dz*(1+n*l)*abs(u(k,:)).^2);
    u(k+1,:)=fft(exp(c*dz/2).*ifft(u2));
    %二次迭代
    for p=1:2
        u2=u1.*exp(i*dz/2*(1+n*l)*(abs(u(k,:)).^2+abs(u(k+1,:)).^2));
        u(k+1,:)=fft(exp(c*dz/2).*ifft(u2));
    end
%缺少迭代仿真将出现很多问题
end 

figure

% waterfall 用于画流水线图样，mesh 用于画彩色曲面，mesh计算速度快
%waterfall(t(1:1:length(t)),z(1:(length(z))/50:length(z)),abs(u(1:(length(z))/50:length(z),:)))
mesh(t(1:1:length(t)),z(1:(length(z))/50:length(z)),abs(u(1:(length(z))/50:length(z),:)))
%注意采样点少点好画，
xlabel('归一化时间')
ylabel('归一化距离','rotation',-65)
zlabel('归一化振幅')
axis([-20,20,0,100,0,1.5]);
title('GVD和KELL同时扰动作用输出光脉冲')