function [varargout] = neos(varargin)
% neos returns <jbs dot B>, sigma_neo, and the Jacobians if requested
% Using Sauter-Angioni model
% Physical formulas based on O. Sauter, C. Angioni, Phys. of Plas.,
% Vol.6, No 7, (1999)
%
% [out1] = neos(stap,geop,model,module_params)
%   params = neos; % get default parameters
%   [jbsB,signeo,djbs_dx,dsigneo_dx] = neos(stap,geop,model,params.neos) % run module

%#codegen
coder.extrinsic('check_gradients');
coder.extrinsic('check_BSfrac_gradients');
coder.extrinsic('param_sens');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing

if nargin == 0,
    %% DEFAULT PARAMETERS
    module_params = struct(...
        'signeo_fact', 1, ...      % scaling factor for neoclassical conductivity
        'jbsB_fact', 1, ...        % scaling factor for bootstrap current
        'sigma_csawtooth', 0, ...  % amplitude of reduction of signeo in the region q<1, off = 0, value lies within [0 1]. 1 is maximal reduction.
        'sigma_wsawtooth', 30,...  % 1/width of sawtooth effect
        'check',false,  ...        % optional self-check of gradients
        'implicit',true...
        );
    
    varargout{1} = module_params;
    return % empty call, probably to get default structures
elseif nargin==4; % change this depending on number of inputs to module
    stap = varargin{1};
    geop = varargin{2};
    model = varargin{3};
    neos_params = varargin{4}; % BS and sigmaneo parameters
else
    error('must call with 0 or 4 inputs');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the actual work here

nrg = numel(model.rgrid.rhogauss);
%assert(nrg<400);

% assign profiles at this time
% q and iota
diota_dx = stap.diota_dx;
q        = stap.q;
q2       = q.^2;
% Psi
drho_dpsi = 1./stap.psip;
dpsip_dx = stap.dpsip_dx;
% Te
te       = stap.te;
dte_dx   = stap.dte_dx;
tep      = stap.tep;
dtep_dx  = stap.dtep_dx;
% Ti
ti       = stap.ti;
dti_dx   = stap.dti_dx;
tip      = stap.tip;
dtip_dx  = stap.dtip_dx;
% Ne
ne       = stap.ne;
dne_dx   = stap.dne_dx;
nep      = stap.nep;
dnep_dx  = stap.dnep_dx;
% Ni
ni       = stap.ni;
dni_dx   = stap.dni_dx;
nip      = stap.nip;
dnip_dx  = stap.dnip_dx;
% Ze
ze       = stap.ze;
% Geometry
F        = geop.F;
epsilon  = geop.epsilon;
Rgeom    = geop.Rgeom;
delta    = geop.delta;

% Calculation of ftrap
ftrap = calc_ftrap(epsilon,delta);

[nuestar,dnuestar_dne,dnuestar_dte,~,dnuestar_dq] = nue_star(ne,te,ze,q,Rgeom,epsilon);
[nuistar,dnuistar_dni,dnuistar_dti,dnuistar_dq] = nui_star(ni,ti,q,Rgeom,epsilon,model.atom.Zi);
  
% avoid negative numbers which could give imag results
q = abs(q); te = abs(te); ti = abs(ti); ne = abs(ne); ni = abs(ni);

%%%%%%%%%%%%%%%%%%%%%%
%% jBS
%%%%%%%%%%%%%%%%%%%%%%
% every time if not constant BS coeff
% collisionality

% BS coefficients
[L31, L32, L34, alfa, dL31_dnuestar, dL32_dnuestar, dL34_dnuestar, dalfa_dnuistar] = ...
    BScoeff(ftrap, ze, nuestar, nuistar);

% 2*pi factor is because the Sauter-Angioni equation is for psi_chease
cc=2*pi*1.602e-19*F;
necoeff = cc.*L31;
nicoeff = cc.*L31;
tecoeff = cc.*(L31 + L32);
ticoeff = cc.*(L31 + alfa.*L34);

lntep = tep./te; % lntep = dlnte/drho = 1/te * dte/drho
lntip = tip./ti;
lnnep = nep./ne;
lnnip = nip./ni;

