function [nres,dnres_dxdot,dnres_dx,dnres_du,dnres_dd,dnres_dw] ...
  = n_equation(x,xdot,~,geop,stap,trap,BCs,model,species)

use_standard_basis = ~( stap.lamH && strcmp(model.hmode.modeltype,'imposed') );
% flag to signal use of standard vs special hmode spline basis 
if use_standard_basis
  mats = model.mats;
else
  mats = model.mats_hmode;
end

nx = model.dims.nx;
nu = model.dims.nu;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ne or ni part

switch species
  case 'e'
    tt = model.ne;
    if ~strcmp(tt.method,'state')
      [nres,dnres_dxdot,dnres_dx,dnres_du,dnres_dd,dnres_dw] = return_empty(model);
      return
    end
    % transport coefficients
    Ds  = trap.dne;
    Vs  = trap.vne;
    Ss  = trap.sne;
    
    Fs  = trap.Gammae;
    
    % boundary condition
    BCrhs = BCs.neValue;
    dBCrhs_dx = BCs.dneValue_dx;
    dBCrhs_du = BCs.dneValue_du;
    
    BCrho = BCs.neRho;
    
    % Density transport coefficients
    dDs_dx    = trap.ddne_dx;
    dDs_dxdot = trap.ddne_dxdot;
    dDs_du    = trap.ddne_du;
    
    dVs_dx    = trap.dvne_dx;
    dVs_dxdot = trap.dvne_dxdot;
    dVs_du    = trap.dvne_du;
    
    dSs_dx    = trap.dsne_dx;
    dSs_dxdot = trap.dsne_dxdot;
    dSs_du    = trap.dsne_du;
    
    dFs_dx    = trap.dGammae_dx;
    dFs_dxdot = trap.dGammae_dxdot;
    dFs_du    = trap.dGammae_du;
  case 'i'
    tt = model.ni;
    if ~strcmp(tt.method,'state')
      [nres,dnres_dxdot,dnres_dx,dnres_du,dnres_dd,dnres_dw] = return_empty(model);
      return
    end
    % transport coefficients
    Ds  = trap.dni;
    Vs  = trap.vpi;
    Ss  = trap.sni;
    
    Fs  = trap.fni;
    
    % boundary condition
    BCrhs = BCs.niValue;
    BCrho = BCs.niRho;
    
    % Density transport coefficients
    dDs_dx    = trap.ddni_dx;
    dDs_dxdot = trap.ddni_dxdot;
    dDs_du    = trap.ddni_du;
    
    dVs_dx    = trap.dvpi_dx;
    dVs_dxdot = trap.dvpi_dxdot;
    dVs_du    = trap.dvpi_du;
    
    dSs_dx    = trap.dsni_dx;
    dSs_dxdot = trap.dsni_dxdot;
    dSs_du    = trap.dsni_du;
    
    dFs_dx    = trap.dfni_dx;
    dFs_dxdot = trap.dfni_dxdot;
    dFs_du    = trap.dfni_du;
  otherwise
    error('unknown species')
end
%% Matrices
rhogauss = model.rgrid.rhogauss;
Vpgauss = geop.Vp;
dVpdtgauss = geop.dVpdt;

Phib    = geop.Phib;
Phibdot = geop.Phibdot;
g0 = geop.g0;
g1 = geop.g1;

mn_varying = Vpgauss;
an_varying = Phibdot/(2*Phib)*Vpgauss.*rhogauss - Vs.*g0;
dn_varying = g1./Vpgauss.*Ds;
hn_varying = -dVpdtgauss;
sn_varying = Vpgauss.* Ss;
fn_varying = Vpgauss.* Fs;

Mn = zeros(tt.nsp);
Mn(mats.Mn.XIindex) = mats.Mn.XI*mn_varying;

An = zeros(tt.nsp);
An(mats.An.XIindex) = mats.An.XI*an_varying;

Dn = zeros(tt.nsp);
Dn(mats.Dn.XIindex) = mats.Dn.XI*dn_varying;

Hn = zeros(tt.nsp);
Hn(mats.Hn.XIindex) = mats.Hn.XI*hn_varying;

sn = mats.Sn.XI*sn_varying;

fn = model.mats.Fn.XI*fn_varying;

% Dn is rank-deficient since nhat(end) must be externally imposed
% (boundary condition)


