function varargout = chi_BgB(varargin)

% Bohm-gyroBohm transport model
%   params = chi_BGB % get default params
%   [chie,chii,dchie_dx,dchii_dx,dchie_dp] = chi_BgB(x,g,v,model,module_params_in) ; % run
%
% In this version, the electron pressure gradient length has been rewritten to
%       gradient lengths of Te and ne

%--------------------------------------------------------------
if nargin == 0,
  varargout{1} = 'chi_e'; % type;
  varargout{2} = [];
  return
elseif nargin==2
   mm.name = mfilename;
   %% DEFAULT PARAMETERS
    
    default_params = struct( ...
        'aeb',      8e-5,    ...      % coefficient of Bohm term of chie [8e-5]
        'aegb',     7e-2,    ...      % coefficient of gyro Bohm term of chie [7e-2]
        'aib',      1.6e-4,  ...      % coefficient of Bohm term of chii [1.6e-4]
        'aigb',     1.75e-2, ...      % coefficient of gyro Bohm term of chii [1.75e-2]
        'cneo',     1e-3,     ...     % Adding constant neo-classical term
        'chi_csawtooth', 0, ...       % amplitude of increased transport in core due to sawtooth,  0 == off
        'chi_wsawtooth', 0.1,...        % width(q) of sawtooth effect
        'chi_ccentral', 2, ...        % amplitude of central extra chie, 0 = off
        'chi_wcentral', 0.1, ...        % half-width(rho) of central extra chie
        'chimax',   20,    ...        % maximum value of chie and chii
        'a1',       2,       ...      % atomic mass of main ion, default D
        'z1',       1,       ...      % charge of main ion, default D
        'implicit', true,...          % choose whether to treat as implicit in solver
        'check',    false   ...       % optional check of derivatives
        );
    
      varargout{1} = mm;
      varargout{2} = default_params;
    return %empty call, probably to get default structures
elseif nargin==6 % change this depending on number of inputs to module
  stap = varargin{1};
  %geop = varargin{2};
  %trap = varargin{3};
  %it = varargin{4};
  model = varargin{5};
  pp = varargin{6}; % distribute inputs
else
    error('must call with 0 or 6 inputs');
end

%%% BgB model

% profiles
te       = stap.te;
dte_dx   = stap.dte_dx;
ne       = stap.ne;
tep      = stap.tep;
dtep_dx  = stap.dtep_dx;
nep      = stap.nep;
diota_dx = stap.diota_dx;
q = stap.q;       
q_squared   = q.^2;

% avoid negative sqrt(imag numbers)
te = abs(te); ne = abs(ne);

% other quantities
Bt      = model.equi.B0;
aeff    = model.equi.epsilon*model.equi.R0*sqrt(model.equi.kappa);     % aeff = a*delta^(1/2)

% gradients
gte     = tep./aeff;                       % Grad(Te)
gne     = nep./aeff;                       % Grad(ne)

% Might need adaption

% rho_star (of main ion)
rho_star = 1.02e-4 .* sqrt(pp.a1 .* te) ./ (pp.z1 .* aeff .* Bt);
gh = 2;

% dimensionless 1/Lpe, i.e. a* grad(pe) / pe

invLte  = abs(gte) ./ te;
invLne  = abs(gne) ./ ne;

% chii, chie bohm

xibte     = te ./ Bt .* q_squared .* gh.*aeff.*invLte;
xibne     = te ./ Bt .* q_squared .* gh.*aeff.*invLne;
xib       = xibte + xibne;

% chii, chie gyro bohm

xigb    = rho_star  ./ Bt .* abs(gte) * aeff;

% transport coefficients: sum of Bohm and gyroBohm parts, limited to maximum

% optional plots
%figure(8); clf;
%plot(model.rgrid.rhogauss,cp.aegb*xigb,'b',model.rgrid.rhogauss,cp.aeb*xib,'r');
%legend('gb','b'); drawnow;
chie0                        = pp.cneo + pp.aegb * xigb + pp.aeb * xib;
chii0                        = pp.cneo + pp.aigb * xigb + pp.aib * xib;

%% Specific fixes
[sawtooth_term,dsawtoothterm_diotagauss] = chi_sawtooth(q,q_squared,model,pp);
chicentral = chi_central(model,pp);

%% total chi
chie = chie0 + sawtooth_term + chicentral;
chii = chii0 + sawtooth_term + chicentral;

%% substitute chi values with max chi when exceeding bound
% indices where normal chie
ie                          = (chie < pp.chimax);
ii                          = (chii < pp.chimax);

chie(~ie)                   = pp.chimax;
chii(~ii)                   = pp.chimax;

varargout{1} = chie;
varargout{2} = chii;

%%
if nargout>2
    % Initialize
    dchii_dx            = zeros(model.rgrid.nrhogauss,model.dims.nx);
    dchie_dx            = zeros(model.rgrid.nrhogauss,model.dims.nx);
    %% dchi_dx:
    % GRADIENTS w.r.t. Tehat and Psihat
    % First, dchi_dTehat:
    dchie_b_Tep_dTehat          =  pp.aeb.*xibte(ie)./gte(ie);
    dchie_b_Te_dTehat           =  pp.aeb.*xibne(ie)./te(ie);
    dchie_gb_Te_dTehat          =  .5*pp.aegb*xigb(ie)./(te(ie));
    dchie_gb_Tep_dTehat         =  pp.aegb*xigb(ie)./gte(ie);
    dchii_b_Tep_dTehat          =  pp.aib.*xibte(ii)./gte(ii);
    dchii_b_Te_dTehat           =  pp.aib.*xibne(ii)./te(ii);
    dchii_gb_Te_dTehat          =  .5*pp.aigb*xigb(ii)./(te(ii));
    dchii_gb_Tep_dTehat         =  pp.aigb*xigb(ii)./gte(ii);
    
    dchii_dx(ii,model.te.xind)           =  bsxfun(@times,dchii_b_Te_dTehat+dchii_gb_Te_dTehat,dte_dx(ii,model.te.xind)) + ...
        bsxfun(@times,dchii_b_Tep_dTehat+dchii_gb_Tep_dTehat,dtep_dx(ii,model.te.xind)./aeff);
    dchie_dx(ie,model.te.xind)           =  bsxfun(@times,sawtooth_term(ie).*(dchie_b_Te_dTehat+dchie_gb_Te_dTehat),dte_dx(ie,model.te.xind)) + ...
        bsxfun(@times,dchie_b_Tep_dTehat+dchie_gb_Tep_dTehat,dtep_dx(ie,model.te.xind)./aeff);
    
    % Second, dchi_dpsihat:
    %dchi_dpsihat:
    dchie0_diota = (-2*q(ie).*pp.aeb.*xib(ie));
    dchii0_diota = (-2*q(ii).*pp.aib.*xib(ii));
    
    dchie_diota  =  dchie0_diota + dsawtoothterm_diotagauss(ie);
    dchii_diota  =  dchii0_diota + dsawtoothterm_diotagauss(ii);
    
    dchie_dpsihat = bsxfun(@times,dchie_diota,diota_dx(ie,model.psi.xind));
    dchii_dpsihat = bsxfun(@times,dchii_diota,diota_dx(ii,model.psi.xind));
    
    dchie_dx(ie,model.psi.xind)      = dchie_dpsihat;
    dchii_dx(ii,model.psi.xind)      = dchii_dpsihat;
  
    varargout{3} =   dchie_dx;
    varargout{4} =   dchii_dx;
end

return