p_e = ne.*te;
p_i = ni.*ti;

% first jBS equation in RAPTOR doc
% NB: necoeff contain 2pi*T*1.602e-19 terms

% default equation with independent ne and Te
jbsB = -neos_params.jbsB_fact.*drho_dpsi.*...
    (p_e.*necoeff.*lnnep +...
     p_i.*nicoeff.*lnnip +...
     p_e.*tecoeff.*lntep +...
     p_i.*ticoeff.*lntip);

varargout{1} = jbsB;

%%%%%%%%%%%%%%%%%%%%%%
%% sigma_neo
%%%%%%%%%%%%%%%%%%%%%%
% Take sawteeth behaviour into account by decreasing the signeo in
% the region where q<1: this will increase q in this region (This results
% in a calculated q profile which is in better agreement with
% observations). signeo decreases -> I_pl decreases -> q increases
if neos_params.sigma_csawtooth ~= 0
    sawteeth_exp = exp(neos_params.sigma_wsawtooth * (0.95 - q) ); % q - 0.95 is done instead of q - 1, to influence the region q>1 less
    one_over_sawteeth_exp_plus_1_squared = 1./((1+sawteeth_exp).^2);
    sawtooth_term = neos_params.sigma_csawtooth*1./( 1 + sawteeth_exp ) + (1 - neos_params.sigma_csawtooth);
    if nargout > 2
        dsawtooth_term_dPsihat =  bsxfun(@times,-neos_params.sigma_csawtooth*neos_params.sigma_wsawtooth*...
            one_over_sawteeth_exp_plus_1_squared.*sawteeth_exp.* q2, diota_dx(:,model.psi.xind));
    end
else
    % must define on all execution paths
    sawteeth_exp = ones(model.rgrid.nrhogauss,1);
    sawtooth_term = ones(model.rgrid.nrhogauss,1);
    one_over_sawteeth_exp_plus_1_squared = ones(model.rgrid.nrhogauss,1);
    if nargout > 2
        dsawtooth_term_dPsihat =  zeros(model.rgrid.nrhogauss,model.psi.nsp);
    end
end

% Compute derivatives of sawtooth term w.r.t. model parameters
if nargout == 7 && ~isempty(neos_params.paramderiv)
    % construct structure signeo_var that contains variables
    % to compute derivatives of signeo w.r.t. model parameters
    signeo_var.dsawtooth_term_dcsawtooth = (1./( 1 + sawteeth_exp ) - 1 );
    signeo_var.dsawtooth_term_dwsawtooth = (-neos_params.sigma_csawtooth * one_over_sawteeth_exp_plus_1_squared .* ...
        sawteeth_exp .* (0.95 - q) );
    signeo_var.sawtooth_term = sawtooth_term;
end

% Spitzer conductivity
NZ = 0.58 + 0.74 ./ (0.76 + ze);
lnLam = 31.3 - log(sqrt(ne)./te);
sigsptz = 1.9012e+04.*te.*sqrt(te)./ze./NZ./lnLam;
Zeff32 = ze.*sqrt(ze); %Zeff^(3/2);
% Neoclassical correction
ft33eff = ftrap./(1. + (0.55-0.1.*ftrap).*sqrt(nuestar) + 0.45.*(1.-ftrap).*nuestar./(Zeff32));
signeo = (1. - ft33eff.*(1.+0.36./ze - ft33eff.*(0.59./ze - 0.23./ze.*ft33eff)));
sigmaneo = neos_params.signeo_fact * sawtooth_term .* sigsptz .*signeo;

varargout{2} = sigmaneo;

% For jBS: two extra output profiles (L31 and f_bse) for use in f_ntm.m
if nargout > 2
    % L_31 is approximately L_bs
    varargout{3} = L31;
    % Calculate the profile of the fraction of bootstrap current due to electrons
    part_e=p_e.*(necoeff.*lnnep+tecoeff.*lntep);
    part_i=p_i.*(necoeff.*lnnep+ticoeff.*lntip);
    f_bse=part_e./(part_e+part_i);
    varargout{4} = f_bse;
