% Natural convection between expansive horizontal plates
function [h, NuS, RaS] = HorzPlateEnclNuT(Tbot,Ttop,S,mat)
  [RaS,Pr,kf] = Ra(Tbot,Ttop,S,mat);                 
  if Ttop > Tbot
    NuS = 1;
  else
    a1=1.44/(1+0.018/Pr+0.00136/Pr^2);               
    Nhat=RaS^(1/3)/75/exp(1.5*Pr^-0.5);              
    Pos=@(x)(abs(x)+x)/2;                            
    NuS=2*(1+Pos(1-1708/RaS)*(a1+2*Nhat^(1-log(Nhat)))+ ...
           Pos((RaS/5830)^(1/3)-1));                 
  end
  h=NuS.*kf/S;
end
