% 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

% direction plot for du/dt = ku

% Here we plot the slopes at the meshpoints: we choose a range of times
% (from 0 to 5) and a range of u's (from 0 to 1), and we calculate du/dt =
% ku for each t-u combination.
k = -0.5;
[t,u] = meshgrid(0:5,0:0.1:1)
figure;plot(t,u,'b.')
xlabel('time  ')
ylabel('u(t)  ')
xlim([-1,6])
ylim([-0.1,1.1])
hold on

du = k.*u
dt = ones(11,6)
quiver(t,u,dt,du,0.2,'r')

% Here we plot 6 different curves for different values of u(0).
t = [-1:0.1:6];
u0 = [0 0.1 0.25 0.5 1 2 4 8];

[t,u0] = meshgrid(t,u0)
u = u0.*exp(k*t)
plot(t(1,:),u',':b')

%% Analysis of 1st order reaction kinetics.

close all
clear, clc

%
% Let's consider a typical equilibrium:
%    k1   k2   
% U1<=>U2<=>U3
%    k-1  k-2  
%
k1 = .8; 
k2 = .4; 
k_1 = .3; 
k_2 = .5; 

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

% The rate matrix is:
row1 = [-k1         k_1         0     ];
row2 = [k1      -(k_1+k2)       k_2   ];
row3 = [0           k2         -(k_2) ];

K = [row1;row2;row3]


% The exact solution is:
dt = 0.3; % stepsize
end_time = 10; % last time point in seconds

[S,D] = eig(K);
D = diag(D);
t_vec = [0:dt:end_time];
nsteps = length(t_vec);

% Here we set up a loop:
u_vec = zeros(nsteps,3);
% Initial values
u_vec_0 = [100;0;0];

% Coordinate of u_vec_0 in S space;
C = S\u_vec_0;

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

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

   % Here we plot progressively   
   plot(t_vec(1:i),u_vec(1:i,1),'-xg',t_vec(1:i),u_vec(1:i,2),...
       '-xb',t_vec(1:i),u_vec(1:i,3),'-xr')
   drawnow expose; % pause(.05)
end

legend('u1','u2','u3')

% plot(t_vec,u_vec(:,1),'-xg',t_vec,u_vec(:,2),'-xb',t_vec,u_vec(:,3),'-xr')
% legend('u1','u2','u3')
% xlabel('time (s)');
% ylabel('Concentration');

%%
% Alternatively, representing directly the initial values as a linear
% combination of the eigenvectors

% Again we can run a loop to calculate values at all intermediate times:

% 1. find the coordinates of [u]_0 in the eigenvector space:
C = S\u_vec_0;

% 2. write [u]_0 as a linear combination of the eigenvectors of K:
S1 = C(1)*S(:,1);
S2 = C(2)*S(:,2);
S3 = C(3)*S(:,3);
S1 + S2 + S3

% Faster , but less transparent
Snm = S*diag(C);
sum(Snm,2)

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

for i = 1:nsteps
   t = t_vec(i);
   % 3. multiply each eigenvector component by the exponent of the product
   % of each eigenvalue by time:
   % D1 = D(1);
   % D2 = D(2);
   % D3 = D(3);
   % u_vec(i,:) =  exp(D1*t)*S1 + exp(D2*t)*S2 + exp(D3*t)*S3;
   
   % Faster , but less transparent
   u_vec(i,:) =  Snm*exp(D*t);
   
   % Here we plot progressively
   plot(t_vec(1:i),u_vec(1:i,1),'-og',t_vec(1:i),u_vec(1:i,2),'-ob',...
       t_vec(1:i),u_vec(1:i,3),'-or')
   drawnow expose; % pause(.1)   
end

legend('u1','u2','u3')

% plot(t_vec,u_vec(:,1),'-og',t_vec,u_vec(:,2),'-ob',t_vec,u_vec(:,3),'-or')
% legend('u1','u2','u3')
% xlabel('time (s)');
% ylabel('Concentration');

%% Direct solution using the exponential of the rate matrix
% Of course we can also solve directly for dU/dt = A*U with U(t) =
% U(o)*exp(At), but in this case the eigen decomposition of A is carried
% out at each cycle of the loop in order to calculate the exponent, and
% thus the operation is much slower.

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

for i = 1:nsteps
   t = t_vec(i);
   u_vec(i,:) = expm(K*t)*u_vec_0;
   
   % Here we plot progressively
   plot(t_vec(1:i),u_vec(1:i,1),'-sg',t_vec(1:i),u_vec(1:i,2),'-sb',...
       t_vec(1:i),u_vec(1:i,3),'-sr')
   drawnow expose; % pause(.1)   
end

legend('u1','u2','u3')

% plot(t_vec,u_vec(:,1),'-og',t_vec,u_vec(:,2),'-ob',t_vec,u_vec(:,3),'-or')
% legend('u1','u2','u3')
% xlabel('time (s)');
% ylabel('Concentration');

%%
% Equilibrium at infinite time:

ss = C(3)*S(:,3)

% We would obtain the same steady state also if we started from a different
% initial concentration of the chemical species as long as the sum of all
% the concentrations is the same.

u_vec_0 = [30;20;50];

% Coordinate of u_vec_0 in S space;
C = S\u_vec_0;

% Steady_state at infinite time:
ss = C(3)*S(:,3)

% We check that the equilibria are what we expected from the equilibrium
% constants:
%
%    k1     k2   
% U1 <=> U2 <=> U3
%    k-1    k-2 
%

k1/k_1
ss(2)/ss(1)
k2/k_2
ss(3)/ss(2)

%%
close all
clear

% Let's consider a more complex equilibrium:
%    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

% The rate matrix is:
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)

