function [varargout] = nue_star(ne,te,zeff,q,Lref,eps)
% Cleaned up version of nustar with derivatives
%
% [nuestar] = nuestar(ne,te,zeff,q,Lref,eps)
% [nuestar,dnuestar_dne,dnuestar_dte,dnuestar_dzeff,dnuestar_q] = ...
%                        nuestar(ne,te,zeff,q,Lref,eps)
%
% computes nuestar 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:
%    nuestar: electron collisionality
%
% Inputs:
%    ne()   : Local electron density
%    te()   : Local electron temperature
%    ni()   : Local ion density
%    ti()   : Local ion temperature
%    Zeff() : effective charge
%    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)
% inputs must be column vectors

if any(ne<=0)
  error('negative ne input in nuestar')
end
if any(te<=0)
  error('negative te input in nuestar')
end

%%
lnLam = 31.3 - 0.5*log(ne) + log(te);
eps32 = eps.*sqrt(eps); % eps^(3/2)
nuestar = 6.921E-18 .* q .* Lref .* ne .* zeff    .* lnLam ./ (te.^2 .* eps32);

varargout{1} = nuestar;

%%
if nargout > 1
    % calculate derivatives
    dlnLam_dne = -1./(2*ne);
    dlnLam_dte = 1./te;
    
    dnuestar_dne    = nuestar./ne + nuestar./lnLam.*dlnLam_dne;
    dnuestar_dte    = nuestar./lnLam.*dlnLam_dte  - 2*nuestar./te;
    dnuestar_dzeff  = nuestar./zeff;
    dnuestar_dq     = nuestar./q;
    
    varargout{2} = dnuestar_dne;
    varargout{3} = dnuestar_dte;
    varargout{4} = dnuestar_dzeff; 
    varargout{5} = dnuestar_dq;

  %%  
    
end

return