function [Bx,By,Bz] = bcoil(coil,h,I,nop)
% BCOIL
%  [Bx,By,Bz] = bcoil(c,h,I,'nop') computes the B-field produced by a current I
%  flowing through the conductor described in c a distance h above
%  the plane of the conductor.  Since the conductor is always drawn
%  on a normalized axis set, the calculation plane is simply evaluated
%  on [-1,1] in both the X and Y directions.
%  Bx,By,Bz are the x,y,z components of the field.
%   The optional argument 'nop' is used to suppress the plot.
%
%  SEE ALSO:  INCOIL, OUTCOIL, RECOIL
%                                                       RJM 4/10/95
%


mu=1;
[x,y]=meshdom(-1:.1:1,-1:.1:1);
[v,w]=size(x);
z = h*ones(v,w);
[n,m]=size(coil);
Bx = zeros(v,w);
By = zeros(v,w);
Bz = zeros(v,w);


for i=1:n,

 xc=coil(i,1);, yc=coil(i,2);, l = coil(i,3);, thc = coil(i,4);


  r = ((x-xc).^2 + (y-yc).^2 + z.^2).^.5;

rx = x-xc;
ry = y-yc;   %  Vector from midpoint to obs point
rz = z;

cox = l/2*cos(thc);
coy = l/2*sin(thc);  % vector from midpoint to end

% Calculate dot product to find angle between vectors

r_dot_co = cox*rx + coy*ry;  % rˇco = ŚrŚŚcoŚcos(phi)
phi = acos(2/l*r_dot_co./r);


%Calculate cross product and normalize
% to find direction of field

B_norm = ((coy.*rz).^2 + (cox.*rz).^2 + (cox.*ry-rx.*coy).^2).^.5;
B_unitx = coy.*rz./B_norm;
B_unity = -cox.*rx./B_norm;
B_unitz = (cox.*ry-rx.*coy)./B_norm;


  Bmag = mu/4/pi*I*l*sin(phi)./(r.^2);

%Calculate each field component

  Bx = Bx + Bmag*B_unitx;
  By = By + Bmag*B_unity;
  Bz = Bz + Bmag*B_unitz;

end

if nargin ~= 4,
clg
subplot(2,2,1),mesh(x,y,Bz),title('Z-component'),shg
subplot(2,2,2),mesh(x,y,Bx),title('X-component')
subplot(2,2,3),mesh(x,y,By),title('Y-component')
subplot(2,2,4),hold on,outcoil(coil),xlabel('X'),ylabel('Y'),hold off

end