% We set the initial value 
u_vec_0 = [0.5 0.2 0.6 0.7]';

% Coordinate of u_vec_0 in S space. We solve the system of linear equations
% using the 'backslash' operator. Again, C is the representation of the
% initial value vector in the eigenvector basis.
C = S\u_vec_0

% We check the 4 components sum up to u_vec_0
S*C

% Reaction time course
dt = 0.05; % stepsize
end_time = 5; % last time point in seconds

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

% Here we set up a loop:
u_vec = zeros(nsteps,4);

figure
xlim([0 end_time]);
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;

   % Here we plot progressively   
   plot(time_vec(1:i),u_vec(1:i,:)); hold on
   drawnow expose; pause(.01)
end

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

% At infinite time the component corresponding to the 3 negative
% eigenvalues will vanish and the remaining components will determine the 
% steady-state :

steady_state = S(:,4)*C(4)

% We check the equilibrium is what we expected from the equilibrium
% constants:
%
%    k1   k2   k3
% U1<=>U2<=>U3<=>U4
%    k-1  k-2  k-3
%
k1/k_1
steady_state(2)/(steady_state(1))
k2/k_2
steady_state(3)/(steady_state(2))
k3/k_3
steady_state(4)/(steady_state(3))

%% Finite difference solutions
%% Forward Euler
close all

dt = 0.05; % stepsize
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

G = dt*K + eye(4);

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')
ylim([0 1]);
grid on

%% Backward Euler

dt = 0.05; % stepsize
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

G = inv(eye(4) - dt*K);

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')
ylim([0 1]);
grid on

%% Stability of Euler methods

%% Forward Euler
close all

[~,D] = eig(K);

Forward_Euler_Stability = figure;

dt = 0.01;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

G = dt*K + eye(4);
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
 
subplot(3,1,1) 
plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

dtD = dt*diag(D)
str = ['\DeltatD = ' num2str(min(dtD),'%7.4f\n')]
text(2,0.5,str)

dt = 0.1;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);
G = dt*K + eye(4);
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
 
subplot(3,1,2) 
plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

dtD = dt*diag(D)
str = ['\DeltatD = ' num2str(min(dtD),'%7.4f\n')]
text(2,0.5,str)

dt = 0.2;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);
G = dt*K + eye(4);
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
 
subplot(3,1,3) 
plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

dtD = dt*diag(D)
str = ['\DeltatD = ' num2str(min(dtD),'%7.4f\n')]
text(2,0.5e5,str)

%% Backward Euler
close all
Backward_Euler_Stability = figure;

dt = 0.1;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

G = inv(eye(4) - dt*K);
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
 
