% Fully developed Nusselt number for circular-sector ducts. 
function [Nu, fReD, Dh, Ph] = CircSectNuInf(Wall,theta,R)
  Dh=2*R.*theta./(theta+1);                                             
  Ph=2*R.*(theta+1);
  n=[1:100]';
  kn=pi/2./theta.*(2*n-1);
  fReD=8*(theta.^2./(1+theta)).^2./sum((kn.*(kn+2)).^-2);               
  switch Wall
  case 'H1'
    denom=sum((kn.^2+7*kn+11)./(kn+2).^3./(kn+6)./(kn.*(kn+4).*(kn+1)).^2);
    Nu=128*(theta.^4./fReD./(1+theta).^3).^2./denom;                    
  case 'H2'
    NuHa=@(tt)-0.662*tt+15.3*tt.^2+10.8*tt.^3-59.6*tt.^4+38.5*tt.^5;
    NuHb=@(tt)0.754+6.55*tt-6.93*tt.^2+3.40*tt.^3-0.804*tt.^4+0.0741*tt.^5;
    Nu=pick_fun(theta,theta<pi/8,NuHa,NuHb);                            
  case 'T'
    NuTa=@(tt)0.945+13.9*tt-85.0*tt.^2+350*tt.^3-723*tt.^4+570*tt.^5;
    NuTb=@(tt)1.70+3.34*tt-2.935*tt.^2+1.357*tt.^3-0.3183*tt.^4+0.0299*tt.^5;
    Nu=pick_fun(theta,theta<pi/8,NuTa,NuTb);                            
  otherwise
    error('Wall= ''T'', ''H1'' or ''H2''');
  end
end
