% Fully developed Nusselt number for elliptical ducts. 
function [Nu, fReD, Dh, Ph] = EllipNuInf(Wall,H,W)
  eta=H./W;
  eta(eta>1)=1./eta(eta >1);
  H=min(H,W);
  s=(1-eta).^2./(1+eta).^2;
  E=pi/4*(eta+1).*(1+3*s./(10+sqrt(4-3*s)));             % Approx.
  E=nthargout(2,@ellipke,1-eta.^2);                      % Exact
  Dh=pi/2*H./E;                                          
  Ph=2.*W.*E;
  fReD=8*(1+eta.^2).*(pi./E).^2;                         

  switch Wall
  case 'T'
    Nu=3.5+2.759*eta-8.987*eta.^2+12.42*eta.^3+ ...
          -8.167*eta.^4+2.132*eta.^5;                    
  case 'H1'
    Nu=9*(pi./E).^2.*(eta.^6+7*eta.^4+7*eta.^2+1)./ ...
                     (17*eta.^4+98*eta.^2+17);           
  case 'H2'
    Nu=-0.9273+17.40*eta-19.65*eta.^2+7.834*eta.^3+ ...
          -0.2862*eta.^4;                                
  otherwise
    error('Wall= ''T'', ''H1'' or ''H2''');
  end
end
