function [tep,dtep_dx] = eval_tep(x,~,v,model,gauss)
% function [tep,dtep_dx] = eval_tep(x,g,model,whichgrid)
% #codegen
if size(v,2)>1 && strcmp(model.hmode.modeltype,'imposed')
  if any(v(model.hmode.vind.activation, :)) && ~all(v(model.hmode.vind.activation, :))
    error('eval_tep for time dependent V not implemented for varying spline basis')
  end
end
use_core_spline_base = hmode_spline_base_active(v,model);
if ~use_core_spline_base
  if ~gauss
    dtep_dx = zeros(model.rgrid.nrho,numel(x));
  else
    dtep_dx = zeros(model.rgrid.nrhogauss,numel(x));
  end
  
  switch model.te.method
    case 'state'
      % te is directly specified by v
      if gauss
        dtep_dtehat = model.te.Lampgauss;
      else
        dtep_dtehat = model.te.Lamp;
      end
      tep = dtep_dtehat*x(model.te.xind,:);
      if nargout>1
        dtep_dx(:,model.te.xind) = dtep_dtehat;
      end
    otherwise
      error('model.te.method = %s not implemented',model.te.method)
  end
else
  % this is the case for special H-mode splines
  % init
  if ~gauss;
    dtep_dx = zeros(model.rgrid.nrho,numel(x));
  else
    dtep_dx = zeros(model.rgrid.nrhogauss,numel(x));
  end
  switch model.te.method
    case 'state'
      % te is specified by x and Te(rho=1)
      if size(v,2)>1 
        if any(v(model.hmode.vind.activation, :)) && ~all(v(model.hmode.vind.activation, :))
          error('eval_tep for time dependent V not implemented for varying spline basis')
        end
      end
      if gauss
        dtep_dxx = model.te.Lampgauss_hmode;
        dtep_dtehat = model.te.Lampgauss_hmode(:, 1:end-1);
      else
        dtep_dxx = model.te.Lamp_hmode;
        dtep_dtehat = model.te.Lamp_hmode(:, 1:end-1);
      end
      xx = [x(model.te.xind,:); v(model.hmode.vind.te_rhoedge,:)./model.te.scal];
      tep = dtep_dxx*xx;
    otherwise
      error('model.te.method = %s not implemented',model.te.method)
  end
  
  if nargout>1
    dtep_dx(:,model.te.xind) = dtep_dtehat;
  end
end

return
