a=1.0;
N=0.75;
vx=0;
vy=0.2;
vz=1;
lamc=1.0;


nspx=301; 
nspy=11;
nspz=61;
nfpx=20;
nfpy=5;
nfpz=4;
xw=20.0*lamc ;
yw=1.0 ;
zw=20.0*lamc ;
bx= zeros(nfpx,nfpy,nfpz);
by= zeros(nfpx,nfpy,nfpz);
bz= zeros(nfpx,nfpy,nfpz);

% focusing on corner where things vary so [-50,y,-250] to some ways in from
% this corner.  Defining [0,0,0] to be x center, y top surface, z front
% e.g., putting the source points from y[-1 to 0], means slab's top surface is
% at y=0.
xs= -50.0            : xw/(nspx-1)  : -50.0 + xw;
ys= -yw              : yw/(nspy-1)  : 0.0;
zs= 0.0              : zw/(nspz-1)  : zw;
%x= -50.0 + xw/nfpx/2 : xw/nfpx      : -50.0 + xw/2 - xw/nfpx/2;
x= -50.0+0.1 : xw/3.0/nfpx : -50.0 + xw/3.0 ;
y= 0.2               : 2.0*yw/(nfpy-1)  : 0.2 + 2.0*yw ;
%z= [0.1 lamc+0.1  2*lamc+0.1  3*lamc+0.1]
z= 0.1 : zw/3.0/nfpz : zw/3.0 ;

size(xs)
size(ys)
size(zs)
size(x)
size(y)
size(z)

%r = [1 2 3 4 5]
%tmp1 = 1/0.1 + 1./r;
%tmp2 = exp(-r/2.0)./r.^2

for i=1:nfpx;
for j=1:nfpy;
for k=1:nfpz;

for is=1:nspx;
for js=1:nspy;
for ks=1:nspz;
r=sqrt((x(i)-xs(is))^2+(y(j)-ys(js))^2+(z(k)-zs(ks))^2);

tmpx = a*N*(vy*(z(k)-zs(ks))-vz*(y(j)-ys(js)))*(1/lamc+1/r)*exp(-r/lamc)/r^2;
bx(i,j,k) = bx(i,j,k) + tmpx*xw*yw*zw/nspx/nspy/nspz;

tmpy = a*N*(vz*(x(i)-xs(is))-vx*(z(k)-zs(ks)))*(1/lamc+1/r)*exp(-r/lamc)/r^2;
by(i,j,k) = by(i,j,k) + tmpy*xw*yw*zw/nspx/nspy/nspz;

tmpz = a*N*(vx*(y(j)-ys(js))-vy*(x(i)-xs(is)))*(1/lamc+1/r)*exp(-r/lamc)/r^2;
bz(i,j,k) = bz(i,j,k) + tmpz*xw*yw*zw/nspx/nspy/nspz;

%if isnan(tmpx)
%     fprintf(' tmpx Nan : %i %i %i %d \n', i,j,k,r)
%end
    
%if abs(tmpz) > 0.001 
%    fprintf(' not small tmpz : %d \n', tmpz)
%end
%if abs(bz) > 0.001 
%    fprintf(' not small Bz : %i %i %i %d \n', i,j,k,bz(i,j,k))
%end

end
end
end
end
end
end


bxx22 = bx(:,2,2);
byx22 = by(:,2,2);
bxx42 = bx(:,4,2);
byx42 = by(:,4,2);
bx1 = squeeze(bx(:,1,:));
bx2 = squeeze(bx(:,2,:));
bx3 = squeeze(bx(:,3,:));
bx4 = squeeze(bx(:,4,:));
bx5 = squeeze(bx(:,5,:));



figure
[yyy,xxx,zzz]=meshgrid(y,x,z);
quiver3(zzz,xxx,yyy,bz,bx,by)
title(['B-vec(x,y,z)  lamc=',num2str(lamc),'mm'])
xlabel('z (mm)')
ylabel('x (mm)')
zlabel('y (mm)')
axis([0 zw/3.0 -50 -50+xw/3.0  -0.5 2.5])
view(-50,10)

figure
subplot(2,2,1)
plot(x,bx(:,2,2))
title(['Bx(x,',num2str(y(2)),',',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('Bx')
subplot(2,2,2)
plot(x,by(:,2,2))
title(['By(x,',num2str(y(2)),',',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('By')
subplot(2,2,3)
plot(x,bz(:,2,2))
title(['Bz(x,',num2str(y(2)),',',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('Bz')

bxxy2 = squeeze(bx(:,:,2));
byxy2 = squeeze(by(:,:,2));
bzxy2 = squeeze(bz(:,:,2));
[yy,xx] = meshgrid(y,x);

figure
subplot(2,2,1)
quiver(xx,yy,bxxy2,byxy2,0.3,'MaxHeadSize',0.1,'AlignVertexCenters','on')
title(['Bxy-vec(x,y,',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

subplot(2,2,2)
contourf(xx,yy,bxxy2)
title(['Bx(x,y,',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

subplot(2,2,3)
contourf(xx,yy,byxy2)
title(['By(x,y,',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

subplot(2,2,4)
contourf(xx,yy,bzxy2)
title(['Bz(x,y,',num2str(z(2)),')  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

%save fielddata.mat
%save xfieldout bx -ASCII