function [V, H]=quadmag(x,z,Jx,Jz)

% [V, H]=quadmag(x,z,Jx,Jz)
% computes the vertical and horizontal magnetic field strengths
% due to a buried quadrilateral element whose vertices NUMBERED CLOCKWISE
% have coordinates x=[x1 x2 x3 x4] and z=[z1 z2 z3 z4]
% and that has a magnetization vector J = [Jx Jz]
% V and H have units of J.
% J.-A. Olive, April 2010, after Talwani and Heirtzler, 1964
%
% modified by J.-A. Olive, October 2015: replacing phase.m by angle.m
% for more recent versions of Matlab


i=sqrt(-1);

% Vertex A
x1=x(1); x2=x(2);
z1=z(1); z2=z(2);
a1=x1+i*z1; t1=angle(a1);
a2=x2+i*z2; t2=angle(a2);
x12=x1-x2;
z21=z2-z1;
r1=sqrt(x1^2+z1^2);
r2=sqrt(x2^2+z2^2);
P=z21^2*(t1-t2)/(z21^2+x12^2)+z21*x12*log(r2/r1)/(z21^2+x12^2);
Q=z21*x12*(t1-t2)/(z21^2+x12^2)+z21^2*log(r2/r1)/(z21^2+x12^2);
Va=2*(Jx*Q-Jz*P);
Ha=2*(Jx*P+Jz*Q);

% Vertex B
x1=x(2); x2=x(3);
z1=z(2); z2=z(3);
a1=x1+i*z1; t1=angle(a1);
a2=x2+i*z2; t2=angle(a2);
x12=x1-x2;
z21=z2-z1;
r1=sqrt(x1^2+z1^2);
r2=sqrt(x2^2+z2^2);
P=z21^2*(t1-t2)/(z21^2+x12^2)+z21*x12*log(r2/r1)/(z21^2+x12^2);
Q=z21*x12*(t1-t2)/(z21^2+x12^2)+z21^2*log(r2/r1)/(z21^2+x12^2);
Vb=2*(Jx*Q-Jz*P);
Hb=2*(Jx*P+Jz*Q);

%Vertex C
x1=x(3); x2=x(4);
z1=z(3); z2=z(4);
a1=x1+i*z1; t1=angle(a1);
a2=x2+i*z2; t2=angle(a2);
x12=x1-x2;
z21=z2-z1;
r1=sqrt(x1^2+z1^2);
r2=sqrt(x2^2+z2^2);
P=z21^2*(t1-t2)/(z21^2+x12^2)+z21*x12*log(r2/r1)/(z21^2+x12^2);
Q=z21*x12*(t1-t2)/(z21^2+x12^2)+z21^2*log(r2/r1)/(z21^2+x12^2);
Vc=2*(Jx*Q-Jz*P);
Hc=2*(Jx*P+Jz*Q);

%Vertex D
x1=x(4); x2=x(1);
z1=z(4); z2=z(1);
a1=x1+i*z1; t1=angle(a1);
a2=x2+i*z2; t2=angle(a2);
x12=x1-x2;
z21=z2-z1;
r1=sqrt(x1^2+z1^2);
r2=sqrt(x2^2+z2^2);
P=z21^2*(t1-t2)/(z21^2+x12^2)+z21*x12*log(r2/r1)/(z21^2+x12^2);
Q=z21*x12*(t1-t2)/(z21^2+x12^2)+z21^2*log(r2/r1)/(z21^2+x12^2);
Vd=2*(Jx*Q-Jz*P);
Hd=2*(Jx*P+Jz*Q);

V=Va+Vb+Vc+Vd;
H=Ha+Hb+Hc+Hd;

