function [Wthi,dWthi_dx] = eval_Wthi(x,g,v,model,gaussgrid)
% function [Wth,dWth_dx] = eval_Wthi(x,g,v,model,gaussgrid)
%#codegen

if gaussgrid
    error('Wthi is an integral quantity, can not be evaluated on gauss grid')
end

if nargout>1
    der = true; % compute derivative
    dWthi_dx = zeros(1,numel(x)); % init
else
    der = false;
end

 [ti,dti_dx] = eval_ti(x,g,v,model,true);
 [ni,dni_dx] = eval_ni(x,g,v,model,true);
 [n1,dn1_dx] = eval_nx(1,x,g,v,model,true,false);
 [n2,dn2_dx] = eval_nx(2,x,g,v,model,true,false);
 [n3,dn3_dx] = eval_nx(3,x,g,v,model,true,false);

geop = geometry_profiles(g,0*g,model);

nitot = ni + n1 + n2 + n3; % all thermal species
dnitot_dx = dni_dx + dn1_dx + dn2_dx + dn3_dx; 

Wthi = int_Vtot(3/2*1.6e-19*(ti.*nitot),geop,model); % total plasma thermal energy
if der
    dWthi_dx = int_Vtot(3/2*1.6e-19*(...
        + bsxfun(@times,nitot,dti_dx) ...
        + bsxfun(@times,ti,dnitot_dx)),geop,model);
end

return
