clear all

nu_i=15.27e-6; k_i=25.74e-3; alpha_i=21.58e-6; pr_i=0.71; beta_i=0.00341;  % for T_f_i=T_inf_i
nu_o=11.71e-6; k_o=22.54e-3; alpha_o=16.3e-6; pr_o=0.719; beta_o=0.00395;  % for T_f_o=T_inf_o
nu=13.49e-6; k=0.0241; alpha=18.9e-6; pr=0.714; beta=0.00367;       % for the enclosure T_f=(T_inf_i+T_inf_o)/2      
H=1; L=0.025; g=9.8;

A_i=sqrt(H*k_i)*0.825; B_i=sqrt(H*k_i)*0.387/(1+(0.492/pr_i)^(9/16))^(8/27)*(g*beta_i*H^3/alpha_i/nu_i)^(1/6);
A_o=sqrt(H*k_o)*0.825; B_o=sqrt(H*k_o)*0.387/(1+(0.492/pr_o)^(9/16))^(8/27)*(g*beta_o*H^3/alpha_o/nu_o)^(1/6);
A=H^2*k/L*0.42*(g*beta*L^3/alpha/nu)^(1/4)*pr^0.012*(H/L)^(-0.3);

% Rearranging three simultaneous eqns. gives one equation about
% x=deltaT_i=T_inf_i-T_s_i
% Evaluate the left hand side f1, and right hand side f2 of the equation
x=[0:1:15]
siz_x=size(x);
for i=1:siz_x(2)
    f1(i)=A*(A_i+B_i*x(i)^(1/6))^2*x(i);
    f2(i)=(A_o+B_o*(40-x(i)-(A_i/sqrt(A)+B_i/sqrt(A)*x(i)^(1/6))^(8/5)*x(i)^(4/5))^(1/6))^2*...
        (40-x(i)-(A_i/sqrt(A)+B_i/sqrt(A)*x(i)^(1/6))^(8/5)*x(i)^(4/5));
end

% plot f1, f2 to estimate the location of the root for the equation
% which gives a range [10 14] for the fzero operation below
figure(1); clf
plot(x,f1,x,f2)
grid on
xlabel('deltaT_i')
legend('LHS','RHS')

% solve for deltaT_i = T_inf_i-T_s_i
options=[];
[dealT_i,fval]=fzero('doublepane',[10 14],options,A,A_i,B_i,A_o,B_o)

