function trap = transport_profiles(stapk,stap,geop,dk,u,it,model,params,inewt,trap)

% function trap = transport_profiles(stapk,stap,geop,it,model,params,inewt,trap)
% stapk: state profiles from previous iteration (for explicit equations);
% stap: state profiles
% geop: geometry profiles

persistent chie0 chii0 dne0 vne0 signeo jbsB jreB ntmfact dwntmdt

ngauss = model.rgrid.nrhogauss;
nx = model.dims.nx;
nu = model.dims.nu;
zz      = zeros(ngauss,1);
dzz_dx  = zeros(ngauss,nx);
dzz_du  = zeros(ngauss,nu);
zn = zeros(model.ntm.n_ntm,1);

if isempty(chie0);   chie0   = zz;    end
if isempty(chii0);   chii0   = zz;    end
if isempty(dne0);    dne0    = zz;    end
if isempty(vne0);    vne0    = zz;    end
if isempty(signeo);  signeo  = zz;    end
if isempty(jbsB);    jbsB    = zz;    end
if isempty(ntmfact); ntmfact = zz;    end
if isempty(jreB);    jreB    = zz;    end
if ~isempty(zn)
  if isempty(dwntmdt), dwntmdt = zn; end
end

% Bootstrap current and Sigma - neoclassical conductivity
if params.neos.implicit
  [jbsB,signeo,trap.Lbs,trap.fbse,djbsB_dx,dsigneo_dx] = neos(stap,geop,model,params.neos);
else % explicit
  if inewt == 0
    [jbsB,signeo,trap.Lbs,trap.fbse] = neos(stapk,geop,model,params.neos);
  end
  djbsB_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dsigneo_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
end
trap.jbsB = jbsB;
trap.djbsB_dx = djbsB_dx;
trap.signeo = signeo;
trap.dsigneo_dx = dsigneo_dx;

%% Auxiliary heating and current drive sources
% EC heating and current drive
[trap.pec, trap.jec,dpec_dx, djec_dx, dpec_du, djec_du] = echcd(stap,geop,u,it,model,params.echcd);
% NB heating and current drive
[trap.pnbe,trap.pnbi,trap.jnb,trap.pfastnb,dpnbe_dx,dpnbi_dx,djnb_dx, ~, dpnbe_du,dpnbi_du,djnb_du] = nbhcd(stap,geop,u,it,model,params.nbhcd);
% LH heating and current drive
[trap.plh, trap.jlh, dplh_dx, djlh_dx, dplh_du, djlh_du] = lhhcd(stap,geop,u,it,model,params.lhhcd);
% IC heating and current drive
[trap.pice,trap.pici,trap.jic,trap.pfastic,...
  dpice_dx,dpici_dx,djic_dx,~,...
  dpice_du,dpici_du,djic_du,~] = ichcd(stap,geop,u,it,model,params.ichcd);

% Runaway electron current
if inewt == 0
  jreB = runaways(stapk,geop,u,it,model,params.runaways);
end
trap.jreB = jreB;

% sum to form paux, jaux
trap.pauxe    = trap.pec + trap.pnbe + trap.plh + trap.pice;
dpauxe_dx = dpec_dx + dpnbe_dx + dplh_dx + dpice_dx;
dpauxe_du = dpec_du + dpnbe_du + dplh_du + dpice_du;

trap.pauxi    = trap.pnbi + trap.pici;
dpauxi_dx = dpnbi_dx + dpici_dx;
dpauxi_du = dpnbi_du + dpici_du;

trap.jauxB     = (trap.jec + trap.jnb + trap.jlh + trap.jic) * model.equi.B0;
trap.djauxB_dx = (djec_dx + djnb_dx  + djlh_dx + djic_dx)    * model.equi.B0;
trap.djauxB_du = (djec_du + djnb_du  + djlh_du + djic_du)    * model.equi.B0;

% disturbance current source
if numel(model.dist.psi.ind)>0
  trap.jdistB = model.dist.psi.Lamgauss*dk(model.dist.psi.ind)*model.equi.B0;
else
  trap.jdistB = zeros(model.rgrid.nrhogauss,1);
end

