agrawalAppendixB.m

   1 distance = input('Enter fiber length in L_D ');
   2 beta2 = input('dispersion: 1 for normal, -1 for anomalous ');
   3 N = input('Nonlinear parameter N = '); % Soliton order
   4 mshape = input('m = 0 for sech, m > 0 for super-Gaussian ');
   5 chirp0 = 0;
   6
   7 % set simulation parameters
   8 nt = 1024; Tmax = 32; % FFT points and window size
   9 step_num = round(20 * distance * N^2); % No. of z steps
  10 deltaz = distance/step_num; % step size in z
  11 dtau = (2*Tmax) / nt; % step size in tau
  12
  13 %% tau and omega arrays
  14 tau = (-nt/2 : nt/2-1) * dtau; % temporal grid
  15 omega = (pi/Tmax) * [(0:nt/2-1) (-nt/2:-1)]; % freq grid
  16
  17 if mshape == 0
  18   uu = sech(tau) .* exp(-0.5j * chirp0 * tau.^2);
  19 else
  20   uu = exp(-0.5 * (1 + 1j * chirp0) .* tau.^(2 * mshape));
  21 end
  22
  23 %% plot input pulse shape and spectrum
  24 temp = fftshift(ifft(uu)) .* (nt * dtau) / sqrt(2 * pi); % spectrum
  25 figure(1); clf; subplot(2,1,1);
  26 plot(tau, abs(uu).^2, '--k'); hold on;
  27 axis([-5 5 0 inf]);
  28 xlabel('Normalized Time');
  29 ylabel('Normalized Power');
  30 title('Input and Output pulse shape and spectrum');
  31
  32 subplot(2, 1, 2);
  33 plot(fftshift(omega)/(2*pi), abs(temp) .^ 2, '--k'); hold on;
  34 axis([-.5 .5 0 inf]);
  35 xlabel('Normaized freq');
  36 ylabel('spectral power');
  37
  38 %% store dispersive phase shifts to speed up code
  39 dispersion = exp(0.5j * beta2 * omega.^2 * deltaz); % [hase factor
  40 hhz = 1j * N^2 * deltaz;
  41
  42 %% begin main loop
  43 %% N/2 -> D -> N/2 first half step nonlinear
  44 temp = uu .* exp(abs(uu) .^ 2 .* hhz / 2);
  45 for n = 1 : step_num
  46   f_temp = ifft(temp) .* dispersion;
  47   uu = fft(f_temp);
  48   temp = uu .* exp(abs(uu) .^ 2 .* hhz);
  49 end
  50 uu = temp .* exp(-abs(uu) .^ 2 .* hhz); % final field
  51 temp = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2 * pi); % final spectrum
  52 %% end of main loop
  53
  54 %% plot output pulse shape and spectrum
  55 subplot(2,1,1);
  56 plot(tau, abs(uu).^2, '-k');
  57 subplot(2,1,2);
  58 plot(fftshift(omega) / (2*pi), abs(temp).^2, '-k');