clear all;
close all;

% Set 'grid' for space discretization
N = 256;

% Initial value initialization
iv = -100+200.*rand(2*N,1);

D = i*diag([-N:(N-1)]);
Dsquare = D*D;
mult = (-Dsquare + eye(2*N))^(-1/2);

iv = mult*iv;
% Choose your destiny: regular or not? Uncomment for more regularity..
% iv = mult*iv;

u0 = iv/norm(iv);

a = pi;
x = zeros(2*N,1);
f = zeros(2*N,1);
sigma=1;
for k=1:2*N
    x(k) = (-N-1+k)*a/N;
    f(k) = 1/sqrt(2*pi*sigma^2) * exp((-x(k)^2)/(2*sigma^2));
end

f = fft(f);
% Use Gaussian distribution as initial value
% u0 = f/norm(f);

scrsz = get(0,'ScreenSize');
fig = figure('OuterPosition',[50 (scrsz(4)/2-50) 2*scrsz(3)/3 scrsz(4)/2]);

disp(['Calculating Marchuk-Strang splitting with stepsize 10^(-3)']);
[sol] = schroedinger_operatorsplitting(N,10^(-3),1,2,u0);
[sol2] = schroedinger_operatorsplitting_nopotential(N,10^(-3),1,2,u0);
%norm(sol-sol2)

g = zeros(length(x),1);
n = length(sol(1,:));

% For each timestep..
for k = 1:n
    % .. calculate the complete solution
    g = 2*N/(sqrt(2*pi))*ifft(sol(:,k));
    g2 = 2*N/(sqrt(2*pi))*ifft(sol2(:,k));
    
    subplot(2,1,1);    
    plot(x,abs(g),'k');        
    xlabel('x');
    ylabel('u(t,x)');
    title(['Time step ' , num2str(k-1) , ' of ' , num2str(n-1) , ' - Potential: 1-cos(x)']);
    axis([-pi pi 0 1]);
    
    subplot(2,1,2)
    plot(x,abs(g2),'b');
    xlabel('x');
    ylabel('u(t,x)');
    title(['Potential: 0']);
    axis([-pi pi 0 1]);
    
%     subplot(2,1,2);
%     plot(k,sum(abs(g).^2)*(pi/N));       
%     title('Overall probability');
%     axis([1 1001 0.8 1.2]);
%     hold on;
    
    pause(0.01);
end