24小时热门版块排行榜

返回列表

lmhdzxx

新虫 (正式写手)

应助: 1 (幼儿园)
金币: 1444.2
散金: 11
帖子: 630
在线: 153小时
虫号: 1693164
注册: 2012-03-15
专业: 激光

[求助] 用分步傅里叶法求解光纤中非线性薛定谔方程

不知这里哪位虫子有用分步傅里叶法求解光纤中非线性薛定谔方程的MATLAB程序啊？我用Agrawal那本非线性光纤光学后面的程序加以改动，但是得出的图效果不理想，有没有哪位同学有程序啊？真心求交流。。。。。
附程序：
%传输距离为1到4，其他参数不变，最后得出不同传输距离下频域光脉冲形状的比较
clear all;
for m=1:4
%---Input Field profile
distance = 10*m; %传输距离
beta2 =-1;%二阶色散参数
chirp = 0; % 初始啁啾参数
peakpower=1;%峰值功率
T0=1;%初始脉宽
r=1  ;  %非线性系数
N=sqrt(r*peakpower*T0^2) ;
alpha=1;%光纤损耗
%---set simulation parameters
nt =1024; Tmax = 10; % FFT points and window size,设置采样点和时间窗口
step_num = round(20*distance*N^2) % No. of z steps to
deltaz = distance/step_num; % step size in z
dtau = (2*Tmax)/nt; % step size in tau
%---tau and omega arrays
tau = (-nt/2:nt/2-1)*dtau; % temporal grid
omega = (pi/Tmax)*[(-nt/2:nt/2-1)]; % frequency grid
%---store dispersive phase shifts to speedup code
dispersion = exp((1i*0.5*beta2*omega.^2-alpha/2)*deltaz);% phase factor
hhz = 1i*N^2*deltaz; % nonlinear phase factor
uu =sqrt(peakpower)*exp(-0.5*tau.^2.*(1+1i*chirp)/T0^2);  % gaussian
temp = uu.*exp((abs(uu).^2).*hhz/2); % note hhz/2
subplot(2,1,1)
if m==1
   plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'r')
   hold on
else if m==2
   plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'y')
   hold on
   else if m==3
         plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'b')
         hold on
      else plot(fftshift(omega)/(2*pi),-log10(abs(temp).^2),'k')
      end
   end
end
%*********[ Beginning of MAIN Loop]***********
% scheme: 1/2N -> D -> 1/2N; first half step nonlinear
for n=1:step_num
f_temp = ifft(temp).*dispersion;
uu = fft(f_temp);
temp = uu.*exp((abs(uu).^2).*hhz/2);
end
uu = temp.*exp(-(abs(uu).^2).*hhz/2);% Final field
temp = fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi);%Final spectrum
subplot(2,1,2)
if m==1
   plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'r')
   hold on
else if m==2
   plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'y')
   hold on
   else if m==3
         plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'b')
         hold on
      else plot(fftshift(omega)/(2*pi),log10(abs(temp)).^2,'k')
      end
   end
end

%***************[ End of MAIN Loop ]**************
%----Plot output pulse shape and spectrum
%TEMP(m,