clear all; close all;

% --- Define constants and initial condition
alpha1 = 0.065;
alpha2 = 0.012;
L = 50.0; % length of domain in x direction in cm
tmax = 1000.0;  % end time in seconds



tNodes = 10000.0;
xNodes = 200.0;
%%
dx = 1.0/xNodes;
dt = 1.0/tNodes;
%%
r1 = (alpha1*dt)/(dx^2);
r2 = 1 - 2*r1;

rStar = (alpha2*dt)/(dx^2);
r2Star = 1 - 2*rStar;
%%
sigma = 10; %FWHM in cm
% --- Loop over time steps
tArray = linspace(0,tmax,tNodes);
xArray = linspace(-L,L,xNodes);


f = @(x) (1/sqrt(2*pi*sigma.^2))*exp(-x.^2/(4*sigma.^2));
u = f(xArray);  % initial condition

%%%%%%%%%%%%%%%%%%%
xSize = length(xArray);
tSize = length(tArray);
uArray = zeros(tSize, xSize);


%Initializing u Array
uArray(1,:) = u(:);

%Creating laplacian matrix
%first, creating r-array for varying coefficient
rArray = ones(1,200);
rArray(1:89)    = r2;
rArray(90:94)  = r2Star;
rArray(95:105)  = r2;
rArray(106:110)  = r2Star;
rArray(111:200) = r2;

rArray = transpose(rArray);

rArray2 = ones(1,200);
rArray2(1:89)    = r1;
rArray2(90:94)  = rStar;
rArray2(95:105)  = r1;
rArray2(106:110)  = rStar;
rArray2(111:200) = r1;

diagR = eye(200);
i = 1;
while i <= 200
    diagR(i,i) = rArray(i);
    i = i + 1;
end

Atemp = diagR;

ASupDiag = zeros(200,200);
ASubDiag = zeros(200,200);


j = 1;
while j <= 100
    ASupDiag(j, j+1) = rArray2(j);
    ASubDiag(j + 1,j) = rArray2(j);
    j = j + 1;
end

% ASubDiag(201,:) = [];
% ASupDiag(:,201) = [];

A = Atemp + ASubDiag + ASupDiag;
A(1,1) = 1;
A(end,end) = 1;
A(end, end - 1) = 0;
A(1,2) = 0;

Prop = A;
%%%% Finished creating propagator


%%

for m=1:tNodes - 1
    uTemp = reshape(uArray(m,:),[],1);

    uArray(m + 1, :) = A*uTemp;
end


figure
colormap hsv
surf(xArray, tArray, uArray,'FaceColor','interp',...
   'EdgeColor','none',...
   'FaceLighting','gouraud')
ylim([0 100])
xlabel('x (cm)')
ylabel('t (seconds)')
zlabel('$\rho$', 'interpreter', 'latex')
camlight left


%% Integration

NormVec = uArray(1,:);
norm = sum(NormVec);
uArray = uArray/norm;


for j = 1:64000
    sumArray = reshape(uArray(j, 45:55), [], 1);
    sum1 = sum(sumArray);
    if sum1 < 0.01
        timestamp = tArray(j);
        break
    else
        continue
    end
end