//NL-regresssion for kLa
clear;clc;format(12);
//Dimensions, physical data
Drx = 0.109; // reactor internal diameter in m
L = 0.500; //reactor length in m
d = 0.050; //impeller diameter in m
rho = 993.33; //density of water in g/cm3 at 37 deg C, 1 atm, CRC Handbook 2015
mu = 0.0006915; //viscosity of water in Pa.s at 37 deg C, 1 atm, CRC Handbook 2015
Nq=0.19; //Tsui & Hu 2008
NpRe=114; //Tsui & Hu 2008
k1=NpRe/Nq; //<-- this value is a constant for higher Re numbers**
di = 0.0032; //internal diameter of air tube in m
d_he = 0.022; //<-- outer diameter of heating element
L_he = 0.292; //<-- length of heating element

//The function and loop below calculate the impeller coverage
function VL=area(hL, D)
    r = D/2;
    if hL < r then
        t = 2.*acos((r-hL)./r);
        AL = ((r.^2).*(t-sin(t)))/2;
    else
        t = 2.*acos((r-(D-hL))./r);
        AL = (%pi.*(r.^2))-(((r.^2).*(t-sin(t)))/2);
    end
    VL = AL.*L;
endfunction

Vhe=(%pi.*(d_he.^2)/4)*L_he; //<-- heating element volume displaced
hL=linspace(0,0.109,500); //<-- all possible liquid heights
VL=area(hL,Drx); //<-- all posssible liquid volumes\
//The loop below extracts liquid heights at integer volumes 1-4 L
Ve=[0.001 0.002 0.003 0.004]; Vf=Ve+Vhe; Vx=[]; H=[]; tol=5e-6;
for i=1:4
    for n=1:length(VL)
        if abs(VL(n)-Vf(i))<tol then
            Vx=[Vx VL(n)];
            H=[H hL(n)];
        end
    end
end
//The vector H now stores the liquid height in m, at Volumes of 1- 4 L
hLd=H-((Drx-d)/2); // <-- height of water covering the impeller
hLd_D=hLd./d; // <-- Impeller immersion ratio

//Extract kLa data
//Data = readxls('kLa_analysis.xls'); D=Data(2);//<--kLa data
D=csvRead('C:\Users\SHRRAJ002\Desktop\Wesley_Andrew_data\Conductivity experiment\Data files\kLa_analysis.csv');
hLdD=D(2:41,2); Vs=D(2:41,3); N=D(2:41,4); kLa=D(2:41,5);
//Generate dimensionless groups
NRe=((rho.*N).*(d.^2))./mu;
Nq=D(2:41,6); //<-- circulation flow number, estimated from mixing time
Np=k1./(NRe./Nq); //<-- predicts power number
NVs=Vs./(N.*d);
M=[hLdD NRe Nq Np NVs];

//Call leastsq function
//Model form based on dimensional analysis
function k=kla(M, I)
    k=I(1).*((M(:,1)).^(I(2))).*(M(:,2).^(I(3))).*(M(:,3).^(I(4))).*(M(:,4).^(I(5))).*(M(:,5).^(I(6)));
endfunction

//Define functions to minimize
function y=func(I, M, kLa)
    y=kla(M, I)- kLa;
endfunction

I0=[0.5 0.1 0.1 0.1 0.1 0.1]'; M2=M(1:20,1:$); kLa2=kLa(1:20); M3=M(21:40,1:$); kLa3=kLa(21:40);
[K1 K2 K3] = leastsq(list(func,M,kLa),I0);
[K2_1 K2_2 K2_3] = leastsq(list(func,M2,kLa2),I0);
[K3_1 K3_2 K3_3] = leastsq(list(func,M3,kLa3),I0);

//Calculate the single variables exponents e.g. d, n, Q etc.
dd=(-1*K2(2))+(-5*K2(5))+(2*K2(3))+(-1*K2(6))+(-3*K2(4));
nn=(-3*K2(5))+(1*K2(3))+(-1*K2(6))+(-1*K2(4));
rr=(-1*K2(5))+(1*K2(3));
KK1=K2(1)*(mu^-K2(3))*(rho^rr);