subplot(3,1,1) 
plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

dt = 0.2;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);
G = inv(eye(4) - dt*K);
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
 
subplot(3,1,2) 
plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

dt = 0.5;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);
G = inv(eye(4) - dt*K);
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
 
subplot(3,1,3) 
plot(time_vec,u_vec)
legend('U1','U2','U3','U4','Location','best')
grid on

%% Trapezoidal method
close all

dt = 0.05;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

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')
ylim([0 1]);
grid on

[S,D] = eig(G)
[~,l_max] = max(diag(D))
l_ind = true(4,1);
l_ind(l_max) = false;
D(l_ind,l_ind) = 0;
Ginf = S*D/S
u_vec_inf = Ginf*u_vec_0

%% Midpoint method
close all

dt = 0.05;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

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')
ylim([0 1]);
grid on

[S,D] = eig(G);
[~,l_max] = max(diag(D));
l_ind = true(4,1);
l_ind(l_max) = false;
D(l_ind,l_ind) = 0;
Ginf = S*D/S;
u_vec_inf = Ginf*u_vec_0

%% 4th order Runge-Kutta (RK4) - original version
close all

dt = 0.05;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

I = eye(4);

P = I + K*dt/2;
% K1 = K*(P*u_vec(:,i));
K1 = K;
Q = I + K1*P*dt/2;
% K2 = K*Q*u_vec(:,i); 
K2 = K;
R = I + K2*Q*dt;
% K3 = K*R*u_vec(:,i);
K3 = K;
G = I + (dt/6)*(K + 2*K1*P + 2*K2*Q + K3*R);
% u_vec(:,i+1) = G*u_vec(:,i);

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

[S,D] = eig(G);
[~,l_max] = max(diag(D));
l_ind = true(4,1);
l_ind(l_max) = false;
D(l_ind,l_ind) = 0;
Ginf = S*D/S;
u_vec_inf = Ginf*u_vec_0

%% 4th order Runge-Kutta (RK4) - rearranged 
close all

dt = 0.05;
end_time = 5; % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

I = eye(4);
K1 = K;
P = I + K*dt/2;
K2 = K*P;
Q = I + K*P*dt/2;
K3 = K*Q;
R = I + K*Q*dt;
K4 = K*R;
G = I + dt*(K1 + 2*K2 + 2*K3 + K4)/6;

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

[S,D] = eig(G);
[~,l_max] = max(diag(D));
l_ind = true(4,1);
l_ind(l_max) = false;
D(l_ind,l_ind) = 0;
Ginf = S*D/S;
u_vec_inf = Ginf*u_vec_0

%% ODE45,15s,23s
Fp2 = @(t,u) K*u;
tspan = [0 10]
[time_vec,u_vec] = ode45(Fp2,tspan,u_vec_0);
figure;plot(time_vec,u_vec,'-x');legend('U1','U2','U3','U4','Location','best');
title('ODE45');grid on
[time_vec,u_vec] = ode15s(Fp2,tspan,u_vec_0);
figure;plot(time_vec,u_vec,'-x');legend('U1','U2','U3','U4','Location','best');
title('ODE15s');grid on
[time_vec,u_vec] = ode23s(Fp2,tspan,u_vec_0);
figure;plot(time_vec,u_vec,'-x');legend('U1','U2','U3','U4','Location','best');
title('ODE23s');grid on

%% A simple kinetic solver
close all
clear

k1 = 5; 
k_1 = 2; 
k2 = 7; 
k_2 = 1; 

Uo = [2.0 2.0 0.0 0.0]';
U1 = Uo(1);
U2 = Uo(2);

row1 = [-k1*U2         0         k_1           0  ];
row2 = [0           -k1*U1       k_1           0  ];
row3 = [k1*U2          0      -(k_1+k2)        k_2];
row4 = [0              0         k2           -k_2];

K = [row1;row2;row3;row4];

time_step = 0.01;
time_vec = [0:time_step:1];
nsteps = length(time_vec);
U = zeros(4,nsteps);
U(:,1) = Uo;

