%%%% solve ODE model to get output time profile
% Program by Peter Halling (Uni Strathclyde), based on an original program by Hui Yu
% Also calls ODE_ProchiralDiacid.m

%clear all
clear variables %following help guidance
close all
clc 
format long

% Input values here may be changed

showplot = false;

S0 = 0.2; % M units
B0 = 55;
Q0 = 0;
Keq = 0.25;  
Evalue = 39;
% Next value here has concn units (M), then first order const
% with s-1, then rest are dimensionless
KMS = 0.1; % M units
k2R = 1000; %s-1
SCRrat = 0.5;
Erat = 0.5;
SCR21rat = 0.02;
SCSrat = 1;
k24rat = 1;
k2rat = 1;
k24Srat = 1;
k12Rrat = 2;
k12Srat = 2;
k34Rrat = 2;
k34Srat = 2;


if k2R < 0  % i.e. never, used to skip bits not wanted in current run
dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, true);   % For test
end
    
  
if k2R < 0

for S0 = [0.05 0.1 0.2 0.4 1] 
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
S0 = 0.2
end
if k2R < 0

for Q0 = [0 0.25 0.5 1 2] 
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
Q0 = 0.5

end
if k2R < 0

for Keq = [0.125 0.25 0.5 1 2 5 10 20 50 100]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
Keq = 0.25

%end
%if k2R < 0 % never
for Evalue = [9 19 39 100]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
Evalue = 39
%end

%if k2R < 0
for KMS = [0.001 0.01 0.05 0.1 0.2 1 10] 
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
KMS = 0.1

end

for k2R = [500 1000 2000]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k2R = 1000
%end

%if k2R < 0

