%%%% solve ODE model to get output time profile
% Modified version (by Peter Halling) of Hui's program 
% This version tries to omit calculation of sensitivities, and focus just
% on progress of state variables.
% Also calls ODE_ProchiralOrdered.m

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

% Input values not normally to be changed
%initrate = 0.005; %sets time scale, used to calculate EobySo
% But now has to be set within function DoOneRun

% Input values here may be changed

showplot = false;

D0 = 0.00333; % M units
Q = 0.00167; % M units
% These first are assumed to remain constant via cofactor recycling
S0 = 0.2; % M units
Keq = 5;
Evalue = 39;
% Next 2 values here are with concn units (based on M), then first order const
% with s-1, then rest are dimensionless
KMD = 0.001;
KMS = 0.001;
k4 = 1000;
SCRrat = 1;
k34rat = 1;
k3rat = 1;
k14rat = 1;
k12Rrat = 1;
k12Srat = 1;

%dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, k3rat, ...
%        k14rat, k12Rrat, k12Srat, true);   % For test

if k4 < 0 % never
for D0 = [0.005 0.01 0.02] % 0.001 gives EobySo too big
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
D0 = 0.01
for Q = [0.001 0.005 0.01 0.02]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
Q = 0.01
end

if k4<0
for D0 = [0.005 0.00667 0.008 0.009]
    Q = 0.01 - D0;
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
for D0 = [0.0025 0.00333 0.004 0.0045]
    Q = 0.005 - D0;
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
for D0 = [0.00125 0.00167 0.002 0.00225]
    Q = 0.0025 - D0;
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
D0 = 0.00333
Q = 0.00167


for S0 = [0.05 0.1 0.2 0.4 1] 
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
S0 = 0.2
%end

for Keq = [1 2.5 5 10 20]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
Keq = 5

end
if k4 < 0 % never
for Evalue = [9 19 39 100]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
Evalue = 19
end

if k4 < 0 % never
for KMD = [0.0005 0.001 0.002 0.01 0.1]
    for KMS = [0.001 0.01 0.1 1]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    end
end
KMD = 0.001
KMS = 0.1
%end

%if k12Srat < 0
for KMD = [0.0005 0.001 0.002 0.01] 
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
KMD = 0.001
%end
%if k34Srat < 0 % never
for KMS = [0.0005 0.001 0.002] 
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
KMS = 0.001

%if k34Srat < 0 % i.e. never, just skip this bit
for k4 = [500 1000 2000]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
k4 = 1000
%end

