% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
% addpath(genpath('../DATABASE'));
% addpath(genpath('../TOOLBOXES'));

%% CHAPTER 13: simulation of chemical reactions

%% Problem 1.

close all
clear

%    k1   k2   k3
% U1<=>U2<=>U3<=>U4
%    k-1  k-2  k-3

k1 = 3; 
k2 = 7; 
k3 = 3; 
k_1 = 5; 
k_2 = 1; 
k_3 = 4;

% dU1/dt = -k1U1 + (k-1)U2 + 0*U3 + 0*U4
% dU2/dt = k1U1 [-(k-1) -(k2)]U2 + (k-2)U3 +0*U4
% dU3/dt = 0*U1 + k2U2 [-(k-2) -k3]U3 +(k-3)U4
% dU4/dt = 0*U1 + 0*U2 + k3U3 -(k-3)U4

% Rate matrix
row1 = [-k1         k_1         0           0  ];
row2 = [k1      -(k_1+k2)       k_2         0  ];
row3 = [0           k2      -(k_2+k3)       k_3];
row4 = [0           0           k3         -k_3];

K = [row1;row2;row3;row4]
[S,D] = eig(K)

u_vec_0 = [0.5 0.2 0.6 0.7]';

%% Exact solution
C = S\u_vec_0
dt = 0.01; % stepsize
end_time = 10; % last time point in seconds

D = diag(D);
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

u_vec = zeros(nsteps,4);

figure
xlim([0 10]);
ylim([0 1]);
xlabel('time (s)');
ylabel('Concentration');
box on
grid on
hold on

for i = 1:nsteps
   t = time_vec(i);
   u_vec(i,:) =  S*diag(exp(D*t))*C;
end

plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')

%% Trapezoidal method
I = eye(4);
G = (I - K*dt/2)\(I + K*dt/2);

u_vec = zeros(4,nsteps);
u_vec(:,1) = u_vec_0;
for i = 1:nsteps-1
    u_vec(:,i+1) = G*u_vec(:,i);
end

figure;plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

%% Midpoint method
I = eye(4);
G = (I + K*dt + K*K*dt*dt/2);
 
u_vec = zeros(4,nsteps);
u_vec(:,1) = u_vec_0;
for i = 1:nsteps-1
    u_vec(:,i+1) = G*u_vec(:,i);
end
 
figure;plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

%% Problem 2a.

% Run a small molecular dynamics simulation with 50 particles at 300 K with
% temperature control using the function 'simple_md' and be ready to
% explain the different steps involved.

trj = simple_md(50,0.1,300,0.0005,1000,10); 

%%  2b. Display the ensemble evolution in time

[~,~,nFrames] = size(trj);
figure
xlim([-4 4])
ylim([-3 3])
zlim([-3 3])

for i = 1:nFrames
    x = trj(1,:,i);
    y = trj(2,:,i);
    z = trj(3,:,i);
    plot3(x,y,z,'o','Linewidth',1.0,'MarkerEdgeColor','r',...
        'MarkerFaceColor','y')
    box('on');grid('on')
    ax = gca;
    ax.BoxStyle = 'full';
    xlabel('X coord')
    ylabel('Y coord')
    zlabel('Z coord')    
    drawnow update ; pause(.1)    
end