%% Sources
% Ohmic power
%eval_jphioR returns jtor/R0,
% since Poh = 1/2pi/R0 * joh*Upl, we get Poh = 1/2pi*<jphi/R>*Upl
% also since jtor = R0*<jphi/R> then Poh = jtor/(2*pi*R0) * Upl
% or: Poh = jtor2pr*Upl with jtor2pr = jphioR/2pi

kiVp =  1/(8*pi^2*4e-7*pi*geop.Phib)./geop.Vp;
psippart = (kiVp .* (geop.F.*geop.g23orp + geop.Fp.*geop.g23or) );
psipppart =(kiVp .* (geop.F.*geop.g23or) );


trap.jphioR = psippart.*stap.psip + psipppart.*stap.psipp;
trap.djphioR_dx = dzz_dx;
trap.djphioR_dx = bsxfun(@times,psippart,stap.dpsip_dx) + ...
  bsxfun(@times,psipppart,stap.dpsipp_dx);

jtor2pr = trap.jphioR/(2*pi);
djtor2pr_dx = trap.djphioR_dx/(2*pi);

jtorPsidot2pir = jtor2pr.*stap.upl;

trap.poh = abs(jtorPsidot2pir);
dpoh_dx    = bsxfun(@times,sign(jtorPsidot2pir).*stap.upl,djtor2pr_dx);
dpoh_dxdot = bsxfun(@times,sign(jtorPsidot2pir).*jtor2pr, stap.dupl_dxdot);

% other pouwer sources
[trap.palphae,trap.palphai,dpalphae_dx,dpalphai_dx]  = palpha(stap,geop,it,model,params.palpha);
[trap.pbrem,~,dpbrem_dx,~]                      = pbrem(stap,geop,it,model,params.pbrem);
[trap.pei,~,  dpei_dx,~]                        = pei(stap,geop,it,model,params.pei);
[trap.prad,~, dprad_dx,~]                       = prad(stap,geop,it,model,params.prad);

% fast ion pressure
trap.pfast = trap.pfastnb + trap.pfastic; % to add: fast ions from alphas

% disturbance power (for Kalman Filter)
if numel(model.dist.te.ind)>0
  trap.pdist = model.dist.te.Lamgauss*dk(model.dist.te.ind);
else
  trap.pdist = zeros(model.rgrid.nrhogauss,1);
end

%%% total powers
% to electrons
trap.petot        = trap.poh + trap.pauxe + trap.palphae  - trap.prad - trap.pbrem  - trap.pei  + trap.pdist;
trap.dpetot_dx    = dpauxe_dx   + dpoh_dx     + dpalphae_dx - dprad_dx - dpbrem_dx - dpei_dx;
trap.dpetot_dxdot = dpoh_dxdot;
trap.dpetot_du    = dpauxe_du;

% to ions
trap.pitot        = trap.pauxi + trap.pei + trap.palphai;
trap.dpitot_dx    = dpauxi_dx + dpei_dx + dpalphai_dx;
trap.dpitot_dxdot = zeros(model.rgrid.nrhogauss,model.dims.nx);
trap.dpitot_du    = dpauxi_du;

%% Heat transport model
if params.chi_e.implicit
  [chie0,chii0,...
    dchie0_dx,dchii0_dx,...
    dchie0_dxdot,dchii0_dxdot,...
    dchie0_du, dchii0_du,...
    ] = chi_e(stap,geop,trap,it,model,params);
else % explicit
  if inewt == 0
    [chie0,chii0] = chi_e(stap,geop,trap,it,model,params);
  end
  dchie0_dx    = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dchie0_dxdot = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dchie0_du    = zeros(model.rgrid.nrhogauss,model.dims.nu);
  dchii0_dx    = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dchii0_dxdot = zeros(model.rgrid.nrhogauss,model.dims.nx);
  dchii0_du    = zeros(model.rgrid.nrhogauss,model.dims.nu);
end

% hmode modification of chie
if params.hmode.active && strcmp(model.te.method,'state')
  [chie,dchie_dx,dchie_dxdot,dchie_du] = hmode(stap,geop,trap,model,params.hmode,'te',chie0,dchie0_dx,dchie0_dxdot,dchie0_du);
