%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Michael Pokojovy                             % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

rho = 19300;
E   = 80E9;

c = sqrt(E/rho);

L = 0.3;
c = sqrt(E/rho);

L = 0.3;

T = 6*L/c;

mu = 2*0.25E-4*c;
mu = 2*c;

U0 = @(x) L/2 - abs(x - L/2);

N = 1000;
M = 10000;

n = N - 2;

xlat = linspace(0, L, N);
dx = xlat(2) - xlat(1);
xlat_int = xlat(2:end-1);

tlat = linspace(0, T, M);
dt = tlat(2) - tlat(1);

V0 = [U0(xlat_int)';
      zeros(length(xlat_int), 1)];

V = zeros(length(xlat_int)*2, length(tlat));
V(:, 1) = V0;

D = sparse(1:n, 1:n,  -2*ones(1,n),n,n);
E = sparse(2:n, 1:n-1, ones(1,n-1),n,n);

A = 1/(dx^2)*(E+D+E');

for i = 2:length(tlat)
    B = [sparse(n, n) speye(n);
         c^2*A  -mu*speye(n)];
    
    V(:, i) = (speye(2*n) - 0.5*dt*B)\((speye(2*n) + 0.5*dt*B)*V(:, i - 1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);

set(gcf, 'PaperUnits', 'centimeters');
xSize = 20; ySize = 10;
xLeft = (21 - xSize)/2; yTop = (30 - ySize)/2;
set(gcf,'PaperPosition', [xLeft yTop xSize ySize]);
set(gcf,'Position',[0 0 xSize*50 ySize*50]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
It = floor(((1:50) - 1)*(M - 1)/(50 - 1)) + 1;
Ix = floor(((1:50) - 1)*(N - 1)/(50 - 1)) + 1;

[T, X] = meshgrid(tlat(It), xlat(Ix));

U = [zeros(1, length(tlat)); 
     V(1:length(xlat_int), :); 
     zeros(1, length(tlat))];

colormap hot;
surf(X, T, U(Ix, It));

view([39 38]);

grid on;

xlabel('x [m]');
ylabel('t [s]');
zlabel('u(t, x) [m]');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);

set(gcf, 'PaperUnits', 'centimeters');
xSize = 20; ySize = 10;
xLeft = (21 - xSize)/2; yTop = (30 - ySize)/2;
set(gcf,'PaperPosition', [xLeft yTop xSize ySize]);
set(gcf,'Position',[0 0 xSize*50 ySize*50]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on;

t = title('L\"osung des ARWP der gezupften Saite', 'interpreter', 'latex');
set(t, 'FontSize', 16);

k = 0;
for i = 1:3
for j = 1:3
    k = k + 1;
    subplot(3, 3, k);
    
    tt = floor((k - 1)*(M - 1)/8) + 1;
    
    plot(xlat, [0 (V(1:n, tt))' 0]);
    
    axis([0 L -L/2 L/2]);
    
    h = title(['$t = $ ' num2str(tlat(tt))], 'interpreter', 'latex');
    set(h, 'FontSize', 16);
    
    zlabel('u(t, x)');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%