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

  phi_Table=pi*[1/6 1/4 1/3 5/12 1/2];
  fRe_Table=r(:).^[0:6]*[96      96      96      96;     % Table 11-2
                    -265.34 -160.06 -155.64 -134.63;
                      58.32 -508.88  192.22  157.28;
                     2161.7  5793.2 -177.11 -74.336;
                    -6475.6  -22093  288.38  15.096;
                       9914   40324 -204.11 -1.8024;
                          0  -28142       0       0];
  switch Wall
  case 'T'
    Nu_Table=r(:).^[0:7]*[7.541   7.541   7.541   7.541; % Table 11-2
                        -39.551 -29.251 -23.175 -20.183;
                         155.04  90.601  40.937  33.483;
                        -1699.6 -535.56 -34.727 -26.286;
                         9621.5  2824.4  18.793  9.9036;
                         -16592 -8045.9 -7.6346 -1.4825;
                              0   11583       0       0;
                              0 -6702.8       0       0];
  case 'H1'
    Nu_Table=r(:).^[0:7]*[8.235   8.235   8.235   8.235; % Table 11-2
                        -36.138 -28.188 -20.951 -17.642;
                         107.54  134.47  39.219  24.711;
                        -1196.2 -1206.2 -82.077 -16.852;
                         6468.3  6854.2  172.33  7.0934;
                        -7824.5  -20679  -183.9 -2.2892;
                        -8280.3   31782  70.105 0.39886;
                              0  -19592       0       0];
  otherwise
    error('Wall= ''T'' or ''H1''');
  end

  rMax=tan(pi*[1/6 1/4 1/3 5/12])/2;         % lower Tri limit
  [NuTri fReTri]=IsoTriNuInf(Wall,rMax,1);   
  fRe_Table=iff(r(:)>rMax,fReTri,fRe_Table);
  Nu_Table=iff(r(:)>rMax,NuTri,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