for i = 2:nsteps
    K = [row1;row2;row3;row4];
    [S,D] = eig(K);
    D = diag(D);    
    Dt = D*time_step;
    eDt = diag(exp(Dt));
    U(:,i) = (S*eDt/S)*U(:,i-1);
    U1 = U(1,i);
    U2 = U(2,i);
    row1 = [-k1*U2         0         k_1           0  ];
    row2 = [0           -k1*U1       k_1           0  ];
    row3 = [k1*U2          0      -(k_1+k2)        k_2]; 
end

figure;plot(time_vec,U,'-x')
legend('U1','U2','U3','U4','Location','best')
title('Eigen solution');grid on

%% A simple kinetic solver using RK4
clear

k1 = 5; 
k_1 = 2; 
k2 = 7; 
k_2 = 1; 

Uo = [2.0 2.0 0.0 0.0]';
U1 = Uo(1);
U2 = Uo(2);

row1 = [-k1*U2         0         k_1           0  ];
row2 = [0           -k1*U1       k_1           0  ];
row3 = [k1*U2          0      -(k_1+k2)        k_2];
row4 = [0              0         k2           -k_2];

K = [row1;row2;row3;row4];

dt = 0.01;
time_vec = [0:dt:1];
nsteps = length(time_vec);
U = zeros(4,nsteps);
U(:,1) = Uo;

for i = 1:nsteps-1
    
    K = [row1;row2;row3;row4];
    
    P = eye(4) + K*dt/2;
    PU = P*U(:,i);
    % update K1 as:     K1 = K(P*U(:,i));
    row1 = [-k1*PU(2)         0         k_1           0  ];
    row2 = [0           -k1*PU(1)       k_1           0  ];
    row3 = [k1*PU(2)          0      -(k_1+k2)        k_2]; 
    K1 = [row1;row2;row3;row4];
    
    Q = eye(4) + K1*P*dt/2;
    QU = Q*U(:,i);
    % update K2 as:     K2 = K(Q*U(:,i)); 
    row1 = [-k1*QU(2)         0         k_1           0  ];
    row2 = [0           -k1*QU(1)       k_1           0  ];
    row3 = [k1*QU(2)          0      -(k_1+k2)        k_2]; 
    K2 = [row1;row2;row3;row4];

    R = eye(4) + K2*Q*dt;
    RU = R*U(:,i);
    % update K3 as:     K3 = K(R*U(:,i));
    row1 = [-k1*RU(2)         0         k_1           0  ];
    row2 = [0           -k1*RU(1)       k_1           0  ];
    row3 = [k1*RU(2)          0      -(k_1+k2)        k_2]; 
    K3 = [row1;row2;row3;row4];

    D = eye(4) + (dt/6)*(K + 2*K1*P + 2*K2*Q + K3*R);

    U(:,i+1) = D*U(:,i);

    % update K as:      K = K(D*U(:,i));    
    U1 = U(1,i+1);
    U2 = U(2,i+1);
    row1 = [-k1*U2         0         k_1           0  ];
    row2 = [0           -k1*U1       k_1           0  ];
    row3 = [k1*U2          0      -(k_1+k2)        k_2]; 
end

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

%% Solution of non-autonomous kinetic systems with MATLAB ODE solvers
tspan = [0 1];
[time_vec,u_vec] = ode23s(@Fp,tspan,Uo);
figure;plot(time_vec,u_vec,'-x');legend('U1','U2','U3','U4','Location','best');
title('ODE23s');grid on
[time_vec,u_vec] = ode15s(@Fp,tspan,Uo);
figure;plot(time_vec,u_vec,'-x');legend('U1','U2','U3','U4','Location','best');
title('ODE15s');grid on

%% An enzymatic reaction
tic

k1 = 10; 
k_1 = 1; 
k2 = 100; 
k_2 = 10; 

Uo = [2.0 0.2 0.0 0.0]';
U1 = Uo(1);
U2 = Uo(2);

row1 = [ -k1*U2         0         k_1            0    ];
row2 = [ 0           -k1*U1     (k_1+k2)      -k_2*U2 ];
row3 = [ k1*U2          0      -(k_1+k2)       k_2*U2 ];
row4 = [ 0              0         k2          -k_2*U2 ];

K = [row1;row2;row3;row4];

