function [Bv,dBv_dx] = eval_Bv(x,g,v,model)
% function Pthr = eval_Bv(x,g,v,model)
% Vertical magnetic field: Wesson "Tokamaks" chapter "Vacuum magnetic field".
% EDPHY "Plasma control" course by FF (2016): Slides-4_Position control.pdf, slide 12.

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

% Formula: abs(Bv) = mu*Ip*(ln(8*R/(a*kappa^0.5)) + betap + li/2 - 1.5)/(4*pi*R)
% Main radius
R0 = model.equi.R0;
% Inverse aspect ratio
epsilonprof = eval_epsilon(x,g,v,model,false);
epsilon = epsilonprof(end);
% Elongation
kappaprof = eval_kappa(x,g,v,model,false);
kappa = kappaprof(end);
% Beta poloidal
[betap,dbetap_dx] = eval_betap(x,g,v,model);
% Plasma inductance
[li3prof,dli3prof_dx] = eval_li3(x,g,v,model);
li3 = li3prof(end);
dli3_dx = dli3prof_dx(end,:);
% Plasma current
[Ipprof,dIpprof_dx] = eval_Ipl(x,g,v,model,true);
Ip = Ipprof(end);
dIp_dx = dIpprof_dx(end,:);

% Get Bv value
Bv = 1.e-7*Ip*(log(8./(epsilon*sqrt(kappa))) + betap + li3*0.5 - 1.5)/R0;

% Get gradient
%dBv_dx = dIp_dx*Bv/Ip + dbetap_dx*mu*Ip/(4*pi*R) + dli3_dx*mu*Ip/(8*pi*R)
if der
    dBv_dx = dIp_dx.*Bv./Ip + dbetap_dx*1.e-7*Ip/(4*pi*R0) + dli3_dx*1.e-7*Ip./(8*pi*R0);
end
return