function [Teres,dTeres_dxdot,dTeres_dx,dTeres_du,dTeres_dd,dTeres_dw] ...
    = te_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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Te part

switch species
  case 'e'
    tt = model.te;
    if ~strcmp(tt.method,'state')
      [Teres,dTeres_dxdot,dTeres_dx,dTeres_du,dTeres_dd,dTeres_dw] = return_empty(model);
      return
    end
    % Total power density
    power   = trap.petot;
    dens    = stap.ne;
    ddens_dx = stap.dne_dx;
    
    densdot = stap.nedot;
    ddensdot_dxdot = stap.dnedot_dxdot;
    
    chi     = trap.chie;
    
    BCrhs = BCs.teValue;
    dBCrhs_dx = BCs.dteValue_dx;
    dBCrhs_du = BCs.dteValue_du;
    BCrho = BCs.teRho;
    
    % Pe source derivatives
    dPower_dxdot  = trap.dpetot_dxdot;
    dPower_dx     = trap.dpetot_dx;
    dPower_du     = trap.dpetot_du;
    dchi_dx       = trap.dchie_dx;
    dchi_dxdot    = trap.dchie_dxdot;
    dchi_du       = trap.dchie_du;
    
  case 'i'
    tt = model.ti; 
    if ~strcmp(tt.method,'state')
      [Teres,dTeres_dxdot,dTeres_dx,dTeres_du,dTeres_dd,dTeres_dw] = return_empty(model);
      return
    end
    power    = trap.pitot;
    dens     = stap.ni;
    ddens_dx = stap.dni_dx;
    
    ddensdot_dxdot = stap.dnidot_dxdot;
    densdot     = stap.nidot;
    chi         = trap.chii;
    
    BCrhs     = BCs.tiValue;
    dBCrhs_dx = BCs.dtiValue_dx;
    dBCrhs_du = BCs.dtiValue_du;
    BCrho     = BCs.tiRho;
    
    % Pe source derivatives
    dPower_dxdot  = trap.dpitot_dxdot;
    dPower_dx     = trap.dpitot_dx;
    dPower_du     = trap.dpitot_du;
    dchi_dx       = trap.dchii_dx;
    dchi_dxdot    = trap.dchii_dxdot;
    dchi_du       = trap.dchii_du;
    
  otherwise
    error('unknown species')
end
%% Matrices

Vpgauss = geop.Vp;
dVpdtgauss = geop.dVpdt;
dVpdrhohatgauss = geop.Vpp;

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

mTe_varying = Vpgauss.* dens;
aTe_varying = dens.* Vpgauss.* Phibdot/Phib;
dTe_varying = g1.* dens.* chi./Vpgauss;
hTe_varying = -5. * dens.* dVpdtgauss - 3.*Vpgauss.*densdot + ...
  Phibdot/Phib * model.rgrid.rhogauss.* dens .* dVpdrhohatgauss;
sTe_varying = Vpgauss.* power;

% MTe
MTe = zeros(tt.nsp);
MTe(mats.MTe.XIindex) =  mats.MTe.XI*mTe_varying;

% ATe
ATe = zeros(tt.nsp);
ATe(mats.ATe.XIindex) = mats.ATe.XI*aTe_varying;

% DTe
DTe = zeros(tt.nsp);
DTe(mats.DTe.XIindex) = mats.DTe.XI*dTe_varying;

% HTe
HTe = zeros(tt.nsp);
HTe(mats.HTe.XIindex) = mats.HTe.XI*hTe_varying;

% sTe
sTe = mats.STe.XI*sTe_varying;

%lTe =
% DTe is rank-deficient since Tehat(end) must be externally imposed
% (boundary condition)


% handle boundary conditions
switch tt.BC.type
  case {'Dirichlet'}
    % set to zero last row of the matrices
    dscal = 1e6; % typical D value for scaling
    if use_standard_basis
      [BClhs,iremove] = get_BC_lhs(model.rgrid.rho,BCrho,tt); % left hand side of BC equation
    else
      % 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
    
    DTe(iremove,:) = 0;
    MTe(iremove,:) = 0;
    ATe(iremove,:) = 0;
    HTe(iremove,:) = 0;
    
    DTe(iremove,:) = dscal*(BClhs/tt.scal); % scale to size of other elements (for conditioning)
    sTe(iremove)   = dscal*(BCrhs/tt.scal); % imposed boundary condition
    
  otherwise
    error('other BC options presently not supported')
