function [ pprime] = eval_pprime(x,~, v, model,gauss)
% #codegen
% function [ pprime] = eval_pprime(x, g, v, model)
% This function calculates  p' = d(pressure)/dpsi
% for integration with equilibrium code
% assumes p = qe*(ni*Ti+ne*Te) 
%
% Outputs:
% pp: dp/dpsi

%% evaluation of in-code quantities
Te = eval_te(x,[],v,model,gauss);
Ti = eval_ti(x,[],v,model,gauss);
ne = eval_ne(x,[],v,model,gauss);
ni = eval_ni(x,[],v,model,gauss);
Tep = eval_tep(x,[],v,model,gauss);
Tip = eval_tip(x,[],v,model,gauss);
nep = eval_nep(x,[],v,model,gauss);
nip = eval_nip(x,[],v,model,gauss);
psip = eval_psip(x,[],v,model,gauss); 

dpdrho = 1.602e-19*(nep.*Te+ nip.*Ti + ne.*Tep + ni.*Tip);
pprime = dpdrho./psip; % pp = dp/dpsi

if ~gauss
    % Spline de l'Hospital
    d2pdrho2_0 = d2pdrho2_0_Hopital(x,v,model, Te(1), ne(1), Ti(1) , ni(1));
    psipp = eval_psipp(x,[],v,model,gauss);
    pprime(1) = d2pdrho2_0/psipp(1); % numerical l'Hopital rule
end

debug_plot = false;
if debug_plot
eval_debugging_plot(mfilename,'psiN',x,model,gauss, Te, 'Te', Ti, 'Ti', ne, 'ne', ni, 'ni', Tep,...
    'Tep', Tip, 'Tip', nep, 'nep', nip, 'nip', psip, 'psip', dpdrho, 'dpdrho', pprime, 'pprime');
end

return

function d2pdrho2_0 = d2pdrho2_0_Hopital(x,v,model, Te0, ne0, Ti0 , ni0)

    % Tep, Tepp
    if strcmp(model.te.method,'state')
        tehat = x(model.te.xind);
    else
        tehat = v(model.te.vind);
    end
    tep0 = model.te.Lamp(1,:)*tehat;
    tepp0 = model.te.Lampp(1,:)*tehat;
    
    % Tip, Tipp
    switch model.ti.method
        case 'state'
            tihat = x(model.ti.xind);
            tip0 = model.ti.Lamp(1,:)*tihat;
            tipp0 = model.ti.Lampp(1,:)*tihat;
        case 'tescal'
            % ti is scaled version of ne.
            % scaling of x coefficients is stored in v
            scal = model.ti.Lam(1,:)*v(model.ti.vind,:);
            tipp0 =  scal.*tepp0;
            tip0 = scal.*tep0;
        case 'direct'
            tipp0 = model.ti.Lampp(1,:)*v(model.ti.vind,:);
            tip0 = model.ti.Lamp(1,:)*v(model.ti.vind,:);
        otherwise
            error('not supported yet')
    end
       
    % nep, nepp
    if strcmp(model.ne.method,'state')
        nehat = x(model.ne.xind);
    else
        nehat = v(model.ne.vind);
    end
    nep0 = model.ne.Lamp(1,:)*nehat;
    nepp0 = model.ne.Lampp(1,:)*nehat;
    
    %nip, nipp
    switch model.ni.method
     case 'nescal'
        % ni is scaled version of ne.
        % scaling of x coefficients is stored in v
        scal = model.ni.Lam(1,:)*v(model.ni.vind,:);
        nipp0 =  scal.*nepp0;
        nip0 = scal.*nep0;

    case 'direct'
        % ni is directly specified by v
            nip0 = model.ni.Lamp(1,:)*v(model.ni.vind,:);
            nipp0 = model.ni.Lampp(1,:)*v(model.ni.vind,:);
    case 'qnze'
           nip = eval_qnze('ni',x,[],v,model,false,1);
           nipp = eval_qnze('ni',x,[],v,model,false,2);
           nip0 = nip(1); nipp0 = nipp(1);
       
    otherwise
        error('not implemented yet')
    end
    
    d2pdrho2_0 = 1.602e-19*(tepp0*ne0+nepp0*Te0 + 2*tep0*nep0...
         + tipp0*ni0 + nipp0*Ti0 + 2*tip0*nip0);
    
    
return

