
%solve i A_t + i/2 alpha A-1/2 beta2 A_xx +gamma |A|^2 =0, or
%        A_t  =-1/2 alpha A -i/2+ beta2 A_xx + i gamma |A|^2
% N=1 soliton
% N=3 soliton
% Chirp N=1 soliton
% Two interacting solitons
% For the analysis of the appropriate case, it should be uncommeted

clear all;
tic;
format long;
fid=fopen('SSFM.txt','wt');
%======Materiall parameters======
%alpha=0.01;                % Fibre loss (dB/km)
alpha=0;
alpha=alpha/4.34;            %Fibre loss (1/km)
beta2=-20e-24;               % GVD /2nd order disp (x^2/km)
gamma=1.317;                 % Fibre nonlinearity {1/W/km)

%===========Time parameter========
T=-512e-12:0.01e-12:511e-12;
deltat = 1.763;

%===== Input pulse and its parameter=====
T0 =7.0902e-12;                 %Initial pulse wdth (s)
%tau=T/T0;  
P0 = 0.30208;                % Input power (Watts);
A0 = (P0);
h=10d0;                       % step size (km)
L_D=((T0^2)/(abs(beta2)));   % Dispersion length in (km)
L_NL= 1/(gamma*P0);
N=1;% soliton number
A = A0*sech(T/T0);           % N= 1 soliton
plotting(:,1)=abs(A);
%write IC to data file
for k = 1: 102301
    if mod (k,100) == 0
        fprintf(fid,'%10.4f %10.4f %e\r\n', T(k)*10^(12), 0, plotting (k,1));
    end
end
z_plot=(0:h:((pi/2)*L_D));   % z-values to plot (km) tally with axis-t in temporal
z_0=pi/2*L_D;                % soliton period (km)
                       

%A = A0*sech(T/T0).*exp(i*pi*T/T0);    
%A = N^2*A0*sech(T/T0);       % N = 3 soliton
%C=1.6  ;                     % chirp parameter
%A=A0*sech(T/T0).*exp((-1i*C/2).*(T/T0).^2);   % N=1 soliton with chirp
%A = A0*sech(T/T0-1.6)+A0*sech(T/T0+1.6)  ;    % interaction soliton

                    % step size (km) tally with dt in temporal soliton
                            %  soliton
%n=1;                        % control stepping

L=max(size(A)); % tally with N in pseudo, in this case is 1024
         % step size (km)


Delta_omega= 1/L/deltat*2*pi;
omega=(-L/2:1:L/2-1)*Delta_omega;  % dimension 1*1024
M=ceil(z_0/h);
for m = 1:1:M                  % Start spatial evolution

  A = exp(h*i*gamma*abs(A).^2).*A;  % Solve nonlinear part of NLS
  c = fftshift(fft(A));  % Take Fourier transform
  %c = fftshift(fft(u))/sqrt(N); 
  %c=fft(u);

  %solve linear part
  c = exp(h*(-alpha/2+i*beta2/2*omega.^2)).*c;                   % Advance in Fourier time
  A = ifft(fftshift(c)); % Return to physical space 
 %u=ifft(ifftshift(c))*sqrt(N);
 %u=ifft(c);

 plot(T*1d12,abs(A));              % simulation 2d
 pause(0.0001)                      % simulation 2d off pause all
 plotting(:,m+1)=abs(A); %simpan untuk 2d slot, dimension punye viarable.
% n=n+1;
%write IC to data file
for k = 1: 102301
    if mod (k,100) == 0
        fprintf(fid,'%10.4f %10.4f %e\r\n', T(k)*10^(12),m*h, plotting (k,m+1));
    end
end

end 
fclose (fid);
toc;    %calculate computer time
cputime=toc;
z_plot=0:h:z_0+h;                    %nk tgk bpe dimension %create distance for plotting [0 1 2 3 4]
length(z_plot);                    % 5
for k =1:1:length(z_plot);         % Choosing 2d plots every 80th step [ 1 2 3 4 ]
    y =z_plot(k)*ones(size(T));      % Spread out along y-axis 1024 columns
    plot3(T*1d12,y,plotting(:,k),'LineWidth',1.5)
    grid on
    hold on
end
xlabel('time (ps)','FontSize', 12);
ylabel('distance(km)', 'FontSize',12);
zlabel('power (W)','FontSize',12);
set(gca,'FontSize',12) ; % size of tick marks on both axes
error=abs(max(plotting(:,1))-max(plotting(:,11)));
figure(2)
plot(T*1d12,plotting(:,1),'LineWidth',1.5);
title('z=0');
xlabel('time (ps)','FontSize', 14);
ylabel('power (w)', 'FontSize',14);
grid on
figure(3)
plot(T*1d12,plotting(:,11),'LineWidth',1.5);
title('z=1');
xlabel('time (ps)','FontSize', 14);
ylabel('power (w)', 'FontSize',14);
grid on
figure(4)
plot(T*1d12,plotting(:,21),'LineWidth',1.5);
title('z=2');
xlabel('time (ps)','FontSize', 14);
ylabel('power (w)', 'FontSize',14);
grid on
figure(5)
plot(T*1d12,plotting(:,31),'LineWidth',1.5);
title('z=3');
xlabel('time (ps)','FontSize', 14);
ylabel('power (w)', 'FontSize',14);
grid on
figure(6)
plot(T*1d12,plotting(:,41),'LineWidth',1.5);
title('z=4');
xlabel('time (ps)','FontSize', 14);
ylabel('power (w)', 'FontSize',14);
grid on






