% Fully developed Nusselt number for hexagonal tubes. 
function [Nu, fReD, Dh, Ph] = HexagNuInf(Wall,H,WW,SS)
  W=max(WW,SS);
  S=min(WW,SS);
  phi=iff(W>S,atan(H./(W-S)),pi/2);
  if min(phi) < 0.5 warning = 'HexagNuInf: phi < pi/6'
  end
  r=H./W;
  Dh=H.*(W+S)./(S+H.*csc(phi));                           
  Ph=2*(S+H./sin(phi));

  phi_Table=pi*[1/6 1/4 1/3 5/12 1/2];
  fRe_Table=r(:).^[0:7]*[96      96      96       96;     % Table 11-3
                    -142.29 -101.68 -103.38  -107.52;
                    -2.4124 -45.989   52.96   109.86;
                     685.24  662.32  169.15  -49.232;
                    -1410.6 -1673.7 -360.29   10.336;
                     1087.6  2268.1  324.07 -0.85724;
                     547.94 -1569.3 -143.19        0;
                    -1079.6  421.02  24.758        0];
  switch Wall
  case 'T'
    Nu_Table=r(:).^[0:8]*[7.541   7.541   7.541    7.541; % Table 11-3
                        -16.981 -13.089 -11.131  -14.012;
                        -35.108 -16.657 -28.576   17.007;
                         367.82  135.66  196.47  -9.7425;
                        -1112.8 -251.14 -434.82   2.8672;
                         1641.4  204.26  507.47 -0.42540;
                        -966.56 -63.553  -331.8 0.025251;
                              0       0  114.53        0;
                              0       0 -16.238        0];
  case 'H1'
    Nu_Table=r(:).^[0:8]*[8.235   8.235   8.235    8.235; % Table 11-3
                        -10.718 -15.758 -10.143  -12.927;
                        -118.93  25.463 -19.369   13.714;
                         939.99 -45.111  126.62  -6.3841;
                        -3139.4  75.780 -262.88   1.3521;
                         5031.3 -62.904  295.86 -0.10905;
                        -3079.8  17.932 -190.22        0;
                              0       0  65.246        0;
                              0       0 -9.2412        0];
  otherwise
    error('Wall= ''T'' or ''H1''');
  end

  rMax=tan(pi*[1/6 1/4 1/3 5/12]);             % lower Rhomb limit
  [NuRhomb fReRhomb]= RhombNuInf(Wall,rMax,1); 
  fRe_Table=iff(r(:)>rMax,fReRhomb,fRe_Table);
  Nu_Table=iff(r(:)>rMax,NuRhomb,Nu_Table);

  [NuRect fReRect]=RectNuInf(Wall,1,r);        
  fRe_Table=[fRe_Table, fReRect(:)];      
  Nu_Table=[Nu_Table, NuRect(:)];

  for i=1:length(r);                           % interpolate
    fReD(i)=interp1(phi_Table,fRe_Table(i,:),phi(i),'linear','extrap');
    Nu(i)=interp1(phi_Table,Nu_Table(i,:),phi(i),'linear','extrap');
  end
end
