function [ti,dti_dx] = eval_ti(x,g,v,model,gauss)
% function [ti,dti_dx] = eval_ti(x,g,v,model,gauss)
% #codegen
use_core_spline_base = hmode_spline_base_active(v,model);

if ~gauss
  dti_dx = zeros(model.rgrid.nrho,numel(x)); % init
else
  dti_dx = zeros(model.rgrid.nrhogauss,numel(x)); % init
end

switch model.ti.method
  case 'state'
    if ~use_core_spline_base
      % ti is directly specified by x
      if gauss
        dti_dtihat = model.ti.Lamgauss;
      else
        dti_dtihat = model.ti.Lam;
      end
      ti = dti_dtihat*x(model.ti.xind,:);
      dti_dx(:,model.ti.xind) = dti_dtihat;
    else
      if size(v,2)>1 
        if any(v(model.hmode.vind.activation, :)) && ~all(v(model.hmode.vind.activation, :))
          error('eval_ti for time dependent V not implemented for varying spline basis')
        end
      end
      if gauss
        dti_dxx = model.ti.Lamgauss_hmode;
      else
        dti_dxx = model.ti.Lam_hmode;
      end
      xx = [x(model.ti.xind,:); v(model.hmode.vind.ti_rhoedge,:)./model.ti.scal];
      ti = dti_dxx*xx;
      dti_dtihat = dti_dxx(:, 1:end-1);
      dti_dx(:,model.ti.xind) = dti_dtihat;
    end
  case 'tescal'
    % ti is scaled version of te.
    % scaling of x coefficients is stored in v
    [te,dte_dx] = eval_te(x,g,v,model,gauss);
    if gauss
      scal = model.ti.Lamgauss*v(model.ti.vind,:);
    else
      scal = model.ti.Lam*v(model.ti.vind,:);
    end
    ti =  scal.*te;
    
    if nargout>1
      dti_dx = bsxfun(@times,scal,dte_dx);
    end
  case 'direct'
    % ti is directly specified by v
    if gauss
      ti = model.ti.Lamgauss*v(model.ti.vind,:);
    else
      ti = model.ti.Lam*v(model.ti.vind,:);
    end
    % dti_dx=0 (already initialized)
  otherwise
    error('not implemented yet')
end

return