% handle boundary conditions
switch tt.BC.type
  case {'Dirichlet','LineAverage','VolumeAverage'}
    switch tt.BC.type % change left hand side of last equation
      case 'Dirichlet'
        if use_standard_basis
          [BClhs,iremove] = get_BC_lhs(model.rgrid.rho,BCrho,tt); % left hand side of BC equation
        else
          % BCs.neRho == -1 indicates the imposed pedestal is active, BC
          % is the amplitude of the outermost core spline
          BClhs = zeros(1,tt.nsp);
          BClhs(tt.nsp) = tt.scal;
          iremove = tt.nsp;
        end
      case 'LineAverage'
        BClhs = sum(bsxfun(@times,model.rgrid.wgauss,tt.Lamgauss));
        iremove = tt.nsp;
      case 'VolumeAverage'
        BClhs = sum(bsxfun(@times,model.rgrid.wgauss.*geop.Vp,tt.Lamgauss));
        iremove = tt.nsp;
      otherwise
        error('unknown BC type');
    end
    
    dscal = 1e22;
    % set to zero last row of the matrices
    Dn(iremove,:) = 0;
    Mn(iremove,:) = 0;
    An(iremove,:) = 0;
    Hn(iremove,:) = 0;
    fn(iremove)   = 0;
    Dn(iremove,:) = dscal * (BClhs/tt.scal); % larger than 1 to be of similar scale as other elements (for conditioning)
    sn(iremove)   = dscal * (BCrhs/tt.scal); % imposed boundary condition
  otherwise
    error('unknown BC option')
end

nhat     = x(tt.xind);
nhatdot  = xdot(tt.xind);

% evaluate ne equation residual
nres = -Mn*nhatdot - An*nhat  - Dn*nhat + Hn*nhat + fn + sn;

