clear all;
close all;

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

% Set stepsize for stepsize (in visualization) :-)
K = 3;

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

% Modify initial value such that we can consider it the DFT of a H^1
% function
D = i*diag([-N:(N-1)]);
Dsquare = D*D;
mult = (-Dsquare + eye(2*N))^(-1/2);

iv = mult*iv;
u0 = iv/norm(iv);

% Calculate norms
u0H1norm = norm((-Dsquare + eye(2*N))^(1/2)*u0);
u0H2norm = norm((-Dsquare + eye(2*N))*u0);

disp(' ');
disp(['----- First initial value:']);

fprintf('H^1 norm of first initial value: %d\nH^2 norm of first initial value: %d\n', round(u0H1norm), round(u0H2norm));

splittingfigure = figure();

for k = [0:3*K] 
    disp([' Now running with stepsize h = 10^(-', num2str(k) ,'/' , num2str(K) ,')...']);
    % Calculate Marchuk-Strang splitting with stepsize h = 10^(-k/K)
    [sol,x,y,z] = schroedinger_operatorsplitting(N,10^(-k/K),1,2,u0); 
    % Plot results logarithmically
    loglog(10^(-k/K),x,'kx',10^(-k/K),y,'rx',10^(-k/K),z,'bx'); 
    if (k==0)
        hold on;
    else
        loglog([10^((-k+1)/K),10^(-(k/K))],[a,x],'k-',[10^((-k+1)/K),10^((-k/K))],[b,y],'r--',[10^((-k+1)/K),10^((-k/K))],[c,z],'b-.'); 
    end
    a = x; b = y; c = z;    
end

% Modify initial value further such that we can see it as the DFT of a
% H^2-function
iv = mult*iv;
u0 = iv/norm(iv);

% Calculate relevant norms
u0H1norm = norm((-Dsquare + eye(2*N))^(1/2)*u0);
u0H2norm = norm((-Dsquare + eye(2*N))*u0);

disp(' ');
disp(['----- Second initial value:']);

fprintf('H^1 norm of second initial value: %d\nH^2 norm of second initial value: %d\n', round(u0H1norm), round(u0H2norm));

for k = [0:3*K] 
    disp([' Now running with stepsize h = 10^(-', num2str(k) ,'/' , num2str(K) ,')...']);
    % Calculate Marchuk-Strang splitting with stepsize h = 10^(-k/K)
    [sol,x,y,z] = schroedinger_operatorsplitting(N,10^(-k/K),1,2,u0); 
    loglog(10^(-k/K),x,'ko',10^(-k/K),y,'ro',10^(-k/K),z,'bo'); 
    if (k~=0) 
        loglog([10^((-k+1)/K),10^(-k/K)],[a,x],'k-',[10^((-k+1)/K),10^(-k/K)],[b,y],'r--',[10^((-k+1)/K),10^(-k/K)],[c,z],'b-.'); 
    end
    a = x; b = y; c = z;   
end

% Bring the figure to the top! :)
figure(splittingfigure);