end

Tehat     = x(tt.xind);
Tehatdot  = xdot(tt.xind);

% evaluate Te equation residual ("G" in FF thesis eq7.40)
Teres = -MTe*Tehatdot - ATe*Tehat  - DTe*Tehat + HTe*Tehat   + sTe(1:tt.nsp);

%% JACOBIANS
if nargout > 1
  % general structure:
  % f = - M*Tedot + A*Te - D*Te + H*Te + s + BC
  % df/dxdot = -dM/dxdot * Tedot - M*(dTedot/dx)  - dA/dxdot * Te           - dD/dxdot * Te           + dH/dxdot *Te            + ds/dxdot
  % df/dx    = -dM/dx    * Tedot                  - dA/dx * Te - A*(dTe/dx) - dD/dx * Te - D*(dTe/dx) + dH/dx * Te + H*(dTe/dx) + ds/dx
  % df/du    = -dM/du    * Tedot                  - dA/du * Te              - dD/du * Te              + dH/du * Te              + ds/du
  
  zx = zeros(tt.nsp,nx); % zero matrix of right size
  zu = zeros(tt.nsp,nu); % zero matrix of right size
  %% Te equation
  dTehat_dx  = zeros(tt.nsp,nx);
  dTehat_dx(:,tt.xind)   = eye(tt.nsp);
  
  % allocate
  dpdist_dd = [zeros(model.rgrid.nrhogauss,numel(model.dist.psi.ind)),model.dist.te.Lamgauss];
  dpdist_dw = [zeros(size(model.dist.te.Lamgaussw)),model.dist.te.Lamgaussw];
  
  dPower_dd     = dpdist_dd;
  dPower_dw     = dpdist_dw;
  
  % FI matrices
  % M_Te
  dMTe_dxdot_Tehatdot = zx; % no state_dot dependence
  dMTe_du_Tehatdot    = zu; % no U dependence
  % A_Te
  dATe_dxdot_Tehat = zx; % no state_dot dependence
  dATe_du_Tehat    = zu; % no U dependence
  % H_Te
  dHTe_du_Tehat = zu; % no U dependence
  
  % D_te
  ddTevarying_dx    = bsxfun(@times,dens.*geop.g1./geop.Vp,dchi_dx);
  ddTevarying_dxdot = bsxfun(@times,dens.*geop.g1./geop.Vp,dchi_dxdot);
  ddTevarying_du    = bsxfun(@times,dens.*geop.g1./geop.Vp,dchi_du);
  % S_Te
  dsTevarying_dxdot = bsxfun(@times,Vpgauss,dPower_dxdot);
  dsTevarying_dx    = bsxfun(@times,Vpgauss,dPower_dx);
  dsTevarying_du    = bsxfun(@times,Vpgauss,dPower_du);
  dsTevarying_dd    = bsxfun(@times,Vpgauss,dPower_dd);
  dsTevarying_dw    = bsxfun(@times,Vpgauss,dPower_dw);
  
  if any(any(ddens_dx~=0)) || any(any(ddensdot_dxdot))
    % D_Te
    ddTevarying_dx = ddTevarying_dx + bsxfun(@times,chi.*geop.g1./geop.Vp,ddens_dx);
    % M_Te
    dmTevarying_dx    =  bsxfun(@times,geop.Vp,ddens_dx);
    % A_Te
    daTevarying_dx    = bsxfun(@times,Vpgauss.* Phibdot/Phib,ddens_dx);
    % H_Te
    dhTevarying_dx    = bsxfun(@times,-5.*ddens_dx,dVpdtgauss) + ...
      Phibdot/Phib * bsxfun(@times,model.rgrid.rhogauss.* dVpdrhohatgauss,ddens_dx);
    dhTevarying_dxdot = -bsxfun(@times,3.*Vpgauss,ddensdot_dxdot);
    
    % Jacobians for FE matrices
    dMTe_dx_Tehatdot = matJac(dmTevarying_dx,mats.MTe,Tehatdot,nx,tt);
    dATe_dx_Tehat    = matJac(daTevarying_dx,mats.ATe,Tehat,nx,tt);
    dHTe_dxdot_Tehat = matJac(dhTevarying_dxdot,mats.HTe,Tehat,nx,tt);
    dHTe_dx_Tehat    = matJac(dhTevarying_dx,mats.HTe,Tehat,nx,tt);
  else
    % Set Jacs to zero: no te dependence
    dMTe_dx_Tehatdot = zx;
    dATe_dx_Tehat    = zx;
    dHTe_dxdot_Tehat = zx;
    dHTe_dx_Tehat    = zx;
  end
  
  % D_Te
  dDTe_dx_Tehat    = matJac(ddTevarying_dx,mats.DTe,Tehat,nx,tt);
  dDTe_dxdot_Tehat = matJac(ddTevarying_dxdot,mats.DTe,Tehat,nx,tt);
  dDTe_du_Tehat    = matJac(ddTevarying_du,mats.DTe,Tehat,nu,tt);
  
  %%
  dsTe_dx    = mats.STe.XI*dsTevarying_dx;
  dsTe_dxdot = mats.STe.XI*dsTevarying_dxdot;
  dsTe_du    = mats.STe.XI*dsTevarying_du;
  dsTe_dd    = mats.STe.XI*dsTevarying_dd;
  dsTe_dw    = mats.STe.XI*dsTevarying_dw;
  
  if strcmp(tt.BC.type,'Dirichlet')
    dsTe_dx(iremove,:) = dscal*dBCrhs_dx/tt.scal; % last element set by boundary condition
    dsTe_dxdot(iremove,:) = 0; % last element set by boundary condition
    dsTe_du(iremove,:) = dscal*dBCrhs_du/tt.scal; % last element set by boundary condition
    dsTe_dd(iremove,:) = 0; % last element set by boundary condition
    dsTe_dw(iremove,:) = 0; % last element set by boundary condition
    
    dDTe_dx_Tehat(iremove,:) = 0; % last element set by boundary condition
    dDTe_dxdot_Tehat(iremove,:) = 0; % last element set by boundary condition
    dDTe_du_Tehat(iremove,:) = 0; % last element set by boundary condition
    
    dMTe_dxdot_Tehatdot(iremove,:) = 0;
    dMTe_dx_Tehatdot(iremove,:) = 0;
    dMTe_du_Tehatdot(iremove,:) = 0;
    MTe(iremove,:) = 0;
    
    dATe_dxdot_Tehat(iremove,:) = 0;
    dATe_dx_Tehat(iremove,:) = 0;
    dATe_du_Tehat(iremove,:) = 0;
    ATe(iremove,:) = 0;
    
    dHTe_dxdot_Tehat(iremove,:) = 0;
    dHTe_dx_Tehat(iremove,:) = 0;
    dHTe_du_Tehat(iremove,:) = 0;
    HTe(iremove,:) = 0;
  end
  
  % jacobians for Te equation
  dTeres_dxdot = -dMTe_dxdot_Tehatdot - MTe*dTehat_dx - dATe_dxdot_Tehat - dDTe_dxdot_Tehat + dHTe_dxdot_Tehat  + dsTe_dxdot;
  dTeres_dx    = -dMTe_dx_Tehatdot    - dATe_dx_Tehat - ATe*dTehat_dx    - DTe*dTehat_dx    - dDTe_dx_Tehat     + dHTe_dx_Tehat +  HTe*dTehat_dx   + dsTe_dx;
  dTeres_du    = -dMTe_du_Tehatdot    - dATe_du_Tehat - dDTe_du_Tehat    + dHTe_du_Tehat    + dsTe_du;
  dTeres_dd    =  dsTe_dd;
  dTeres_dw    =  dsTe_dw;
end
end

function [Teres,dTeres_dxdot,dTeres_dx,dTeres_du,dTeres_dd,dTeres_dw] = return_empty(model)
% return empty matrices of right size
Teres = [];
dTeres_dxdot = zeros(0,model.dims.nx);
dTeres_dx    = zeros(0,model.dims.nx);
dTeres_du    = zeros(0,model.dims.nu);
dTeres_dd    = zeros(0,model.dims.nd);
dTeres_dw    = zeros(0,model.dims.nw);
end