%% JACOBIANS
if nargout > 1
  % general structure:
  % f = - M*ndot - A*n - D*n + H*n + s + BC
  % df/dxdot = -dM/dxdot * ndot - M*(dndot/dx)  - dA/dxdot * n           - dD/dxdot * n           + dH/dxdot *n            + ds/dxdot
  % df/dx    = -dM/dx    * ndot                  - dA/dx * n - A*(dn/dx) - dD/dx * n - D*(dn/dx) + dH/dx * n + H*(dn/dx) + ds/dx
  % df/du    = -dM/du    * ndot                  - dA/du * n              - dD/du * n              + dH/du * n              + ds/du
  
  zx = zeros(tt.nsp,nx); % zero matrix of right size
  zu = zeros(tt.nsp,nu); % zero matrix of right size
  % Initialize
  dnhat_dx  = zeros(tt.nsp,nx);
  dnhat_dx(:,tt.xind)   = eye(tt.nsp);
  
  % allocate
  dsdist_dd = [zeros(model.rgrid.nrhogauss,numel(model.dist.psi.ind)),model.dist.te.Lamgauss];
  dsdist_dw = [zeros(size(model.dist.te.Lamgaussw)),model.dist.te.Lamgaussw];
  
  dSs_dd     = dsdist_dd;
  dSs_dw     = dsdist_dw;
  
  % Dn
  ddnvarying_dx    = bsxfun(@times,geop.g1./geop.Vp,dDs_dx);
  ddnvarying_dxdot = bsxfun(@times,geop.g1./geop.Vp,dDs_dxdot);
  ddnvarying_du    = bsxfun(@times,geop.g1./geop.Vp,dDs_du);
  
  % An
  danvarying_dx    = -bsxfun(@times,g0,dVs_dx);
  danvarying_dxdot = -bsxfun(@times,g0,dVs_dxdot);
  danvarying_du    = -bsxfun(@times,g0,dVs_du);
  
  % Sn
  dsnvarying_dxdot = bsxfun(@times,Vpgauss,dSs_dxdot);
  dsnvarying_dx    = bsxfun(@times,Vpgauss,dSs_dx);
  dsnvarying_du    = bsxfun(@times,Vpgauss,dSs_du);
  dsnvarying_dd    = bsxfun(@times,Vpgauss,dSs_dd);
  dsnvarying_dw    = bsxfun(@times,Vpgauss,dSs_dw);
  
  % Fn
  dfnvarying_dxdot = bsxfun(@times,Vpgauss,dFs_dxdot);
  dfnvarying_dx    = bsxfun(@times,Vpgauss,dFs_dx);
  dfnvarying_du    = bsxfun(@times,Vpgauss,dFs_du);
  
  % Jacobians for FE matrices
  
  % Set Mn Jacs to zero: no ne dependence
  dMn_dxdot_nhatdot = zx;
  dMn_dx_nhatdot    = zx;
  dMn_du_nhatdot    = zu;
  
  % Set Hn Jacs to zero: no ne dependence
  dHn_dxdot_nhat = zx;
  dHn_dx_nhat    = zx;
  dHn_du_nhat    = zu;
  
  % D_ne
  dDn_dx_nhat    = matJac(ddnvarying_dx,mats.Dn,nhat,nx,tt);
  dDn_dxdot_nhat = matJac(ddnvarying_dxdot,mats.Dn,nhat,nx,tt);
  dDn_du_nhat    = matJac(ddnvarying_du,mats.Dn,nhat,nu,tt);
  % A_ne
  dAn_dx_nhat    = matJac(danvarying_dx,mats.An,nhat,nx,tt);
  dAn_dxdot_nhat = matJac(danvarying_dxdot,mats.An,nhat,nx,tt);
  dAn_du_nhat    = matJac(danvarying_du,mats.An,nhat,nu,tt);
  
  
  
  %%
  dsn_dx    = mats.Sn.XI*dsnvarying_dx;
  dsn_dxdot = mats.Sn.XI*dsnvarying_dxdot;
  dsn_du    = mats.Sn.XI*dsnvarying_du;
  dsn_dd    = mats.Sn.XI*dsnvarying_dd;
  dsn_dw    = mats.Sn.XI*dsnvarying_dw;
  
  dfn_dx    = model.mats.Fn.XI*dfnvarying_dx;
  dfn_dxdot = model.mats.Fn.XI*dfnvarying_dxdot;
  dfn_du    = model.mats.Fn.XI*dfnvarying_du;
  %%
  
  dsn_dx(iremove,:) = dscal*dBCrhs_dx/tt.scal; % last element set by boundary condition
  dsn_dxdot(iremove,:) = 0; % last element set by boundary condition
  dsn_du(iremove,:) = dscal*dBCrhs_du/tt.scal; % last element set by boundary condition
  dsn_dd(iremove,:) = 0; % last element set by boundary condition
  dsn_dw(iremove,:) = 0; % last element set by boundary condition
  
  dfn_dx(iremove,:) = 0; % last element set by boundary condition
  dfn_dxdot(iremove,:) = 0; % last element set by boundary condition
  dfn_du(iremove,:) = 0; % last element set by boundary condition
  
  dDn_dx_nhat(iremove,:) = 0; % last element set by boundary condition
  dDn_dxdot_nhat(iremove,:) = 0; % last element set by boundary condition
  dDn_du_nhat(iremove,:) = 0; % last element set by boundary condition
  
  dMn_dxdot_nhatdot(iremove,:) = 0;
  dMn_dx_nhatdot(iremove,:) = 0;
  dMn_du_nhatdot(iremove,:) = 0;
  Mn(iremove,:) = 0;
  
  dAn_dxdot_nhat(iremove,:) = 0;
  dAn_dx_nhat(iremove,:) = 0;
  dAn_du_nhat(iremove,:) = 0;
  An(iremove,:) = 0;
  
  dHn_dxdot_nhat(iremove,:) = 0;
  dHn_dx_nhat(iremove,:) = 0;
  dHn_du_nhat(iremove,:) = 0;
  Hn(iremove,:) = 0;
  
  % jacobians for n equation
  dnres_dxdot = -dMn_dxdot_nhatdot - Mn*dnhat_dx - dAn_dxdot_nhat - dDn_dxdot_nhat + dHn_dxdot_nhat + dfn_dxdot + dsn_dxdot;
  dnres_dx    = -dMn_dx_nhatdot    - dAn_dx_nhat - An*dnhat_dx    - Dn*dnhat_dx    - dDn_dx_nhat    + dHn_dx_nhat + Hn*dnhat_dx + dfn_dx + dsn_dx;
  dnres_du    = -dMn_du_nhatdot    - dAn_du_nhat                  - dDn_du_nhat    + dHn_du_nhat    + dfn_du    + dsn_du;
  dnres_dd    =  dsn_dd; 
  dnres_dw    =  dsn_dw;
end

end % end of

function [nres,dnres_dxdot,dnres_dx,dnres_du,dnres_dd,dnres_dw] = return_empty(model)
% return empty matrices of right size
nres = [];
dnres_dxdot = zeros(0,model.dims.nx);
dnres_dx    = zeros(0,model.dims.nx);
dnres_du    = zeros(0,model.dims.nu);
dnres_dd    = zeros(0,model.dims.nd);
dnres_dw    = zeros(0,model.dims.nw);
end