end
    
%%%%%%%%%%%%%%%%%%%%%%
%% Jacobians
%%%%%%%%%%%%%%%%%%%%%%
if nargout > 4
    % jBS gradients: djbs/dTe * dTe/dTehat + djbs/dTeprime *
    % dTeprime/dTehat etc
    djbsB_dx = zeros(nrg, model.dims.nx); % init
    
    % bootstrap coefficients
	dnecoeff_dnuestar = bsxfun(@times,cc,dL31_dnuestar);
    dnicoeff_dnuestar = bsxfun(@times,cc,dL31_dnuestar);
	dtecoeff_dnuestar = bsxfun(@times,cc,(dL31_dnuestar+dL32_dnuestar));
	dticoeff_dnuestar = bsxfun(@times,cc,dL31_dnuestar) + bsxfun(@times,cc.*alfa,dL34_dnuestar);
	dticoeff_dnuistar = bsxfun(@times,cc.*L34,dalfa_dnuistar);
    
    % Over nue_star
    djbsB_dnuestar =  -neos_params.jbsB_fact.*drho_dpsi.*...
        (p_e.*dnecoeff_dnuestar.*lnnep + ...
         p_i.*dnicoeff_dnuestar.*lnnip + ...
         p_e.*dtecoeff_dnuestar.*lntep + ...
         p_i.*dticoeff_dnuestar.*lntip);
    % Over nui_star
    djbsB_dnuistar = -neos_params.jbsB_fact.*drho_dpsi.*(p_i.*dticoeff_dnuistar.*lntip);
    
    % dsigneo_dx:
    dsigmaneo_dx = zeros(model.rgrid.nrhogauss,model.dims.nx); % init
    % Calculation derivatives you will need to compute dsigneo_dx
    
    dsigneo_dfteff   = -1 +(-0.36+(1.18-0.69.*ft33eff).*ft33eff)./ze;
    dfteff_dnuestar   = -ft33eff.^2./ftrap.*...
        ((0.55-0.1.*ftrap)./2./sqrt(nuestar) + 0.45*(1.-ftrap)./Zeff32);
    
    %%%%%%%%%%%%%%%%%
	if strcmp(model.psi.method,'state')
        % jBS
        diota_dpsihat = diota_dx(:,model.psi.xind);        
        dnuestar_dpsihat = bsxfun(@times,dnuestar_dq.*-q.^2,diota_dpsihat);
        dnuistar_dpsihat = bsxfun(@times,dnuistar_dq.*-q.^2,diota_dpsihat);
        % djbs/dpsihat
        % element-wise: djbs(i)/dpsih(j) where jbs(i) = f(i)/(C_i^T*psih)
        % = f(i)/(C_i^T*psih)^2 * Cij = jbs(i)/(C_i^T*psih) .* Cij
        % (+ nustar terms)       
        % Over x
        djbsB_dx(:,model.psi.xind) = bsxfun(@times,-jbsB.*drho_dpsi,dpsip_dx(:,model.psi.xind)) + ...
            bsxfun(@times,djbsB_dnuestar,dnuestar_dpsihat) + ...
            bsxfun(@times,djbsB_dnuistar,dnuistar_dpsihat);
        
        % sigma_neo
        dsigneo_dq = dsigneo_dfteff.*dfteff_dnuestar.*dnuestar_dq;        
        % Over x
        dsigmaneo_dx(:,model.psi.xind) = ...
            bsxfun(@times,-q.^2.*sigmaneo./signeo .*dsigneo_dq, diota_dx(:,model.psi.xind)) + ...
            bsxfun(@times,sigmaneo./sawtooth_term,dsawtooth_term_dPsihat);
    end
    
    %%%%%%%%%%%%%%%%
	if strcmp(model.te.method,'state')
        % jBS
        % Over Te
        djbsB_dTe = -neos_params.jbsB_fact.*drho_dpsi.*(ne.*necoeff.*lnnep + ...
            ne.*tecoeff.*lntep - ...
            p_e.*tecoeff.*tep./te.^2);
        % Over dTe/drho
        djbsB_dTeprime = -neos_params.jbsB_fact.*drho_dpsi.*(p_e.*tecoeff./te);
        % Over nue_star
        dnuestar_dtehat  = bsxfun(@times,dnuestar_dte,dte_dx(:,model.te.xind));
        % Over x
        djbsB_dx(:,model.te.xind) = bsxfun(@times,djbsB_dTe,dte_dx(:,model.te.xind))+...							   
                                   bsxfun(@times,djbsB_dTeprime,dtep_dx(:,model.te.xind))+...
                                   bsxfun(@times,djbsB_dnuestar,dnuestar_dtehat);
                               
        % sigma_neo
        dsigneo_dte   = dsigneo_dfteff.*dfteff_dnuestar.*dnuestar_dte;   
        dsigsptz_dte  = sigsptz.*(3./2./te - 1./lnLam./te);
        dsigmaneo_dte = sigmaneo./sigsptz.*dsigsptz_dte + sigmaneo./signeo .* dsigneo_dte;        
        % Over x
        dsigmaneo_dx(:,model.te.xind)  = bsxfun(@times,dsigmaneo_dte,dte_dx(:,model.te.xind));
    end
    
    %%%%%%%%%%%%%%%%
    % Ti scaled from Te	or state	
    if strcmp(model.ti.method,'tescal')||strcmp(model.ti.method,'state')
        % jBS
        % Over Ti
        djbsB_dTi = -neos_params.jbsB_fact.*drho_dpsi.*(ni.*nicoeff.*lnnip + ...
            ni.*ticoeff.*lntip - p_i.*ticoeff.*tip./ti.^2);
        % Over dTi_drho
        djbsB_dTiprime = -neos_params.jbsB_fact.*drho_dpsi.*(p_i.*ticoeff./ti);
        % Over nui_star
        dnuistar_dx  = bsxfun(@times,dnuistar_dti,dti_dx);
        % Over x
        djbsB_dx = djbsB_dx +...
            bsxfun(@times,djbsB_dTi,dti_dx) + ...
            bsxfun(@times,djbsB_dTiprime,dtip_dx) + ...
            bsxfun(@times,djbsB_dnuistar,dnuistar_dx);	
    end
    
    %%%%%%%%%%%%%%%%
	if strcmp(model.ne.method,'state')
        % jBS
        % Over Ne
        djbsB_dNe = -neos_params.jbsB_fact.*drho_dpsi.*(te.*necoeff.*lnnep + ...
            te.*tecoeff.*lntep - p_e.*necoeff.*nep./ne.^2);
        % Over dNe/drho
        djbsB_dNeprime = -neos_params.jbsB_fact.*drho_dpsi.*(p_e.*necoeff./ne);
        % Over nue_star
        dnuestar_dnehat  = bsxfun(@times,dnuestar_dne,dne_dx(:,model.ne.xind));
        % Over x
        djbsB_dx(:,model.ne.xind) = bsxfun(@times,djbsB_dNe,dne_dx(:,model.ne.xind))+...							   
                                   bsxfun(@times,djbsB_dNeprime,dnep_dx(:,model.ne.xind))+...
                                   bsxfun(@times,djbsB_dnuestar,dnuestar_dnehat);
       
        % sigma_neo
        dsigneo_dne = dsigneo_dfteff.*dfteff_dnuestar.*dnuestar_dne;
        dsigsptz_dne = sigsptz./(lnLam.*2.*ne);
        dsigmaneo_dne = sigmaneo./sigsptz.*dsigsptz_dne + sigmaneo./signeo .* dsigneo_dne;        
        % Over x
        dsigmaneo_dx(:,model.ne.xind) = bsxfun(@times,dsigmaneo_dne,dne_dx(:,model.ne.xind));
    end
    
    %%%%%%%%%%%%%%%%
    % Ni scaled from Ne	or state	
    if (strcmp(model.ne.method,'state')&&strcmp(model.ni.method,'nescal'))||...
       (strcmp(model.ne.method,'state')&&strcmp(model.ni.method,'qnze'))||...
        strcmp(model.ni.method,'state')
        % jBS
        % Over Ni
        djbsB_dNi = -neos_params.jbsB_fact.*drho_dpsi.*(ti.*nicoeff.*lnnip + ...
            ti.*ticoeff.*lntip - p_i.*nicoeff.*nip./ni.^2);
        % Over dNi_drho
        djbsB_dNiprime = -neos_params.jbsB_fact.*drho_dpsi.*(p_i.*nicoeff./ni);
        % Over nui_star
        dnuistar_dx  = bsxfun(@times,dnuistar_dni,dni_dx);
        
        % Over x
        djbsB_dx = djbsB_dx+...
            bsxfun(@times,djbsB_dNi,dni_dx) + ...
            bsxfun(@times,djbsB_dNiprime,dnip_dx) + ...
            bsxfun(@times,djbsB_dnuistar,dnuistar_dx);
    end
  
    varargout{5} = djbsB_dx;
    varargout{6} = dsigmaneo_dx;   
      
    % dsigneo_dZeff
    % Computation of sigma_neo derivative w.r.t Zeff
    if nargout==7
        dfteff_dZeff = dfteff_dnuestar.*dnuestar_dZeff + ...
            ftrap./(1+(0.55-0.1.*ftrap).*sqrt(nuestar) + 0.45.*(1-ftrap).*nuestar./Zeff32).^2.* ...
            (1.5*(0.4591-ftrap).*nuestar./(Zeff32.*ze));
        dsigneo_dZeff    = -dfteff_dZeff + ...
            (-0.36 + 0.59*ft33eff-0.23*ft33eff.^2).*((dfteff_dZeff.*ze-ft33eff)./ze.^2) ...
            + ft33eff./ze.*(0.59.*dfteff_dZeff-0.46.*ft33eff.*dfteff_dZeff);
        dsigsptz_dZeff   = -sigsptz.*(1./ze - 1./NZ.*(0.74./(0.76 +ze).^2));
        dsigmaneo_dZeff  = sigmaneo./sigsptz .* dsigsptz_dZeff + sigmaneo./signeo .* dsigneo_dZeff;

        varargout{7} = dsigmaneo_dZeff;

    end
    
    % Calculation derivatives signeo w.r.t model parameters
    if nargout == 7 && ~isempty(sn.paramderiv)
        % Calculation dsigneo_dZeff
        signeo_var.dsigma_dZeff = dsigmaneo_dZeff;
        varargout{7} = param_sens(sigmaneo,signeo_var,model,sn,x,g,v,it);
    end
