www.pudn.com > NLS.rar > NLS.m, change:2012-03-20,size:2742b


% This code solves the NLS equation with the split-step method 
% idu/dz - sgn(beta2)/2 d^2u/d(tau)^2 + N^2*|u|^2*u = 0 
% Specify input parameters 
clear all;  
distance=10;   % fiber length 
beta2 =  -1;   % dispersion of fober  
N = 1;         % nonlinear parameter    
mshape = 1;    % pulse shape mshap=0 for sech, mshape>0 for super-Gausssian 
chirp0 = 0;    % input pulse chirp (default value) 
A = 300;   % input pulse power  
a = 1.263;   % 修正常数 
 
%%%%%%%%%%%     set simulation parameters 
nt = 1024; Tmax = 32;                % 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(叠代距离h) 
dtau = (2*Tmax)/nt;                  % step size in tau 
%%%%%%%%%    tau and omega arrays 
tau = (-nt/2:nt/2-1)*dtau;           % temporal grid 
omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]; % frequency grid 
%%%%%%%%%    Input Field profile 
if mshape==0  % 输入脉冲生成 
uu = sech(tau).*exp(-0.5i*chirp0*tau.^2); % soliton 
else                                      % super-Gaussian 
uu = exp(-0.5*(1+i*chirp0).*tau.^(2*mshape)); 
end 
 
%%%%%%%%% Plot input pulse shape and spectrum 
temp = fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi); % spectrum 
figure(1); 
subplot(2,1,1); 
plot(tau, abs(uu).^2, 'r','linewidth',2); hold on; % plot(tau, A*a^2*abs(uu).^2) 输出高斯脉冲时域 
axis([-5 5 0 inf]); 
xlabel('Normalized Time'); 
ylabel('Normalized Power'); 
subplot(2,1,2); 
plot(fftshift(omega)/(2*pi), abs(temp).^2, 'b','linewidth',2); % plot(fftshift(omega)/(2*pi), A*a^2*abs(temp).^2) 输出高斯脉冲频域 
hold on; 
axis([-.5 .5 0 inf]); 
xlabel('Normalized Frequency'); 
ylabel('Spectral Power'); 
%%%%%%%%%  store dispersive phase shifts to speedup code 
dispersion = exp(i*0.5*beta2*omega.^2*deltaz); % phase factor 
hhz = i*N^2*deltaz;                             % nonlinear phase factor 
%%%%%%%%     [ Beginning of MAIN Loop]    %%%%%%%%%%% 
%%%%%%%% scheme: 1/2N -> D -> 1/2N; first half step nonlinear 
temp = uu.*exp(abs(uu).^2.*hhz/2);             % note hhz/2 
for n=1:step_num 
f_temp = ifft(temp).*dispersion; 
uu = fft(f_temp); 
temp = uu.*exp(abs(uu).^2.*hhz); 
end 
uu = temp.*exp(-abs(uu).^2.*hhz/2);            % Final field 
temp = fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi); %Final spectrum 
%%%%%%%%   [ End of MAIN Loop ]          %%%%%%%%%%%%% 
%%%%%%%%    Plot output pulse shape and spectrum 
figure(2);   
subplot(2,1,1); 
plot(tau, abs(uu).^2 ,'c','linewidth',2); 
hold on; 
axis([-5 5 0 inf]); 
xlabel('Normalized Time'); 
ylabel('Normalized Power'); 
subplot(2,1,2); 
plot(fftshift(omega)/(2*pi), abs(temp).^2,'k','linewidth',2); 
hold on; 
axis([-.5 .5 0 inf]); 
xlabel('Normalized Frequency'); 
ylabel('Spectral Power');