分步傅里叶算法求解GNLSE

%光脉冲参数定义Ppeak=20;%输入脉冲峰值功率，单位：WFWHM=300*10^-15;%输入脉冲半高全宽，单位：sT0=FWHM/(2*log(1+sqrt(2)));%输入双曲正割脉冲1/e脉宽的一半，单位：sC=0;%输入脉冲啁啾参量lamda=1550*10^-9;%输入脉冲中心波长，单位：mc=3*10^8;%真空中的光速，单位：m/sw0=2*pi*c/lamda;%输入脉冲中心角频率，单位：rad/s%时频窗口定义num=2^16;%时频域内的离散点数Ttotal=500*(2*T0);%时间窗口宽度，单位：sdeltaT=Ttotal/(num-1);%时间分辨率，单位：sT=(-num/2:num/2-1)*deltaT;%时间向量，单位：sdeltaw=2*pi/Ttotal;%角频率分辨率，单位：rad/sw=deltaw*[(0:num/2-1)(-num/2:-1)];%角频率向量，单位：rad/s%输入脉冲时频分析Et=sqrt(Ppeak)*sech(T/T0).*exp(-j*C/2*(T/T0).^2);%双曲正割脉冲的时域电场分布Ew=ifft(Et);%双曲正割脉冲的频域电场分布%输入脉冲时频特性图形显示figure(1)subplot(2,1,1)plot(T*10^12,(abs(Et)).^2,'k','linewidth',2)title('输入光脉冲功率在时域和频域内的分布')xlabel('时间（ps）')ylabel('功率（W）')axis([-10100inf])subplot(2,1,2)plot(2*pi*c./(fftshift(w)+w0)*10^9,(abs(fftshift(Ew))).^2/max((abs(fftshift(Ew))).^2),'k','linewidth',2)xlabel('波长（nm）')ylabel('归一化光功率谱密度')axis([152015800inf])%光纤参数的定义beta2=-3.8237e-27;%GVD参数，单位s2/mbeta3=6.3147e-41;%TOD参数，单位s3/mgama=10e-3;%非线性系数，单位1/(W.m)Tr=3e-15;%拉曼响应函数关于时间的一阶矩，单位：salfa=1.5e-3/4.343;%光纤损耗系数，单位：1/mL=1000;%光纤长度，单位：mh=1;%空间离散步长，单位：mN=L/h+1;%空间离散点数%计算中所用到的矩阵定义Etold=zeros(1,num);%用于存储电场随时间的分布Ewold=zeros(1,num);%用于存储电场随角频率的分布Raman=zeros(1,num);%用于存储时域内电场模的平方对时间的导数Self=zeros(1,num);%用于存储时域内电场模的平方与电场乘积对时间的导数%光纤输入端（z=0）的电场时频分布Etold=Et;%z=0处时域电场分布Ewold=ifft(Etold);%z=0处角频率域电场分布%前h/2色散影响的计算dispersion_half=exp(j*(beta2*(w.^2)*(h/2)/2+beta3*(w.^3)*(h/2)/6)-alfa*(h/2)/2);dispersion=exp(j*(beta2*(w.^2)*h/2+beta3*(w.^3)*h/6)-alfa*h/2);Ewold=Ewold.*dispersion_half;%计算第一个半步长色散的影响Etold=fft(Ewold);%N-2段周期性h非线性和h色散影响的计算fornumber=1:N-2%求h长度内SPM、SST和Raman作用Etsqure=(abs(Etold)).^2;Etold_conj=conj(Etold);Raman=gradient(Etsqure,deltaT);Self=Raman+Etold_conj.*gradient(Etold,deltaT);Etold=Etold.*exp(h*j*gama*(Etsqure-Tr*Raman+j/w0*Self));Ewold=ifft(Etold);%求h长度光纤损耗、GVD和TOD作用Ewold=Ewold.*dispersion;%逆傅里叶变换求时域电场分布Etold=fft(Ewold);end%最后一个h非线性和h/2色散影响的计算%求最后一个h长度中非线性影响Etsqure=(abs(Etold)).^2;Etold_conj=conj(Etold);Raman=gradient(Etsqure,deltaT);Self=Raman+Etold_conj.*gradient(Etold,deltaT);Etold=Etold.*exp(h*j*gama*(Etsqure-Tr*Raman+j/w0*Self));Ewold=ifft(Etold);%求最后一个h/2长度中色散影响Ewold=Ewold.*dispersion_half;Etold=fft(Ewold);%输出脉冲时频特性图形显示figure(2)subplot(2,1,1)plot(T*10^12,(abs(Etold)).^2,'k','linewidth',2)title('输出光脉冲功率在时域和频域内的分布')xlabel('时间（ps）')ylabel('功率（W）')axis([20400inf])subplot(2,1,2)plot(2*pi*c./(fftshift(w)+w0)*10^9,(abs(f