end

return

%function [L31, L32, L34, alfa] = BScoeff(ft,varargin)
function [varargout] = BScoeff(ft, zeff, nuestar, nuistar)

% FROM NEOS FUNCTION

sqnuestar=sqrt(nuestar);
ft2 = ft.*ft;
ft6 = ft2.*ft2.*ft2;

ft31eff = ft ./ (1.+(1.-0.1.*ft).*sqnuestar + 0.5.*(1.-ft).*nuestar./zeff);
ft32ee_eff = ft ./ (1. + 0.26.*(1.-ft).*sqnuestar + 0.18.*(1.-0.37.*ft).*nuestar./sqrt(zeff));
ft32ei_eff = ft ./ (1. + (1.+0.6.*ft).*sqnuestar + 0.85.*(1.-0.37.*ft).*nuestar.*(1.+zeff));
ft34eff = ft ./ (1.+(1.-0.1.*ft).*sqnuestar + 0.5.*(1.-0.5.*ft).*nuestar./zeff);
alfa0 = - 1.17.*(1.-ft) ./ (1.-0.22.*ft-0.19.*ft.^2);

% quick calculation for often-used products
ft32ee_eff2 = ft32ee_eff.*ft32ee_eff;
ft32ee_eff4 = ft32ee_eff2.*ft32ee_eff2;
ft32ei_eff2 = ft32ei_eff.*ft32ei_eff;
ft32ei_eff3 = ft32ei_eff2.*ft32ei_eff;
ft32ei_eff4 = ft32ei_eff2.*ft32ei_eff2;

