||
如下Matlab 程序计算 论文Youfa Wang and Wenfeng Wang, Study of ultrafast pulse coupling dynamics considering retarded
nonlinear response and self-steepening effects, Joutnal of lightwave Technology, Vol.24, No.2,
pp1041-1047,中的 图2。 做一点简单改动,文中其他图,也可以计算。
% nonlinear equation including nonlinear delay, self-steepening, and higher order dispersion
% Youfa Wang and Wenfeng Wang, Study of ultrafast pulse coupling dynamics considering retarded
% nonlinear response and self-steepening effects, Joutnal of lightwave Technology, Vol.24, No.2,
% pp1041-1047,
MN3=0;
lcd=0.5; %lcd>0 corresponding to beita2<0, lcd<0 corresponds
%lcd=1000;
ldd=-200*lcd;
lmd=-0.00;
fr=0.18;
st=0.028; %st=1/(w0*t0)
%st=0;
J=400;
N =1024; % Number of Fourier modes (Time domain sampling points)
M3=6000;
dz =1.570796326/M3; % space step, make sure nonlinear<0.05
T =20; % normalized time, t=T*T0, it can't be too big ot too small, it affect accuracy
delt=-0.0; % normal mismatch=(b1-b2)/2c
T0=28.4;
dt = T/N; % time step
n = [-N/2:1:N/2-1]'; % Indices
t = n.*dt;
%----Ranman response in fiber---------------------------
n1=[0:1:N/2]'; % the response is 0, if t1<0. it is not easy to use t. the result in funcation no difference
t1=n1*dt;
ht=(32*32+12.2*12.2)/(32*32*12.2)*T0*exp(-t1*T0/32).*sin(t1*T0/12.2);
%----------------display convolution-------------
y=[-N/2:1:N-1]';
t2=dt*y;
%--------------------------------------------------------------
ww = 4*n.*n*pi*pi/T/T; % Square of frequency. Note i^=-1.
w=2*pi*n./T;
www=w.*ww;
g1=i*(delt/T0)*w-i*ww./(2*lcd)-i*www./(6*ldd);
g2=i*(delt/T0)*w-i*ww./(2*lcd)-i*www./(6*ldd); %w=2*pi*f*n./N, f=1/dt=N/T,so w=2*pi*n./TP=0;
p=1.2;
s10=p.*sech(t);
s1=s10;
s20=0*s10;
s2=s20;
S1=s1.*0; S2=s1.*0;SC1=s1.*0; SC2=s1.*0;
p10=dt*(sum(abs(s10').*abs(s10'))-0.5*(abs(s10(N,1)*s10(N,1))+abs(s10(1,1)*s10(1,1)))); %energy in waveguide 1
for m3 = 1:1:M3 % Start space evolution
s1squre=abs(s1.*s1);
s1cubic=s1squre.*s1;
s2squre=abs(s2.*s2);
s2cubic=s2squre.*s2;
chta=conv(ht,s1squre);
cht1=dt*chta(1:N);
chtb=conv(ht,s2squre);
cht2=dt*chtb(1:N);
s1(1)=s1(1)+4*(1-fr)*i*dz*s1cubic(1)-4*(1-fr)*st*dz/dt*s1cubic(1)+4*fr*i*dz*s1(1).*cht1(1)-...
4*fr*st/dz*s1(1).*cht1(1);
s1(2:N)=s1(2:N)+4*(1-fr)*i*dz*s1cubic(2:N)-4*(1-fr)*st*dz/dt*(s1cubic(2:N)-...
s1cubic(1:N-1))+4*fr*i*dz*s1(2:N).*cht1(2:N)-4*fr*st*dz/dt*(s1(2:N).*cht1(2:N)-s1(1:N-1).*cht1(1:N-1));
s2(1)=s2(1)+4*(1-fr)*i*dz*s2cubic(1)-4*(1-fr)*st*dz/dt*s2cubic(1)+4*fr*i*dz*s2(1).*cht2(1)-...
4*fr*st/dz*s2(1).*cht2(1);
s2(2:N)=s2(2:N)+4*(1-fr)*i*dz*s2cubic(2:N)-4*(1-fr)*st*dz/dt*(s2cubic(2:N)-...
s2cubic(1:N-1))+4*fr*i*dz*s2(2:N).*cht2(2:N)-4*fr*st*dz/dt*(s2(2:N).*cht2(2:N)-s2(1:N-1).*cht2(1:N-1));
sca1 = fftshift(fft(s1)); % Take Fourier transform
sca2 = fftshift(fft(s2));
sc2=exp(g2.*dz).*(sca2+i*(1+lmd*w).*sca1.*dz);
sc1=exp(g1.*dz).*(sca1+i*(1+lmd*w).*sca2.*dz); % frequency domain phase shift
s2 = ifft(fftshift(sc2)); % Return to physical space
s1 = ifft(fftshift(sc1));
% if rem(m3,J) == 0 % Save output every J steps.
% S1 = [S1 s1]; % put solutions in U array
% S2=[S2 s2];
% SC1 = [SC1 sc1]; % put solutions in U array
% SC2=[SC2 sc2];
% MN3=[MN3 m3];
% z3=dz*MN3'; % output location
%end
m3,a=max(abs(cht1)),
end
p1=dt*(sum(abs(s1').*abs(s1'))-0.5*(abs(s1(N,1)*s1(N,1))+abs(s1(1,1)*s1(1,1))));
p2=dt*(sum(abs(s2').*abs(s2'))-0.5*(abs(s2(N,1)*s2(N,1))+abs(s2(1,1)*s2(1,1))));
p1+p2-p10,
%figure(3)
%waterfall(w',z3',abs(SC1').*abs(SC1')) %t' is 1xn, z' is 1xm, and U1' is mxn
%figure(4)
%waterfall(w',z3',abs(SC2').*abs(SC2')) %t' is 1xn, z' is 1xm, and U1' is mxn
figure(1)
plot(t,abs(s1'.*s1'), t,abs(s2'.*s2'),'.');
figure(2)
plot(t,abs(s2'.*s2'),'.');
%figure(2)
%plot(w,abs(sc1'.*sc1'), w,abs(sc2'.*sc2'));
%-----------display convolution---
%plot(t2,chta)
%------------------------------------
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 07:06
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社