% Fully developed Nusselt number for rectangular ducts
function [Nu, fReD, Dh, Ph] = RectNuInf(Wall,H,W)
  r=(H(:)./W(:))';
  r(r>1)=1./r(r>1);
  Dh=2*W.*H./(W+H);                                                     
  Ph=2*(W+H);
  Bn=pi*([1:500]-0.5)';
  F = (1-(6.*r).*sum(tanh(Bn./r)./Bn.^5));
  fReD = 96./F./(1+r).^2;                                               
  fReD(r<.001)=96;
  switch Wall
  case 'T'
    Nu=7.541-19.442*r+35.786*r.^2-34.597*r.^3+16.560*r.^4-2.8659*r.^5;  
  case 'H1'
    Nu=8.235-16.819*r+25.327*r.^2-20.079*r.^3+8.3064*r.^4-1.3622*r.^5;  
  case 'H2'
    Nu=2.9235-0.3737*r+2.4553*r.^2-3.5494*r.^3+2.0489*r.^4-0.4135*r.^5; 
  otherwise
    error('Wall= ''T'', ''H1'' or ''H2''');
  end
end
