% Function to mix laminar and turbulent convection correlations 
function Nu = MixNu(flag,NuLam,NuTurb,Rc,R,Pr,n)
  Nu=zeros(size(R));
  j=find(R < Rc);            % find laminar cases
  if ~isempty(j)
    Prj=Pick(Pr,j);                                           
    Nu(j) = NuLam(R(j),Prj);
  end
  j=find(R >= Rc);           % find mixed turbulent cases
  if ~isempty(j)
    Prj=Pick(Pr,j);                                           
    if flag=='T'
      Nu(j)=NuTurb(R(j),Prj)+...
        (NuLam(Rc,Prj)-NuTurb(Rc,Prj)).*(Rc./R(j)).^n;        
    elseif flag=='H'
      Nu(j)=1./(1./NuTurb(R(j),Prj)+...
        (1./NuLam(Rc,Prj)-1./NuTurb(Rc,Prj)).*(Rc./R(j)).^n); 
    else
      error('Unknown flag, ''T'' or ''H'' required');
    end
  end
end
