function [hR, hZ, flux] = hRZ_espira(R, Zb, Rw, Zw)
% [hR, hZ, phi] =HRZ_ESPIRA(R, Z, Rp)
% Calculo do campo H poloidal e do fluxo  provocado por uma espira na coordeenada R=Rw, Z=Zw
% nos pontos (R, Z) , I = 1 Amp
% R  -- Array- Coordenadas R dos pontos 
% Z  -- Array- Coordenadas Z dos pontos 
% Rw -- escalar- Raio da espira 
% Zw -- escalar- Coordenada Z da espira
if nargin < 4
    Zw=0;
end
if nargin < 3
    Rw=1.;
end
Z = Zb - Zw;
ss = (Rw + R).^2 + Z.^2; 		% (Rw+R)^2  + Z^2
sqrt_ss = sqrt(ss);
num_r = Rw^2 + R.^2 + Z.^2;		% Rw^2 +R^2 + Z^2
num_z = 2 * Rw^2 - num_r; 		% Rw^2 -R^2 - Z^2
%num_z = R.^2 - Rw^2 - Z.^2; 	% Rw^2 -R^2 - Z^2
den= ss - 4*Rw * R;             % (Rw-R)^2 + Z^2
m = 4 * Rw * R ./ ss;
[K, E] = ellipke(m);
flux = sqrt_ss  .* ((1- m/2) .* K -E) ;
hR  = Z  ./ R .* (num_r ./ den .* E - K) ./ sqrt_ss  / 2 / pi;
hZ = (K + num_z ./ den .* E ) ./ sqrt_ss  / 2 / pi ;
%disp(['num_r = ' num2str(num_r)]);
%disp(['num_z = ' num2str(num_z)]);
%disp(['   R = ' num2str(R)]);
%disp(['   Z = ' num2str(Zb)]);
%disp(['  hR = ' num2str(hR)]);
%disp(['  hZ = ' num2str(hZ)]);