%[Bx,By,Bz]=CurrentLoop(x,y,z,px,py,pz,fx,fy,fz,r,I) - A. Maurer
%x y z:     loop center location
%px py pz:  loop orientation vector components (loop boresight)
%fx fy fz:  field measurement location
%r:         loop radius
%I:         loop current

function B=CurrentLoop(x,y,z,px,py,pz,fx,fy,fz,r,I)

%one loop frame basis vector is "on axis"
L3=[px py pz];

%choose arbitrary orthogonal axis (since loops are symmetrical)
if pz~=0
    L1=[1 1 (-px-py)/pz];  
elseif py~=0
    L1=[1 (-px-pz)/py 1];   
elseif px~=0
    L1=[(-py-pz)/px 1 1]; 
end
L1=L1/norm(L1);        

%round it out with a right handed set
L2=cross(L1,L3);

%rock out with an inertial frame that is intuitive
I1=[1 0 0];
I2=[0 1 0];
I3=[0 0 1];

%transformation matrices, Loop to Inertial, visa versa
TLI=[dot(I1,L1) dot(I2,L1) dot(I3,L1); dot(I1,L2) dot(I2,L2) dot(I3,L2); dot(I1,L3) dot(I2,L3) dot(I3,L3);];
TIL=transpose(TLI);

%Transform the measurement coordinations into the loop frame, centered on
%loop center
T=[fx,fy,fz]*TIL-[x,y,z]*TIL;

B=[0, 0,0];
%octogon "loops"
n=8
for i=0:1:n-1
    theta=i*pi/(n/2);
    %find element direction
    dl= [-sin(theta)    cos(theta)     0];
    %find element location
    l=  [cos(theta)*r   sin(theta)*r   0];
    %Biot-Savart goodness
    B=B+(1.26E-6*I/(4*pi))*cross(dl,T-l)/(norm(T-l)^3); 
end

%transform the field vector back to the inertial frame
B=B*TLI;