dft31eff_dnuestar = -0.5 .* ft31eff.*ft31eff .*((1. - 0.1*ft) ./ sqnuestar + (1. - ft) ./ zeff) ./ ft;
dft32ee_eff_dnuestar = -ft32ee_eff2 .* (0.13 .* (1  -     ft) ./ sqnuestar + 0.18 .* (1. - 0.37 .*ft) ./ sqrt(zeff)) ./ ft;
dft32ei_eff_dnuestar = -ft32ei_eff2 .* (0.5 .*  (1. + 0.6*ft) ./ sqnuestar + 0.85 .* (1. - 0.37 .* ft) .* (1. + zeff)) ./ ft;
dft34eff_dnuestar = -0.5 .* ft34eff.*ft34eff .*((1. - 0.1*ft) ./ sqnuestar + (1. - 0.5 .* ft) ./ zeff) ./ ft;

% coefficients
zeffp1 = zeff+1.;

L31 = ft31eff .* ( (1.+1.4./zeffp1) ...
    + ft31eff .* (-1.9./zeffp1 + ft31eff .* (0.3./zeffp1 + 0.2./zeffp1 .* ft31eff)));
dL31_dnuestar  = (1. + (1.4 + ft31eff .* (-3.8 + ft31eff .* (0.9 + 0.8 .* ft31eff))) ./ zeffp1) .* dft31eff_dnuestar;