end_time = 3.0;
time_step = 0.01;
time_vec = [0:time_step:3];
nsteps = length(time_vec);
U = zeros(4,nsteps);
U(:,1) = Uo;

for i = 2:nsteps
    K = [row1;row2;row3;row4];
    [S,D] = eig(K);
    D = diag(D);
    Dt = D*time_step;
    eDt = diag(exp(Dt));
    U(:,i) = (S*eDt/S)*U(:,i-1);
    U1 = U(1,i);
    U2 = U(2,i);
    row1 = [ -k1*U2         0         k_1            0    ];
    row2 = [ 0           -k1*U1     (k_1+k2)      -k_2*U2 ];
    row3 = [ k1*U2          0      -(k_1+k2)       k_2*U2 ];
    row4 = [ 0              0         k2          -k_2*U2 ];
end

toc

figure;plot(time_vec,U,'-x')
legend('S','E','ES','P','Location','best')
title('Eigen solution');grid on

%% ODE15s and ODE23s
tspan = [0 3];
[time_vec,u_vec] = ode23s(@Fp_enz,tspan,Uo);
figure;plot(time_vec,u_vec,'-x');legend('S','E','ES','P','Location','best');
title('ODE23s');grid on
tic
[time_vec,u_vec] = ode15s(@Fp_enz,tspan,Uo);
toc
figure;plot(time_vec,u_vec,'-x');legend('S','E','ES','P','Location','best');
title('ODE15s');grid on

%% Chapter 13: Special Topics

%% Conversione of higher order differential equations into 1st order systems
 
%% Dumped spring
clear, clc
close all

m = 3;
b = 1; 
k = 3;

u_vec_0 = [0.1;1]
K = [-b/m -k/m;...
      1    0];

dt = 0.05;      % stepsize
end_time = 30;  % last time point in seconds
time_vec = [0:dt:end_time];
nsteps = length(time_vec);

%% Using RK4
I = eye(2);
K1 = K;
P = I + K*dt/2;
K2 = K*P;
Q = I + K*P*dt/2;
K3 = K*Q;
R = I + K*Q*dt;
K4 = K*R;
G = I + dt*(K1 + 2*K2 + 2*K3 + K4)/6;

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

udot_vec = (K*u_vec);
u_vec2 = [udot_vec(1,:) ; u_vec];

Dumped_oscillation_RK4 = figure;
plot(time_vec,u_vec2,'.-')
legend('Acceleration','Velocity','Displacement','Location','best')
title('RK4');xlim([0 end_time]);ylim([-1.2 1.2]);grid on

%% Using ODE15s
Fp_spring = @(t,u) K*u;
tspan = [0 end_time];

