function [Wthe,dWthe_dx] = eval_Wthe(x,g,v,model,gaussgrid)
% function [Wth,dWth_dx] = eval_Wthe(x,g,v,model,gaussgrid)
%#codegen

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

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

 [te,dte_dx] = eval_te(x,g,v,model,true);
 [ne,dne_dx] = eval_ne(x,g,v,model,true);

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

Wthe = int_Vtot(3/2*1.6e-19*(te.*ne),geop,model); % total plasma thermal energy
if der
    dWthe_dx = int_Vtot(3/2*1.6e-19*(...
          bsxfun(@times,ne,dte_dx) ...
        + bsxfun(@times,te,dne_dx)) ,geop,model);
end

return