C1 = (0.05 + 0.62 .* zeff) ./ zeff ./ (1. + 0.44 .* zeff);
C2 = 1. ./ (1. + 0.22 .* zeff);
%C3 = -1.2 .* (1 + 0.5 .* zeff);
%F32ee = ft32ee_eff .* (C1 + ft32ee_eff .* ( C2 + ft32ee_eff .* (-1.2 .* C2 + ft32ee_eff .* (-C1 +0.2 .* C2 +C3))));
% dF32ee_dft32ee_eff = C1 + ft32ee_eff .* ( 2 .* C2 + ft32ee_eff .* (-3.6 .* C2 + ft32ee_eff .* 4. .* (- C1 + 0.2 .* C2 + C3)));

F32ee = C1.*(ft32ee_eff-ft32ee_eff2.*ft32ee_eff2) + ...
    C2.*ft32ee_eff2.*(1 - 1.2*ft32ee_eff + 0.2*ft32ee_eff2);
dF32ee_dft32ee_eff = C1.*(1 - 4*ft32ee_eff2.*ft32ee_eff) + ...
    C2.*(2*ft32ee_eff  -3.6*ft32ee_eff2 + 0.8*ft32ee_eff2.*ft32ee_eff) ;

C4 = - (0.56 + 1.93 .* zeff) ./ (zeff .* (1. + 0.44 .* zeff));
C5 = 4.95 ./ (1. + 2.48 .* zeff);
C6 = 1.2 ./ (1.+0.5.*zeff);
%F32ei = ft32ei_eff .* (C4 + ft32ei_eff .* (C5 + ft32ei_eff .* (-0.55 .* C5 + ft32ei_eff .* (-C4 -0.45 .* C5 -C3))));
%dF32ei_dft32ei_eff = C4 + ft32ei_eff .* (2. .* C5 + ft32ei_eff .* (-1.65 .* C5 + ft32ei_eff .* 4. .* (-C4 -0.45 .* C5 -C3)));

F32ei = C4 .* (ft32ei_eff-ft32ei_eff4) ...
    + C5.* ft32ei_eff2.*(1.-0.55.*ft32ei_eff-0.45.*ft32ei_eff2) ...
    + C6 .* (ft32ee_eff4-ft32ei_eff4);

dF32ei_dft32ei_eff = C4.*(1-4*ft32ei_eff3) ...
    + C5.*(2*ft32ei_eff - 1.65*ft32ei_eff2 - 1.8*ft32ei_eff3) ...
    + C6 .* (-4*ft32ei_eff3);

dF32ei_dft32ee_eff = C6.*(4*ft32ee_eff.^3);

