myfile = "kinetic_datafit_regression.sce";//determine the optimal values of the reaction rates k1 to k6 using the nonlinear datafit function

function [cmatrix] = kinetic_predicted_ode(c0k, k1k6k, xtime2)//predicted values using the ode integration function
    //c0k = initial molal concentration of the various constituents arranged as a column vector
    //t0k = initial time (seconds)
    //ttk = time nodes at which the intrations are required (seconds)
    //k1k6k = raction rate constants ()
    //xzero vector of zeros
    //disp("Pause at line 9 of function kinetic_predicted_ode");pause    
    //for j = 1:ynum
    t0k = xtime2(1);//use first node as initial time
    ntime = size(xtime2, "*");//number of time nodes
    ttk = xtime2(2:ntime);//times at which the integration shall be performed
    //disp("Pause at line 12 of function kinetic_predicted_ode.sci");pause
        cmatrix = ode(c0k, t0k, ttk, list(kinetic_cdot,  k1k6k));//integrate the function   for the current constituen    
endfunction

function [ee] = kinetic_gap_ode(k1k6a, c0k, xtime1, cdata)//Gap function to compute the difference between the measured and predicted molal concentration for use with the ode function
    //k1k6a = reaction rate constants
    //C01 = measures initial molal concentration values as a single column vector
    CInit =  c0k(:, 1);//vector of the molal concentrations at start of the reaction
    //xmatrixInit = matrix of molal concentration products obtained using the molal concentrations at starting time only
    //xtime1 = time nodes at which the solutions are to be determined
    //t01 = starting time of the reactions
        //disp("Pause at line 2 of function kinetic_gap");pause
        //xmatrix2 = values of molal concentration products from measured data
        //ym = derivatives of molal concentration from measured data       
    [cmatrix] = kinetic_predicted_ode(CInit, k1k6a, xtime1);//predicted values using the ode integration function
    ee = cdata - cmatrix;//difference in the given and calculated values of the molal concentrations
endfunction


function[xmatrix, yvector] = kinetic_xmatrix(ctgx, cdgx, cmgx, cglx, cex, cax, xtimex)//compute the matrix of products of molal concentrations
    //ctgx = concentration of tg species as single column vector
    //cdgx = concentration of dg species as single column vector
    //cmgx = concentration of mg species as single column vector
    //cglx = concentration of gl species as single column vector
    //cex = concentration of e (ester) species as single column vector
    //cax = concentration of a (alcohol) species as single column vector
    //xzero = zeros with dimensions corresponding with those of ctg and others as single column vector
    //yvector = vector of rates of change with time as single column vector
    //xtimex = vector of times at which experimental data are measured
    [mrow, mcol] = size(ctgx);//determine the number of rows
    xzerox = zeros(mrow, 1);//set dimensions to correspond to those of ctgx and others
x1 = [-(ctgx .* cax) (cdgx .* cex) xzerox  xzerox  xzerox xzerox];//row 1 of X
x2 = [(ctgx .* cax)  -(cdgx .* cex)  -(cdgx .* cax)  (cmgx .* cex)  xzerox  xzerox];//row 2 of x
x3 = [xzerox  xzerox  (cdgx .* cax)  -(cmgx .* cex)  -(cmgx .* cax)   (cglx .* cex)];//row 3 of x
x4 = [ xzerox  xzerox  xzerox  xzerox  (cmgx .* cax)  -(cglx .* cex)];//row 4 of x
x5 = [(ctgx .* cax)  -(cdgx .* cex)  (cdgx .* cax)  -(cmgx .* cex)  (cmgx .* cax) -(cglx .* cex)];//row 5 of x
x6 =  - x5;//row 6 of x
    xmatrix = [x1; x2;  x3;  x4;  x5;  x6];//arrange as a matrix 
    if (argn(2) >= 7) then
        y1 = splin(xtimex, ctgx);//rate of change of ctgx with time
        y2 = splin(xtimex, cdgx);//rate of change of cdgx with time        
        y3 = splin(xtimex, cmgx);//rate of change of cmgx with time        
        y4 = splin(xtimex, cglx);//rate of change of cglx with time        
        y5 = splin(xtimex, cex);//rate of change of cex with time        
        y6 = splin(xtimex, cax);//rate of change of cax with time        
        yvector = [y1; y2; y3; y4; y5; y6];//complete vector of rates of change with time
    end
    //disp("Pause at line 28 of function kinetic_xmatrix");pause
endfunction

function [cdot] = kinetic_cdot(t, cmatrix, k1k6d)//compute dy/dt
    //t = time (seconds)
    //cmatrix = matrix of molal concentration for all the constituents at all time nodes
    [krow, kcol] = size(k1k6d);//number of rows nd columns
    if (krow > 1) then
        k1k6b = k1k6d';//transpose
    else
        k1k6b = k1k6d;//leave unchanged
    end
    //[trow, tcol] = size(cmatrix);//number of rows and columns at current integration step
   //disp("Pause at line 13 of function kinetic_cdot.sci");pause
            xmatrix1 = kinetic_xmatrix(cmatrix(1),cmatrix(2), cmatrix(3), cmatrix(4), cmatrix(5), cmatrix(6));//compute the matrix of coefficients for the current time node t
    for i = 1:6//consider all the constituents   
        xm = [];//initialize
        for k = 1:6
            xm(k) = xmatrix1(i, k) * k1k6d(k);
        end
        cdot(i) = sum(xm);//predicted rate of reaction for constituent i at time node  j   
    end
    //disp("Pause at line 15 of function kinetic_cdot");pause// 
    //disp("urrent integration step ="); disp(tcol);
endfunction


xtime = [0	2	5	20	40	60	80]  * 60;//time (seconds)
tg =  [0.0869565	0.0492455	0.0277857	0.0087255	0.0002543	0.0002543	0.0001265];//mol/litre
dg= [0.0056604	0.0020696	0.0017374	0.0004606	0.0018868	0.0017374	0.0016718];//(mol/litre)
mg = 	[0.006383	0.0024572	0.0007586	0.0008941	0.001753	0.001358	0.001281];//(mol/litre)
g = [0.0003365	0.0008764	0.0012316	0.0011395	0.0056585	0.0060582	0.0084919];//(mol/litre)
e = [0.0004409	0.0433677	0.0507625	0.084096	0.0899657	0.0923184	0.0923184];//(mol/litre)
a = [0.0164359	0.0108777	0.0074697	0.0022768	0.0005042	0.0005042	0.0022768];//(mol/litre)
//p0 = [0.167 0.126 1.461 0.022 1.803 0.017]';//use values nitially obtained by Eze
b0 = ones(6, 1);//initial values
blow = zeros(6, 1);// + 0.001;// lower bounds
bhigh = ones(6, 1) * 2;//upper bound
xzero = zeros(7, 1);//column vector of 7 rows and 1 column filled with zeros
    cdata = [tg; dg; mg; g; e; a];//matrix of molal concentrations of each constrituent for all the measured times
    CInit1 = cdata(:, 1);//vector of initial molal concentrations
    iprint = 2;//print results at every iteration
    disp("Pause at line 103 of main script");pause
    [dmin1, popt, status] = datafit(iprint, list(kinetic_gap_ode, xtime, cdata), CInit1, "b", blow, bhigh, b0);//do the non-linear least squares regression data analysis for the current temperature of transesterification in degree Celcius

