function [Gammae,dGammae_dx,dGammae_dxdot,dGammae_du] = n_transport(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0
  %% DEFAULT PARAMETERS
  varargout{1} = RAPTORmodule('vpdn_manual');
  return %empty call, to get default structures
elseif nargin == 1
  error('deprecated use, use RAPTORmodules instead')
elseif nargin==6
  stap   = varargin{1};
  geop   = varargin{2};
  trap   = varargin{3};
  it     = varargin{4};
  model  = varargin{5};
  params = varargin{6};
else
  error('must call with 0 or 6 inputs');
end


% Initialize
dzz_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dzz_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
% Set to zero matrices
if nargout>2
  dDne_dx = dzz_dx;
  dDne_dxdot = dzz_dx;
  dDne_du = dzz_du;
  dVpe_dx = dzz_dx;
  dVpe_dxdot = dzz_dx;
  dVpe_du = dzz_du;
  
  dGammae_dxdot = dzz_dx;
  dGammae_du = dzz_du;
end

switch model.vpdn_e.name
  case 'vpdn_manual'
    % Dne is prescribed
    % Vpe is prescribed
    [Dne, Vpe] = vpdn_manual(stap,geop,trap,it,model,params.vpdn_e);
    calc_gamma = true;
  case 'vpdn_chiescal'
    % Dne is lineraly proportional to chi_e
    % Vpe is lineraly proportional to Dne
    [Dne, Vpe] = vpdn_chiescal(stap,geop,trap,it,model,params.vpdn_e);
    calc_gamma = true;
    
  case 'vpdn_benchmark'
    % From Y.S. Na ITPA-IOS benchmark effort
    [Dne, Vpe] = vpdn_benchmark(stap,geop,trap,it,model,params.vpdn_e);

    calc_gamma = true;
    
  case 'vpdn_MS'
    if nargout<=2
      [Dne,Vpe] = vpdn_MS(stap,geop,trap,it,model,params.vpdn_e);
    else
      [Dne,Vpe,dDne_dx,dDne_dxdot,dDne_du,dVpe_dx,dVpe_dxdot,dVpe_du] = vpdn_MS(stap,geop,trap,it,model,params.vpdn_e);
    end
    
    calc_gamma = true;
    
  case 'vpdn_QLKNN4Dkin'
    [tcoeff,dtcoeff_dx] = vpdn_QLKNN4Dkin(stap,geop,trap,it,model,params.vpdn_e);
    Dne = tcoeff.de;  Vpe = tcoeff.ve;
    dDne_dx = dtcoeff_dx.dde_dx; dVpe_dx = dtcoeff_dx.dve_dx;
    
    calc_gamma = true;
  case 'vpdn_QLKNN10D'
    [tcoeff,dtcoeff_dx] = vpdn_QLKNN10D(stap,geop,trap,it,model,params.vpdn_e);
    Dne = tcoeff.de;  Vpe = tcoeff.ve;
    dDne_dx = dtcoeff_dx.dde_dx; dVpe_dx = dtcoeff_dx.dve_dx;
%     Gammae = tcoeff.Gammae;  
%     dGammae_dx = dtcoeff_dx.dGammae_dx;

    Gammae = Vpe .* stap.ne;
    dGammae_dx =  bsxfun(@times, stap.dne_dx, Vpe) + bsxfun(@times, stap.ne, dVpe_dx);
    Vpe = zeros(size(tcoeff.ve));
    dVpe_dx = zeros(size(dtcoeff_dx.dve_dx));
    tcoeff.ve = Vpe;
    dtcoeff_dx.dve_dx = dVpe_dx;
    calc_gamma = false; % Gamma given directly
    
    %% check - recalculate Gamma from D,V
    [Gammae2] = calculate_Gamma(stap,geop,model,tcoeff.de,tcoeff.ve,...
      dtcoeff_dx.dde_dx,dtcoeff_dx.dve_dx,'e');
    
    disp(norm(Gammae2-Gammae)/norm(Gammae))
    plot(model.rgrid.rhogauss,[Gammae Gammae2]); 
    title('these should be the same');
    legend('Gamma direct','Gamma from D,V')
    drawnow
  otherwise
    error('not supported')
end

if calc_gamma
  % calculate fluxes from transport coefficients
  [Gammae,dGammae_dx,dGammae_dxdot,dGammae_du] = calculate_Gamma(stap,geop,model,Dne,Vpe,dDne_dx,dVpe_dx,'e');
end

end

function [Gamma,dGamma_dx,dGamma_dxdot,dGamma_du] = calculate_Gamma(stap,geop,model,Dn,Vp,dDn_dx,dVp_dx,species)
% calculates particle flux from transport coefficients

g1oVp2 = geop.g1./geop.Vp.^2;
g0oVp = geop.g0./geop.Vp;

facnep = -g1oVp2.*Dn;
facne  = g0oVp.*Vp;

switch species
  case 'e'
    nn = stap.ne;
    nnp = stap.nep;
    dnnp_dx = stap.dnep_dx;
    dnn_dx = stap.dne_dx;
  otherwise
    error('not yet implemented');
end
Gamma = facnep.*nnp + facne.*nn;

% Initialize
dzz_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dzz_du = zeros(model.rgrid.nrhogauss,model.dims.nu);
% Set to zero matrices
if nargout>1
  dfacnep_dx = bsxfun(@times,-g1oVp2,dDn_dx);
  dfacne_dx  = bsxfun(@times,g0oVp,dVp_dx);
  
  dGamma_dx = bsxfun(@times,facnep,dnnp_dx) + bsxfun(@times,nnp,dfacnep_dx)...
    + bsxfun(@times,facne,dnn_dx) + bsxfun(@times,nn,dfacne_dx);
  dGamma_dxdot = dzz_dx;
  dGamma_du = dzz_du;
end

end