function [Pthr] = eval_PLHthr(x,g,v,model)
% function Pthr = eval_PLHthr(x,g,v,model)
% H-mode Pthreshold scaling from Y.Martin JP Conf.series 123 (2008)
% #codegen

ne = eval_neLHL(x,g,v,model,true);
ne20lav = sum(bsxfun(@times,model.rgrid.wgauss,ne))*1e-20; % int ne drhotor
B0 = model.equi.B0;
R0 = model.equi.R0;
a  = R0*model.equi.epsilon;
% ne20lav : line averaged density in units of 10^20. 
%     a : minor radius 
%    B0 : toroidal field
%
%   Pthr: LH threshold power in [W]. All inputs in SI

ne20lavmin = 0.4;
Pthr = 2.15e+6 * exp(1)^0.107 * ne20lav.^0.782 * B0^0.772 * a^0.975 * R0^0.999;
Pthr = Pthr.*bsxfun(@max,1,(ne20lavmin./ne20lav).^0.782.*exp(-20*(ne20lav-ne20lavmin)));

return