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


nspx=21;
nspy=11;
nspz=31;
nfpx=21;
nfpy=6;
nfpz=31;
xw=50.0 ;
yw=1.0 ;
zw=100.0 ;
bx= zeros(nfpx,nfpy,nfpz);
by= zeros(nfpx,nfpy,nfpz);
bz= zeros(nfpx,nfpy,nfpz);

xs= -xw/2 : xw/(nspx-1) :xw/2;
x= -xw/2 : xw/(nfpx-1) :xw/2;
ys= -yw/2 : yw/(nspy-1) :yw/2;
y= 0 : 2.0*yw/(nfpy-1) :2.0*yw;
zs= -zw/2 : zw/(nspz-1) :zw/2;
z= -zw/2 : zw/(nfpz-1) :zw/2;

%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 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,:));


[zz,xx] = meshgrid(zz,xx);
figure
subplot(2,2,1)
surf(xx,zz,bx1)
title(['Bx(x,1,z)   lamc=',num2str(lamc),'mm'])
subplot(2,2,2)
surf(xx,zz,bx3)
title(['Bx(x,3,z)  lamc=',num2str(lamc),'mm'])
subplot(2,2,3)
surf(xx,zz,bx5)
title(['Bx(x,5,z)  lamc=',num2str(lamc),'mm'])


bxxy10 = squeeze(bx(:,:,10));
byxy10 = squeeze(by(:,:,10));
[yy,xx] = meshgrid(y,x);

figure
subplot(2,2,1)
[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([-25 25 -30 30 -0.5 2.5])
view(-50,10)

subplot(2,2,2)
quiver(xx,yy,bxxy10,byxy10,0.3,'MaxHeadSize',0.01,'AlignVertexCenters','on')
title(['Bxy-vec(x,y,10)  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

subplot(2,2,3)
contourf(xx,yy,bxxy10)
title(['Bx(x,y,10)  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')
subplot(2,2,4)
contourf(xx,yy,byxy10)
title(['By(x,y,10)  lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

%save fielddata.mat
%save xfieldout bx -ASCII