for SCRrat = [0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
SCRrat = 1

for k34rat = [0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
k34rat = 1
for k3rat = [0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
k3rat = 1
%end

%if k14rat < 0
for k14rat = [0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
k14rat = 1
%end
%if k34Srat < 0 % i.e. never
for k12Rrat = [0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
k12Rrat = 1
for k12Srat = [0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
end
k12Srat = 1

end

if k4<0
for Keq = [1 2 5 10 20]
    for S0 = [0.05 0.1 0.2 0.5 1 2]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    end
end
KMD = 0.001
KMS = 0.001
end

if k4 < 0
for bigrat = [0.0001 0.0002 0.0005 0.001 0.002 0.005 0.01]
    for KMS = [0.0005 0.001 0.002]
        for SCRrat = [1 2 5 10]
            for k34rat = [0.5 1 2]
                k3rat = bigrat*k34rat/(KMS*SCRrat);
                if (k3rat > 0.09)
                    if (k3rat < 11)
                        dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
                    end
                end
            end
        end
    end
end

end

if k4 < 0
for SCRrat = [0.1 0.2 0.5 1 2 5 10]
    for KMS = [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 (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    end
end
SCRrat = 1
KMS = 0.001

for k3rat = [0.1 0.2 0.5 1 2 5 10]
    for k34rat = [0.1 0.2 0.5 1 2 5 10]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    end
end
SCRrat = 1
KMS = 0.001

%end

%KMS = 0.1
for SCRrat = [0.1 0.2 0.5 1 2 5 10]
    for k3rat = [0.2 0.5 1 2 5]
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    end
end
SCRrat = 1
k3rat = 1

end
if k4 < 0

for S0 = [1] %[0.05 0.1 0.2 0.5 1]
    for Evalue = [19] % 39 99 199] 
    dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    end
end
S0 = 0.2
Evalue = 39

end

    dummy = DoOneRun (D0, Q, 0.5, Keq, 19, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);
    dummy = DoOneRun (D0, Q, 1, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ...
        k3rat, k14rat, k12Rrat, k12Srat, false);





function dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, ...
    k34rat, k3rat,  k14rat, k12Rrat, k12Srat, showplot)

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

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

k3R=k34rat*k4;
k3S=k3R/(k3rat*Evalue);
thisdenom = KMS*(2*Evalue*k3R*k3S + Evalue*k3R*k4 + Evalue*k3S*k4 + k3R*k3S + k3R*k4);
SCRf=k4*k3R*k3S*Evalue/thisdenom;
k1=SCRf*(Evalue+1)*KMS/(Evalue*KMD);
SCSf=SCRf/Evalue;
SCRb=SCRf/SCRrat;
SCSb=SCRb/Evalue;
km4=k1*k14rat;
km1 = k1*k4*SCRf/(Keq*km4*SCRb);
km2R=k12Rrat*km1;
km2S = k12Srat*km1/Evalue;
k2R = SCRf*(k3R + km2R)/k3R;
km3R = SCRb*(k3R + km2R)/km2R;
k2S = SCSf*(k3S + km2S)/k3S;
km3S = SCSb*(k3S + km2S)/km2S;

term1 = (1/k4 + 1/k3S)/Evalue + 1/k4 + 1/k3R + (Q*km4 + k4)*(1 + 1/Evalue)/(k1*D0*k4);
term2 = (1 + km1*(Q*km4/k4 + 1)/(k1*D0))/(S0*SCRf);
EobySo = initrate*(term1 + term2);

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

% Remove dimensions
k1n = S0*k1;
k2Rn = S0*k2R;
k2Sn = S0*k2S;
km3Rn = S0*km3R;
km3Sn = S0*km3S;
km4n = S0*km4;


Sn0 = 1; % By definition of non-dimensional basis
En0 = 1; % Same
Dn0 = D0/S0; 
Qn0 = Q/S0;
% All other state variables start at 0

% Set up the values of state variables and simulation parameters
% initial condition of state variables
Xinit = [En0 0 0 0 0 Dn0 Sn0 0 0 Qn0]';
% nominal parameter values
k = [k1n k2Rn k2Sn k3R k3S k4 km1 km2R km2S km3Rn km3Sn km4n]';

% set up the simulation time and sampling interval
%tspan = 0:0.001:0.1; % start with transient
%tspan = 0:1:100; % initial rate period
%tspan = 0:20:1000; % further progress
tspan = 0:100:5000; % to get better conversion on slow reactions
%tspan = 0:200:10000; % to get better conversion on slow reactions


if ~toobig

% 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;
%AbsTol_value=[1e-11;1e-11;1e-13;1e-11;1e-12;(1e-6).*ones(5,1);(1e-13).*ones(11,1);...
%    (1e-14).*ones(11,1);(1e-15).*ones(11,1);(1e-13).*ones(11,1);(1e-14).*ones(11,1);...
%    (1e-9).*ones(55,1)]; %??????????
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_ProchiralOrdered(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(8,:);
PSn = x(9,:);
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 = ["D0", "Q", "S0", "Keq", "Eval", "KMD", "KMS", "k4", "SCRrat", ...
    "k34rat", "k3rat", "k14rat", "k12Rrat", "k12Srat"];
values = [D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, ...
    k34rat, k3rat, k14rat, k12Rrat, k12Srat];
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

    %dim = [0.15 .6 .3 .3];
    %annotation('textbox',dim,'String',lbl_text,'FitBoxToText','on');
    %dim = [0.15 .53 .3 .3];
    %annotation('textbox',dim,'String',lbl_text2,'FitBoxToText','on');
    %dim = [0.15 .46 .3 .3];
    %annotation('textbox',dim,'String',lbl_text3,'FitBoxToText','on');
    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,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)
%rowno = startrow;
%foundit = false;
%while ~foundit
%    cellref = MakeCellRef(int8(rowno), int8(colno))
%    entry = readmatrix(filename, 'Sheet', 1, 'Range', strcat(cellref, ':', cellref))
    % Doesn't work, reads series of rows, then type mismatch on next
    % statement
%    if entry == ''
%        foundit = true
%    end
%    rowno = rowno + 1;
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