L32 = (F32ee + F32ei);
dL32_dnuestar = (dF32ee_dft32ee_eff + dF32ei_dft32ee_eff) .* dft32ee_eff_dnuestar ...
    + dF32ei_dft32ei_eff .* dft32ei_eff_dnuestar;

L34 = ft34eff.* ( (1.+1.4./zeffp1) ...
    - ft34eff.*(1.9./zeffp1-ft34eff.*(0.3./zeffp1+0.2./zeffp1.*ft34eff)) );
dL34_dnuestar = (1. + (1.4 + ft34eff .* (-3.8 + ft34eff .* (0.9 + 0.8 .* ft34eff))) ./ zeffp1) .* dft34eff_dnuestar;

sqnui = sqrt(nuistar);
nui2ft6 = nuistar.^2 .* ft6;
A1 = (alfa0 + 0.25 .* (1. - ft2) .* sqnui) ./ (1. + 0.5 .* sqnui);
dA1_dnuistar = ((1 - ft2) - 2 .* alfa0) ./ ( 8 .* sqnui .* (1 + 0.5 .* sqnui) .^2);
A2 = 0.315 .* nui2ft6;
dA2_dnuistar = 0.63 .* ft6 .* nuistar;
A3 =  1 ./ (1. + 0.15 .* nui2ft6);
dA3_dnuistar = - A3 .^2 .* (0.3 .* nuistar .* ft6);
alfa = (A1 + A2) .* A3;
dalfa_dnuistar = (dA1_dnuistar + dA2_dnuistar) .* A3 + (A1 + A2) .* dA3_dnuistar;

varargout{1} = L31;
varargout{2} = L32;
varargout{3} = L34;
varargout{4} = alfa;

if nargout == 8
    varargout{5} = dL31_dnuestar;
    varargout{6} = dL32_dnuestar;
    varargout{7} = dL34_dnuestar;
    varargout{8} = dalfa_dnuistar;
    
    % optional check
%     check = false;
%     if check
%         check_BSfrac_gradients(ft,zeff,nuestar,nuistar,L31,L32,L34,alfa,...
%             dL31_dnuestar,dL32_dnuestar,dL34_dnuestar,dalfa_dnuistar);
%     end
end

return

%%
function check_BSfrac_gradients(ftrap,Zeff,nuestar,nuistar,L31,L32,L34,alfa,...
    dL31_dnuestar,dL32_dnuestar,dL34_dnuestar,dalfa_dnuistar)


figname = 'BSfrac_check';
hf = findobj('name',figname);
if isempty(hf)
    hf = figure('name',figname);
end

% checkf or BSfrac function
dnuestar = 1e-6*rand;
dnuistar = 0e-6*rand;
[L31_d, L32_d, L34_d, alfa_d] = ...
    BScoeff(ftrap,Zeff, nuestar+dnuestar, nuistar);

%%
dL31 = L31_d-L31;
dL32 = L32_d-L32;
dL34 = L34_d-L34;
dalfa = alfa - alfa_d;

subplot(4,2,1,'parent',hf); title('dL31/dnuestar')
plot([dL31,dL31_dnuestar*dnuestar]);
subplot(4,2,2,'parent',hf);
plot([dL31-dL31_dnuestar*dnuestar]);

subplot(4,2,3,'parent',hf); title('dL32/dnuestar')
plot([dL32,dL32_dnuestar*dnuestar]);
subplot(4,2,4,'parent',hf);
plot([dL32-dL32_dnuestar*dnuestar]);

subplot(4,2,5,'parent',hf); title('dL34/dnuestar')
plot([dL34,dL34_dnuestar*dnuestar]);
subplot(4,2,6,'parent',hf);
plot([dL34-dL34_dnuestar*dnuestar]);

subplot(4,2,7,'parent',hf); title('dalfa/dnuestar')
plot([dalfa,dalfa_dnuistar*dnuistar]);
subplot(4,2,8,'parent',hf);
plot([dalfa-dalfa_dnuistar*dnuistar]);
drawnow;
%%
return

return