else
  [chie,dchie_dx,dchie_dxdot,dchie_du] = chie_defaults(chie0,dchie0_dx,dchie0_dxdot,dchie0_du);
end
  
if params.hmode.active && strcmp(model.ti.method,'state')
  [chii,dchii_dx,dchii_dxdot,dchii_du] = hmode(stap,geop,trap,model,params.hmode,'ti',chii0,dchii0_dx,dchii0_dxdot,dchii0_du);
else
  [chii,dchii_dx,dchii_dxdot,dchii_du] = chie_defaults(chii0,dchii0_dx,dchii0_dxdot,dchii0_du);
end

% ntm modification of chie
if params.ntm.active
  if params.ntm.implicit
    [ntmfact,dntmfact_dx] = ntm_chie(stap,geop,u,model,params);
  else % explicit
    if inewt == 0
      % treat NTM effect only in explicit form: use xk
      [ntmfact] = ntm_chie(stapk,geop,u,model,params);
    end
    dntmfact_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
  end
  % further factor for chie due to ntm
  chie = ntmfact.*chie;
  dchie_dx = bsxfun(@times,ntmfact,dchie_dx) + ...
    bsxfun(@times,chie,dntmfact_dx);
end

trap.chie = chie;
trap.dchie_dx = dchie_dx;
trap.dchie_dxdot = dchie_dxdot;
trap.dchie_du = dchie_du;

trap.chii     = chii;
trap.dchii_dx = dchii_dx;
trap.dchii_dxdot = dchii_dxdot;
trap.dchii_du = dchii_du;

%% Particle transport model
if strcmp(model.ne.method,'state')
  % Particle transport: electrons
  if ~isfield(params.vpdn_e,'implicit') || params.vpdn_e.implicit
    % implicit
    [dne0, vne0,...
      ddne0_dx, dvne0_dx,...
      trap.ddne_dxdot,  trap.dvne_dxdot,...
      trap.ddne_du, trap.dvne_du] =...
      vpdn_e(stap,geop,trap,it,model,params);
  else
    if params.numerics.LoDestro || inewt == 0
      % explicit
      [dne0,vne0] =...
        vpdn_e(stapk,geop,trap,it,model,params);
    end
    ddne0_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
    dvne0_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
  end
  [trap.sne, trap.dsne_dx, trap.dsne_dxdot, trap.dsne_du] = sn_e(it,model,params);
  
  % hmode modification of Dne,Vne
  if params.hmode.active
    [trap.dne,trap.vne,trap.ddne_dx,trap.dvne_dx] = ...
      hmode(stap,geop,trap,model,params.hmode,'ne',dne0,ddne0_dx,vne0,dvne0_dx);
  else
    trap.dne = dne0;
    trap.vne = vne0;
    trap.ddne_dx = ddne0_dx;
    trap.dvne_dx = dvne0_dx;
  end
end

%% full flux calculation
if isfield(params,'testgamma') && params.testgamma
  [trap.Gammae,trap.dGammae_dx,trap.dGammae_dxdot,trap.dGammae_du] = ...
    n_transport(stap,geop,trap,it,model,params);
  trap.dne = zz;
  trap.vne = zz;
  trap.ddne_dx = dzz_dx;
  trap.ddne_dxdot = dzz_dx;
  trap.dvne_dx = dzz_dx;
  trap.dvne_dxdot = dzz_dx;
  trap.dvne_dxu = dzz_du;
else
  trap.Gammae = zz;
  trap.dGammae_dx = dzz_dx;
  trap.dGammae_dxdot = dzz_dx;
  trap.dGammae_du = dzz_du;
end

%% Fast ions
% place holder

%% Neoclassical Tearing Mode Growth rate
if params.ntm.active && strcmp(model.ntm.method,'state')
  if params.ntm.implicit
    error('not supported')
  else % explicit
    if inewt == 0 % calculate once
      [dwntmdt] ...
        = mre_terms(stapk,trap,geop,u,model,params,it);
    end
    trap.dwntmdt = dwntmdt;
  end
end

end

function [chi,dchi_dx,dchi_dxdot,dchi_du] = chie_defaults(chi,dchi_dx,dchi_dxdot,dchi_du)
% pass directly (shorthand for saving space);
end