Tm2=[72 41.333333 27.666667 20.666667 20.333333]'; //<-- experimentally recorded mixing times at 2L
Tm3=[70.666667 24.666667 19.666667 17.666667 16.333333]'; //<-- experimentally recorded mixing times at 3L
Tm=[Tm2;Tm2;Tm2;Tm2;Tm3;Tm3;Tm3;Tm3];

//Models to use for predicting kLa
//Single model 2L
function kLa=model2L(n, Vs);
    NRe=((rho.*n).*(d.^2))./mu;
    //Tm=109.6.*(n.^(-0.834)); //<-- model for mixing time
    Nq=(%pi./(n.*Tm)).*(L./Drx).*((Drx./d).^3);
    Np=600.*(Nq./NRe);
    NVs=Vs./(n.*d);
    kLa = K2_2(1).*((hLd_D(2)).^K2_2(2)).*(NRe.^K2_2(3)).*(Nq.^K2_2(4)).*(Np.^K2_2(5)).*(NVs.^K2_2(6));
endfunction

//Different models 2L
function kLa=model2L_comb(n, Vs);
    NRe=((rho.*n).*(d.^2))./mu;
    //Tm=109.6.*(n.^(-0.834)); //<-- model for mixing time
    Nq=(%pi./(n.*Tm)).*(L./Drx).*((Drx./d).^3);
    Np=600.*(Nq./NRe);
    NVs=Vs./(n.*d);
    kLa = K2(1).*((hLd_D(2)).^K2(2)).*(NRe.^K2(3)).*(Nq.^K2(4)).*(Np.^K2(5)).*(NVs.^K2(6));
endfunction

//Single model 3L
function kla=model3L(n, Vs);
    NRe=((rho.*n).*(d.^2))./mu;
    //Tm=93.77.*(n.^(-0.899)); //<-- model for mixing time
    Nq=(%pi./(n.*Tm)).*(L./Drx).*((Drx./d).^3);
    Np=600.*(Nq./NRe);
    NVs=Vs./(n.*d);
    kla = K3_2(1).*((hLd_D(2)).^K3_2(2)).*(NRe.^K3_2(3)).*(Nq.^K3_2(4)).*(Np.^K3_2(5)).*(NVs.^K3_2(6));
endfunction

//Different models 3L
function kla=model3L_comb(n, Vs);
    NRe=((rho.*n).*(d.^2))./mu;
    //Tm=93.77.*(n.^(-0.899)); //<-- model for mixing time
    Nq=(%pi./(n.*Tm)).*(L./Drx).*((Drx./d).^3);
    Np=600.*(Nq./NRe);
    NVs=Vs./(n.*d);
    kla = K2(1).*((hLd_D(2)).^K2(2)).*(NRe.^K2(3)).*(Nq.^K2(4)).*(Np.^K2(5)).*(NVs.^K2(6));
endfunction

//Calculate coefficient of determination
//single model
y_bar=(sum(kLa)/length(kLa)); SS_total=0; //<-- mean kLa
for c=1:length(kLa)
    SS_total=SS_total+((kLa(c)-y_bar).^2)
end
R2=1-(K1/SS_total);

//Different models
y_bar2=(sum(kLa2)/length(kLa2)); SS_2total=0; //<-- mean kLa
for c=1:length(kLa2)
    SS_2total=SS_2total+((kLa2(c)-y_bar2).^2)
end
R2_2=1-(K2_1/SS_2total);

y_bar3=(sum(kLa3)/length(kLa3)); SS_3total=0; //<-- mean kLa
for c=1:length(kLa3)
    SS_3total=SS_3total+((kLa3(c)-y_bar3).^2)
end
R2_3=1-(K3_1/SS_3total);