for SCRrat = [0.25 0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
SCRrat = 0.5

for Erat = [0.25 0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
Erat = 0.5

for SCR21rat = [0.01 0.02 0.04 0.1 0.2 0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
SCR21rat = 0.02

for SCSrat = [0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
SCSrat = 1

for k24rat = [0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k24rat = 1
for k2rat = [0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k2rat = 1


for k24Srat = [0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k24Srat = 1



%end

%if k2R < 0 
for k12Rrat = [1 2 4]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k12Rrat = 2
for k12Srat = [1 2 4]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k12Srat = 2

for k34Rrat = [1 2 4]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k34Rrat = 2
for k34Srat = [1 2 4]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
end
k34Srat = 2

%end

if k2R<0
for Keq = [1 2 5 10 20]
    for S0 = [0.05 0.1 0.2 0.5 1 2]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
    end
end
Keq = 1
S0 = 0.2
%end

for SCRrat = [0.1 0.2 0.5 1 2 5 10] %[0.1 0.2 0.5 1 2 5 10]
    for SCR21rat = [0.2 0.5 1 2 5 10] % 1 2 5 10] %[0.0005 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
    end
end
SCRrat = 0.5
SCR21rat = 1

for Erat = [0.1 0.2 0.5 1 2 5 10] %[0.1 0.2 0.5 1 2 5 10]
    for SCSrat = [0.1 0.2 0.5 1 2 5 10] % 1 2 5 10] %[0.0005 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5]
       dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, false);  
    end
end
Erat = 0.5
SCSrat = 1

end

% Now the function, which doesn't actually return anything useful, just
% does some jobs

function dummy = DoOneRun (S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat, showplot)

initrate = 0.005; %sets time scale, used to calculate EobySo

% Now calculate rate constants, at first all on dimensioned basis

k4R = k2R*k24rat;
k2S = k2R/(k2rat*Evalue);
k4S = k24Srat*k2S;

denom1 = B0*KMS*SCR21rat*(k2rat + 1 + k2rat/k24Srat + 1/k24rat);
SCRf = k2R*(B0*SCR21rat - Erat*KMS - KMS)/denom1;
if SCRf < 0
    error('SCRf negative')
end

SCRb = SCRf/SCRrat;
SCSf = SCRf/Evalue;
SCR2f = SCR21rat*SCRf;
SCS2f = SCR2f/(Evalue*Erat);
SCR2b = SCRf*SCR2f/(SCRb*Keq);
SCSb = SCSf/SCSrat;
SCS2b = SCRb*SCR2b*SCSf*SCS2f/(SCSb*SCRf*SCR2f);

km1R = k12Rrat*k2R;
k1R = SCRf*(1+k12Rrat);
km2R = SCRb*(1+1/k12Rrat);

km1S = k12Srat*k2S;
k1S = SCSf*(1 + k12Srat);
km2S = SCSb*(1 + 1/k12Srat);
km3R = k34Rrat*k4R;
k3R = SCR2f*(1 + k34Rrat);
km4R = SCR2b*(1 + 1/k34Rrat);
km3S = k34Srat*k4S;
k3S = SCS2f*(1 + k34Srat);
km4S = SCS2b*(1 + 1/k34Srat);


EobySo = initrate*(1/(S0*SCRf)+1/k4R+1/k2R+1/(k4S*Evalue)+1/(k2S*Evalue)+ (1+Erat)/(B0*SCR21rat*SCRf))

if EobySo > 0.01
    toobig = true;
    %error ('EobySo > 0.01')
else
    toobig = false;
end

% Remove dimensions
k1Rn = S0*k1R;
k1Sn = S0*k1S;
km2Rn = S0*km2R;
km2Sn = S0*km2S;
k3Rn = S0*k3R;
k3Sn = S0*k3S;
km4Rn = S0*km4R;
km4Sn = S0*km4S;


Sn0 = 1; % By definition of non-dimensional basis
En0 = 1; % Same
Bn0 = B0/S0;
Qn0 = Q0/S0;

% All other state variables start at 0

% Set up the values of state variables and simulation parameters
% initial condition of state variables

% En  ESRn  ESSn  EsRn  EsSn  EPRn  EPSn  Sn  Qn  Bn  PRn  PSn

Xinit = [En0 0 0 0 0 0 0 Sn0 Qn0 Bn0 0 0]';
% values of all rate constants used in integrating model
k = [k1Rn k1Sn k2R k2S k3Rn k3Sn k4R k4S km1R km1S km2Rn km2Sn km3R km3S km4Rn km4Sn]';

% set the simulation time and sampling interval
%tspan = 0:0.001:0.1; % start with transient
%tspan = 0:0.4:20; % initial rate period
tspan = 0:20:1000; % further progress
%tspan = 0:20:3000; % further progress

if ~toobig % otherwise skip this trial and wait for another

% With normalised non-dimensional values for all variables and parameters,
% probably similar AbsTol for all, so try keeping default 1e-6.
AbsTol_value = 1e-6;
options = odeset('AbsTol',AbsTol_value);
% Run the ODE simulation
%[t, x_s] = feval(@ode15s, @ODE_Prochiral, tspan, xs_init, options, k, EobySo); 
% This doesn't work any more when have extra parameter EobySo, but managed
% to avoid feval call completely after viewing ode15s help
%[t,y] = ode15s(@(t,y) odefcn(t,y,A,B), tspan, y0);  Example from help
[t, xout] = ode15s(@(t,x) ODE_ProchiralDiacid(t,x,k,EobySo), tspan, Xinit, options);

% Help gives "[t,y] = ode15s(odefun,tspan,y0). Each row in the solution array
%y corresponds to a value returned in column vector t."
%disp(xout)
%Yes, xout has a column for each state variable, and a row for each time, so
x = transpose(xout);

% x has one row for each state variable, one col for each time
PRn = x(11,:);
PSn = x(12,:);
diff = PRn - PSn;
conversion = PRn + PSn;
ees = diff./conversion;

else % is toobig
    showplot = false;
    t = '*';
    conversion = 'EobySo too big';
    ees = 'EobySo too big';
end

% Prepare for output
names = ["S0", "Q0", "B0", "Keq", "Eval", "KMS", "k2R", "SCRrat", "Erat", ...
    "SCR21rat", "SCSrat", "k24rat", "k2rat", "k24Srat", "k12Rrat", "k12Srat", ...
    "k34Rrat", "k34Srat"];
values = [S0, Q0, B0, Keq, Evalue, KMS, k2R, SCRrat, Erat, SCR21rat, ...
    SCSrat, k24rat, k2rat, k24Srat, k12Rrat, k12Srat, k34Rrat, k34Srat];
nonames = length(names);
if length(values) ~= nonames 
    error ('Different length names and values')
end



if showplot
    nolines = fix((nonames-1)/5) + 1;
    for lineno = 1:nolines
      firstctr = 5*(lineno-1) + 1;
      lastctr = 5*(lineno-1)+ 5;
      if lastctr > nonames
        lastctr = nonames;
      end
      lbl(lineno)="";
      for ctr = firstctr:lastctr
        lbl(lineno) = strcat(lbl(lineno), names(ctr), "=", num2str(values(ctr)), " ");
        %thislbl = strcat(lbl(lineno), names(ctr), "=", num2str(values(ctr)), " ")
      end
  end
 
    figure (1)
    %plot(t,x(1,:),'b--',t,x(2,:),'r-',t,x(3,:),'g-.',t,x(4,:),'m-',t,x(5,:),'c-.','LineWidth',2.5);
    %plot(t,x(6,:),'b--',t,x(7,:),'r-',t,x(8,:),'g-.',t,x(9,:),'m-',t,x(10,:),'c-.','LineWidth',2.5);
    plot(t,conversion,'b--',t,ees,'r-',t,x(1,:),'g-.','LineWidth',2.5);
    hold on;
    %plot(t(b),a,'g.','MarkerSize',24);
    %text(t(b),a+0.05,max_text);xlabel('\fontsize{16}time/s');

    % From help
    %dim = [.2 .5 .3 .3];
    %annotation('textbox',dim,'String',str,'FitBoxToText','on');
    % By trial and error, x and y run from origin in bottom left

    for ctr = 1:nolines
        dim = [0.15 0.67-0.07*ctr .3 .3];
        annotation('textbox',dim,'String',lbl(ctr), 'FitBoxToText', 'on');
    end

    ylabel('\fontsize{16}non-dimensional concentrations');
    %title('\fontsize{14}time profile of state variables');
    %legend('En','EDn','EtRn','EtSn','EQn');
    %legend('Dn','Qn','Sn','PRn','PSn');
    legend('Conv','ee_P','En');
end

% Output_progress

filename = 'C:\Users\ph\My Documents\A Edited\ProgressCurves.xlsx';
%rowno = int8(5)
startcol = 5;

firstdatacol = startcol + nonames + 2;
rowno = FindBlankRow(filename, 5, firstdatacol);

for ctr = 1:nonames
    writematrix(names(ctr), filename, 'Sheet', 1, 'Range', MakeCellRef(rowno, startcol+ctr-1));
    writematrix(values(ctr), filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+ 1, startcol+ctr-1));
end

writematrix('time',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno, firstdatacol));
writematrix('conv',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+1, firstdatacol));
writematrix('ee_P',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+2, firstdatacol));

%nextref = MakeCellRef(rowno,startcol+15);
nextref = MakeCellRef(rowno,firstdatacol+1);
writematrix(transpose(t),filename,'Sheet',1,'Range',nextref);
nextref = MakeCellRef(rowno+1, firstdatacol+1);
writematrix(conversion,filename,'Sheet',1,'Range',nextref);
nextref = MakeCellRef(rowno+2, firstdatacol+1);
writematrix(ees,filename,'Sheet',1,'Range',nextref);

dummy = 1; %to prevent warning

end


function rowno = FindBlankRow(filename, startrow, colno)

cellref = MakeCellRef(int8(startrow), int8(colno));
entries = readmatrix(filename, 'Sheet', 1, 'Range', strcat(cellref, ':', cellref));
norows = length(entries);
rowno = startrow + norows;
end

function cellref = MakeCellRef(rowno, colno)
firstletno = idivide(colno - 1, int8(26));
secondletno = mod(colno - 1, int8(26)) + 1;
if firstletno > 0 
    cellref = char(firstletno + 64); 
else
    cellref = '';
end
cellref = strcat(cellref,char(secondletno+64));
cellref = strcat(cellref,num2str(rowno));
end