| 查看: 1256 | 回复: 0 | ||
[求助]
用分步傅里叶法求解光纤中非线性薛定谔方程
|
|
不知这里哪位虫子有用分步傅里叶法求解光纤中非线性薛定谔方程的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,  =temp;end %for m=1:4 plot(fftshift(omega)/(2*pi),-20*(log10(abs(TEMP(m, ).^2)));% plot(fftshift(omega)/(2*pi),-20*(log10(abs(TEMP(m, ).^2)));%end |
回复此楼» 猜你喜欢
科研孤儿太难了
已经有9人回复
-
时间戳又变了
已经有15人回复
-
时间戳变了，能看出什么问题？
已经有23人回复
-
只有每年这种时候来逛逛小木虫
已经有22人回复
-
今天系统多次维护，明天很可能放榜！
已经有10人回复
-
放榜前的不淡定
已经有12人回复
-
言之凿凿，明天官网会通知放榜？
已经有16人回复
-
重要消息，中午系统在维护
已经有12人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有20人回复
-
今天基金会出结果吗？20260819
已经有15人回复
点击这里搜索更多相关资源
40