[time_vec,u_vec] = ode15s(Fp_spring,tspan,u_vec_0);
udot_vec = (K*u_vec')';
u_vec2 = [udot_vec(:,1) u_vec];
Dumped_oscillation_ODE15s = figure;
plot(time_vec,u_vec2,'.-');legend('Acceleration','Velocity','Displacement','Location','best');
title('ODE15s');
grid on;xlim([0 end_time]);ylim([-1.2 1.2]);

%% Leapfrog and velocity Verlet methods

%% A simple oscillator
close all
clear

x0 = 1;
v0 = 0;
xv0 = [x0; v0];
h = 2*pi/32
nsteps = 33
Circle_step = [0:pi/16:2*pi]
Circle_x = cos(Circle_step);
Circle_y = sin(Circle_step);

%% Forward Euler
xv = zeros(2,32);
xv(:,1) = xv0;
G = [1 h;-h 1];

for i = 2:nsteps
    xv(:,i) = G*xv(:,i-1);
end
xv = xv';    
figure;plot(Circle_x,Circle_y,':b')
hold on
plot(xv(:,1),xv(:,2),'-or')
axis equal
xlim([-2 2]);ylim([-1.4 1.8]);
% axescenter
hold off
vline([0],'y')
hline([0],'y')
xlabel('x');ylabel('v');
title('Forward Euler')

%% Backward Euler
xv = zeros(2,32);
xv(:,1) = xv0;
G = inv([1 -h;h 1]);

for i = 2:nsteps
    xv(:,i) = G*xv(:,i-1);
end
xv = xv';    
figure;plot(Circle_x,Circle_y,':b')
hold on
plot(xv(:,1),xv(:,2),'-or')
axis equal
xlim([-1.2 1.2]);ylim([-1.2 1.2]);
hold off
vline([0],'y')
hline([0],'y')
xlabel('x');ylabel('v');
title('Backward Euler')

%% Trapezoidal
xv = zeros(2,32);
xv(:,1) = xv0;
G = [1 -h/2;h/2 1]\[1 h/2;-h/2 1];

for i = 2:nsteps
    xv(:,i) = G*xv(:,i-1);
end
xv = xv';    
figure;plot(Circle_x,Circle_y,':b')
hold on
plot(xv(:,1),xv(:,2),'-or')
axis equal
xlim([-1.2 1.2]);ylim([-1.2 1.2]);
hold off
vline([0],'y')
hline([0],'y')
xlabel('x');ylabel('v');
title('Trapezoidal')
E = 0.5*(xv(:,1).^2+xv(:,2).^2);
figure; plot([0:32],E)
xlabel('step');ylabel('E');
title('Energy')

%% Velocity Verlet
xv = zeros(2,32);
xv(:,1) = xv0;
h2 = h*h;
h3 = h*h*h;
G = [1-0.5*h2 h;-h+0.25*h3 1-0.5*h2];

for i = 2:nsteps
    
% Original formulation    
%     xv(1,i) = xv(1,i-1) + h*xv(2,i-1) - 0.5*h2*xv(1,i-1);
%     xv(2,i) = xv(2,i-1) - 0.5*h*xv(1,i-1) - 0.5*h*xv(1,i);        

% Rearranged formulation
%     xv(1,i) = (1 - 0.5*h2)*xv(1,i-1) + h*xv(2,i-1) ;
%     xv(2,i) = (-h + 0.25*h3)*xv(1,i-1) + (1 - 0.5*h2)*xv(2,i-1) ;        

% Growth matrix
xv(:,i) = G*xv(:,i-1);
        
end
xv = xv';    
figure;plot(Circle_x,Circle_y,':b')
hold on
plot(xv(:,1),xv(:,2),'-or')
axis equal
xlim([-1.2 1.2]);ylim([-1.2 1.2]);
hold off
vline([0],'y')
hline([0],'y')
xlabel('x');ylabel('v');
title('Velocity Verlet')
E = 0.5*(xv(:,1).^2+xv(:,2).^2);
figure; plot([0:32],E)
xlabel('step');ylabel('E');
title('Energy')

%% Symbolic representation of the Stoichiometric, Gradient, and Jacobian matrices

syms x1 x2 x3 x4 k_1 k_m1 k_2 k_m2 t real
x = [x1 x2 x3 x4]'

S = [-1  1  0  0;
     -1  1  0  0;
      1 -1 -1  1;
      0  0  1 -1]

v1 = k_1*x1*x2
v2 = k_m1*x3
v3 = k_2*x3
v4 = k_m2*x4
v = [v1 v2 v3 v4]'

dxdt = S*v

G = [(1/2)*diff(v1,x1)  (1/2)*diff(v1,x2)  diff(v1,x3) diff(v1,x4);
           diff(v2,x1)        diff(v2,x2)  diff(v2,x3) diff(v2,x4);
           diff(v3,x1)        diff(v3,x1)  diff(v3,x3) diff(v3,x4);
           diff(v4,x1)        diff(v4,x1)  diff(v4,x3) diff(v4,x4)]

J = S*G 
dxdt = J*x

% Time scale decomposition of a small network of reactions
close all

k_1 = 1; 
k_m1 = .6; 
k_2 = .2; 
k_m2 = .02; 
 
Xo = [2.0 2.0 2.0 2.0]';
X1 = Xo(1);
X2 = Xo(2);

row1 = [-k_1*X2         0            k_m1             0  ];
row2 = [0             -k_1*X1        k_m1             0  ];
row3 = [(1/2)*k_1*X2  (1/2)*k_1*X1 -(k_m1+k_2)       k_m2];
row4 = [0              0             k_2            -k_m2];

J = [row1;row2;row3;row4];

time_step = 0.1;
time_vec = [0:time_step:20];
nstep = length(time_vec);

X = zeros(4,nstep);
W1 = zeros(4,nstep);
W2 = zeros(4,nstep);
W3 = zeros(4,nstep);
W4 = zeros(4,nstep);
X(:,1) = Xo;

    [V,D] = eig(J)
    tc = inv(D(1:3,1:3))
    w = V\X(:,1);    

    W1(:,1) = w(1)*V(:,1);
    W2(:,1) = w(2)*V(:,2);
    W3(:,1) = w(3)*V(:,3);
    W4(:,1) = w(4)*V(:,4);

for i = 2:nstep
    J = [row1;row2;row3;row4];
    [V,D] = eig(J);
    
    D = diag(D);    
    Dt = D*time_step;
    eDt = diag(exp(Dt));

    % Direct eigen solution
    % X(:,i) = (V*eDt/V)*X(:,i-1);
    
    % Separation of normal modes
    w = V\X(:,i-1);    

    W1(:,i) = w(1)*V(:,1)*eDt(1,1);
    W2(:,i) = w(2)*V(:,2)*eDt(2,2);
    W3(:,i) = w(3)*V(:,3)*eDt(3,3);
    W4(:,i) = w(4)*V(:,4)*eDt(4,4);
    X(:,i) = W1(:,i) + W2(:,i) + W3(:,i) + W4(:,i);
    
    X1 = X(1,i);
    X2 = X(2,i);
    
row1 = [-k_1*X2         0            k_m1             0  ];
row2 = [0             -k_1*X1        k_m1             0  ];
row3 = [(1/2)*k_1*X2  (1/2)*k_1*X1 -(k_m1+k_2)       k_m2];

end


figure;plot(time_vec,X,'-')
legend('X1','X2','X3','X4','Location','best')
title('Eigen solution');grid on
xlabel('Time (s)')
ylabel('[Compound]')

figure;plot(time_vec(1:end),W1 ,'-')
legend('U1','U2','U3','U4','Location','best')
title('Normal Mode 1 ');grid on
xlabel('Time (s)')
ylabel('[Compound component]')

figure;plot(time_vec(1:end),W2 ,'-')
legend('U1','U2','U3','U4','Location','best')
title('Normal Mode 2 ');grid on
xlabel('Time (s)')
ylabel('[Compound component]')

figure;plot(time_vec(1:end),W3 ,'-')
legend('U1','U2','U3','U4','Location','best')
title('Normal Mode 3 ');grid on
xlabel('Time (s)')
ylabel('[Compound component]')

figure;plot(time_vec(1:end),W4 ,'-')
legend('U1','U2','U3','U4','Location','best')
title('Normal Mode 4 ');grid on
xlabel('Time (s)')
ylabel('[Compound component]')

%% Phase portraits
close all

figure;plot(X(1,:),X(3,:),'-',X(1,:),X(4,:),'-',X(3,:),X(4,:),'-')
legend('X1 versus X3','X1 versus X4','X3 versus X4','Location','best')
title('Eigen solution ');grid on
xlabel('[Species 1]')
ylabel('[Species 2]')

figure;plot(W1(1,:),W1(3,:),'-',W1(1,:),W1(4,:),'-',W1(3,:),W1(4,:),'-')
legend('U1 versus U3','U1 versus U4','U3 versus U4','Location','best')
title('Normal Mode 1 ');grid on
xlabel('[Species 1]')
ylabel('[Species 2]')

figure;plot(W3(1,:),W3(3,:),'-',W3(1,:),W3(4,:),'-',W3(3,:),W3(4,:),'-')
legend('U1 versus U3','U1 versus U4','U3 versus U4','Location','best')
title('Normal Mode 3 ');grid on
xlabel('[Species 1]')
ylabel('[Species 2]')

figure;plot(W4(1,:),W4(3,:),'-',W4(1,:),W4(4,:),'-',W4(3,:),W4(4,:),'-')
legend('U1 versus U3','U1 versus U4','U3 versus U4','Location','best')
title('Normal Mode 4 ');grid on
xlabel('[Species 1]')
ylabel('[Species 2]')

