分步傅里葉算法求解GNLSE

%光脈沖參數(shù)定義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%時(shí)頻窗口定義num=2^16;%時(shí)頻域內(nèi)的離散點(diǎn)數(shù)Ttotal=500*(2*T0);%時(shí)間窗口寬度，單位：sdeltaT=Ttotal/(num-1);%時(shí)間分辨率，單位：sT=(-num/2:num/2-1)*deltaT;%時(shí)間向量，單位：sdeltaw=2*pi/Ttotal;%角頻率分辨率，單位：rad/sw=deltaw*[(0:num/2-1)(-num/2:-1)];%角頻率向量，單位：rad/s%輸入脈沖時(shí)頻分析Et=sqrt(Ppeak)*sech(T/T0).*exp(-j*C/2*(T/T0).^2);%雙曲正割脈沖的時(shí)域電場分布Ew=ifft(Et);%雙曲正割脈沖的頻域電場分布%輸入脈沖時(shí)頻特性圖形顯示figure(1)subplot(2,1,1)plot(T*10^12,(abs(Et)).^2,'k','linewidth',2)title('輸入光脈沖功率在時(shí)域和頻域內(nèi)的分布')xlabel('時(shí)間（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])%光纖參數(shù)的定義beta2=-3.8237e-27;%GVD參數(shù)，單位s2/mbeta3=6.3147e-41;%TOD參數(shù)，單位s3/mgama=10e-3;%非線性系數(shù)，單位1/(W.m)Tr=3e-15;%拉曼響應(yīng)函數(shù)關(guān)于時(shí)間的一階矩，單位：salfa=1.5e-3/4.343;%光纖損耗系數(shù)，單位：1/mL=1000;%光纖長度，單位：mh=1;%空間離散步長，單位：mN=L/h+1;%空間離散點(diǎn)數(shù)%計(jì)算中所用到的矩陣定義Etold=zeros(1,num);%用于存儲(chǔ)電場隨時(shí)間的分布Ewold=zeros(1,num);%用于存儲(chǔ)電場隨角頻率的分布Raman=zeros(1,num);%用于存儲(chǔ)時(shí)域內(nèi)電場模的平方對時(shí)間的導(dǎo)數(shù)Self=zeros(1,num);%用于存儲(chǔ)時(shí)域內(nèi)電場模的平方與電場乘積對時(shí)間的導(dǎo)數(shù)%光纖輸入端（z=0）的電場時(shí)頻分布Etold=Et;%z=0處時(shí)域電場分布Ewold=ifft(Etold);%z=0處角頻率域電場分布%前h/2色散影響的計(jì)算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;%計(jì)算第一個(gè)半步長色散的影響Etold=fft(Ewold);%N-2段周期性h非線性和h色散影響的計(jì)算fornumber=1:N-2%求h長度內(nèi)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;%逆傅里葉變換求時(shí)域電場分布Etold=fft(Ewold);end%最后一個(gè)h非線性和h/2色散影響的計(jì)算%求最后一個(gè)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);%求最后一個(gè)h/2長度中色散影響Ewold=Ewold.*dispersion_half;Etold=fft(Ewold);%輸出脈沖時(shí)頻特性圖形顯示figure(2)subplot(2,1,1)plot(T*10^12,(abs(Etold)).^2,'k','linewidth',2)title('輸出光脈沖功率在時(shí)域和頻域內(nèi)的分布')xlabel('時(shí)間（ps）')ylabel('功率（W）')axis([20400inf])subplot(2,1,2)plot(2*pi*c./(fftshift(w)+w0)*10^9,(abs(f