//plot 3d surfaces for kLa and experimental points
//For 2L:
scf(0);clf;
Nm=N(1:5); Vm=[Vs(1) Vs(6) Vs(11) Vs(16)];
//plot experimetnal results
kLa2m=[kLa2(1:5) kLa2(6:10) kLa2(11:15) kLa2(16:20)];
plot3d(Nm,Vm,kLa2m); //plot x=N , y=Vs , z= kla
title('2L kLa model vs. experimental','fontsize', 4);
xlabel('n [1/s]','fontsize', 3);
ylabel('Vs [m/s]','fontsize', 3);
zlabel('kLa [1/h]', 'fontsize', 3);
xgrid(0);
a=gca(); b=a.children(1);
b.surface_mode="off"; b.mark_mode="on"; b.mark_size=1.0; b.mark_style=0;
//plot model predictions
//kmodel2=model2L(N, Vs); //<-- enable to use separate models
kmodel2=model2L_comb(N, Vs); //<-- enable to use single model
Km2=[kmodel2(1:5) kmodel2(6:10) kmodel2(11:15) kmodel2(16:20)];
plot3d(Nm, Vm, Km2);
f=gcf();
f.color_map = graycolormap(128);
e=gca(); g=e.children(1);
g.color_flag=1;
//plot lines linking data points to surface
scf(0);
for m=1:length(Nm)
    for p=1:length(Vm)
        param3d([Nm(m) Nm(m)], [Vm(p) Vm(p)], [Km2(m,p) kLa2m(m,p)]);
    end
end

//For 3L:
scf(1);clf;
//plot experimetnal results
kLa3m=[kLa3(1:5) kLa3(6:10) kLa3(11:15) kLa3(16:20)];
plot3d(Nm,Vm,kLa3m); //plot x=N , y=Vs , z= kla
title('3L kLa model vs. experimental','fontsize', 4);
xlabel('n [1/s]','fontsize', 3);
ylabel('Vs [m/s]','fontsize', 3);
zlabel('kLa [1/h]', 'fontsize', 3);
xgrid(0);
a=gca(); b=a.children(1);
b.surface_mode="off"; b.mark_mode="on"; b.mark_size=1.0; b.mark_style=0;
//plot model predictions
//kmodel3=model3L(N, Vs); //<-- enable to use separate models
kmodel3=model3L_comb(N, Vs); //<-- enable to use single model
Km3=[kmodel3(21:25) kmodel2(26:30) kmodel2(31:35) kmodel2(36:40)];
plot3d(Nm, Vm, Km3);
f=gcf();
f.color_map = graycolormap(128);
e=gca(); g=e.children(1);
g.color_flag=1;
scf(1);
for m=1:length(Nm)
    for p=1:length(Vm)
        param3d([Nm(m) Nm(m)], [Vm(p) Vm(p)], [Km3(m,p) kLa3m(m,p)]);
    end
end

//Report results of least squares regression
//Single model
disp('Single model for kLa at all volumes is');
format(5);
disp('kLa = '+string(K2(1))+'.(hLd/d^'+string(K2(2))+').(NRe^'+string(K2(3))+').(Nq^'+string(K2(4))+').(Np^'+string(K2(5))+').(NVs^'+string(K2(6))+')');
disp('R^2 coefficient of determination for combined model is '+string(R2));
disp('');
//Different models
disp('Model for kLa at 2L volume is');
format(5);
disp('kLa = '+string(K2_2(1))+'.(hLd/d^'+string(K2_2(2))+').(NRe^'+string(K2_2(3))+').(Nq^'+string(K2_2(4))+').(Np^'+string(K2_2(5))+').(NVs^'+string(K2_2(6))+')');
disp('R^2 coefficient of determination for 2L model is '+string(R2_2));
disp('');
disp('Model for kLa at 3L volume is');
format(5);
disp('kLa = '+string(K3_2(1))+'.(hLd/d^'+string(K3_2(2))+').(NRe^'+string(K3_2(3))+').(Nq^'+string(K3_2(4))+').(Np^'+string(K3_2(5))+').(NVs^'+string(K3_2(6))+')');
disp('R^2 coefficient of determination for 3L model is '+string(R2_3));
format(12);

 
