
gA=1E-8;  % assuming gA^2=1E-16, which should give typical Bfields of pT
Ngl=0.75;
Ncu=1.9;
vx=0.0*6.6E5;   % 660m/s=6.6E5 mm/s for 6Angstrom neutron
vy=0.0*6.6E5;
vz=1.0*6.6E5;
lamc=1.0;   % units of mm

alpha = 0.004564;   % 4564. s T/m^2  =0.0287 sT/mm^2

nspx=101; 
nspy=101;
nspz=51;
nfpx=20;
nfpy=5;
nfpz=4;
xw=20.0*lamc ;
yw=10.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= -25.0            : xw/(nspx-1)  : -25.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= -25.0+0.1 : xw/3.0/nfpx : -25.0 + xw/3.0 ;
y= 0.2               : 2.0/(nfpy-1)  : 0.2 + 2.0 ;
%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 = alpha*gA^2*Ncu*(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*1E12;  
%  gives B in units of picoT

tmpy = alpha*gA^2*Ncu*(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*1E12;  
% gives B in units of picoT

tmpz = alpha*gA^2*Ncu*(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*1E12; 
% gives B in units of picoT

%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)(pT) lamc=',num2str(lamc),'mm'])
xlabel('z (mm)')
ylabel('x (mm)')
zlabel('y (mm)')
axis([0 zw/3.0 -26 -26+xw/3.0  -0.5 2.5])
view(-50,10)

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

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)),')(pT) 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)),')(pT) 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)),')(pT) 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)),')(pT) lamc=',num2str(lamc),'mm'])
xlabel('x (mm)')
ylabel('y (mm)')

%save fielddata.mat
%save xfieldout bx -ASCII
