% Annular fin resistance 
function [Rfin, Vol] = AnnFinApprox(h,k,ri,ro,w,tip)
  Bi=h.*w./k;
  if max(Bi) > 0.1
    warning = 'AnnFinApprox: Bi > 0.1'
    Bi
  end
  m=sqrt(2*h.*ri^2./k./w);
  mro=m*ro/ri;
  if tip=='c'   % convective tip
    D1=2*ri.*besselk(1,mro)-m.*w.*besselk(0,mro);
    D2=2*ri.*besseli(1,mro)+m.*w.*besseli(0,mro);
    denom=2*ri.*(besseli(1,mro).*besselk(0,m)+besseli(0,m).*besselk(1,mro))+ ...
          m.*w.*(besseli(0,mro).*besselk(0,m)-besseli(0,m).*besselk(0,mro));
  else % tip=='a' % adiabatic tip
    D1=besselk(1,mro);
    D2=besseli(1,mro);
    denom=besseli(0,m).*besselk(1,mro)+besseli(1,mro).*besselk(0,m);
  end
  dT1=m.*(D1.*besseli(1,m)-D2.*besselk(1,m))./denom;   
  Rfin=1/(2*pi*k*w*(-dT1));
  V1=2*pi*ri^2*w./m*(D1*((ro/ri)*besseli(1,mro)-besseli(1,m))- ...
                     D2*((ro/ri)*besselk(1,mro)-besselk(1,m)))./denom;
  Vol=pi*(ro.^2-ri.^2)*w;
  Vol=[V1,Vol-V1];
end
