function [varargout] = nui_star(ni,ti,q,Lref,eps,zion)
% nuistar with derivatives
%
% [nuistar] = nui_star(ni,ti,q,Lref,eps,zion)
% [nuistar,dnuistar_dni,dnuistar_dti,dnuistar_dq] = ...
%                 nuistar(ni,ti,q,Lref,eps,zion)
%
% computes nuistar according to formulas from O. Sauter et al, Phys. Plasmas 7 (1999) 2834.
% used then for neoclassical conductivity and Bootstrap current.
% Lref is used as generic length scale normalization (equal to R in Sauter's case)
%
% Outputs:
%    nuistar: ion collisionality
%
% Inputs:
%    ni()   : Local ion density
%    ti()   : Local ion temperature
%    q()    : Local safety factor
%    Lref() : Length scale. In Sauter paper Lref = Rgeom, geometrical center of local flux surface
%    eps()  : local inverse aspect ratio (a/R)
%    zion   : charge of main ion
% inputs must be column vectors

assert(all(ti>0),'negative ti input in nuistar')
assert(all(ni>0),'negative ni input in nuistar')

zion2 = zion.*zion;
lnLami = 30. - log(zion2.*zion) - 0.5*log(ni) + 1.5*log(ti);
eps32 = eps.*sqrt(eps); % eps^(3/2)
nuistar = 4.900E-18 .* q .* Lref .* ni .* (zion2.*zion2) .* lnLami ./ (ti.^2 .* eps32);

varargout{1} = nuistar;

%%
if nargout > 1
  % calculate derivatives
  
  dlnLami_dti = 3./(2*ti);
  dlnLami_dni = -1./(2.*ni);
  
  dnuistar_dti    = nuistar./lnLami.*dlnLami_dti - 2*nuistar./ti; % no dependence
  dnuistar_dni    = nuistar./ni + nuistar./lnLami.*dlnLami_dni;
  dnuistar_dq     = nuistar./q;
  %     dnuistar_dZion  = 4*nuistar./zion;
  
  varargout{2} = dnuistar_dni;
  varargout{3} = dnuistar_dti;
  varargout{4} = dnuistar_dq;
  
end

return