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=201;   % 5 per lambda givse accruate value for lamc=1
nspy=31;
nspz=201;
nfpx=20;
nfpy=6;
nfpz=4;
xw=20.0*lamc ;
yw=1.0 ;
zw=20.0*lamc ;
bx1= zeros(nfpx,nfpy,nfpz);
by1= zeros(nfpx,nfpy,nfpz);
bz1= zeros(nfpx,nfpy,nfpz);
bx2= zeros(nfpx,nfpy,nfpz);
by2= zeros(nfpx,nfpy,nfpz);
bz2= zeros(nfpx,nfpy,nfpz);
bx3= zeros(nfpx,nfpy,nfpz);
by3= zeros(nfpx,nfpy,nfpz);
bz3= zeros(nfpx,nfpy,nfpz);
bx4= zeros(nfpx,nfpy,nfpz);
by4= zeros(nfpx,nfpy,nfpz);
bz4= 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;
ys1= -yw              : yw/(nspy-1)  : 0.0;  % layer of Cu just below vacuum gap (top is y=0).
zs= 0.0              : zw/(nspz-1)  : zw;
ys2= -2.0*yw              : yw/(nspy-1)  : -yw;  % layer of glass below Cu
ys3= 2.0*yw              : yw/(nspy-1)  : 3.0*yw;  % layer of glass above vacuum gap
ys4= 3.0*yw              : yw/(nspy-1)  : 4.0*yw;  % layer of cu above top glass


%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.98 ;
%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;
r1=sqrt((x(i)-xs(is))^2+(y(j)-ys1(js))^2+(z(k)-zs(ks))^2);

tmpx = alpha*gA^2*Ncu*(vy*(z(k)-zs(ks))-vz*(y(j)-ys1(js)))*(1/lamc+1/r1)*exp(-r1/lamc)/r1^2;
bx1(i,j,k) = bx1(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/r1)*exp(-r1/lamc)/r1^2;
by1(i,j,k) = by1(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)-ys1(js))-vy*(x(i)-xs(is)))*(1/lamc+1/r1)*exp(-r1/lamc)/r1^2;
bz1(i,j,k) = bz1(i,j,k) + tmpz*xw*yw*zw/nspx/nspy/nspz*1E12; 
% gives B in units of picoT


r2=sqrt((x(i)-xs(is))^2+(y(j)-ys2(js))^2+(z(k)-zs(ks))^2);

tmpx = alpha*gA^2*Ngl*(vy*(z(k)-zs(ks))-vz*(y(j)-ys2(js)))*(1/lamc+1/r2)*exp(-r2/lamc)/r2^2;
bx2(i,j,k) = bx2(i,j,k) + tmpx*xw*yw*zw/nspx/nspy/nspz*1E12;  
%  gives B in units of picoT

tmpy = alpha*gA^2*Ngl*(vz*(x(i)-xs(is))-vx*(z(k)-zs(ks)))*(1/lamc+1/r2)*exp(-r2/lamc)/r2^2;
by2(i,j,k) = by2(i,j,k) + tmpy*xw*yw*zw/nspx/nspy/nspz*1E12;  
% gives B in units of picoT

tmpz = alpha*gA^2*Ngl*(vx*(y(j)-ys2(js))-vy*(x(i)-xs(is)))*(1/lamc+1/r2)*exp(-r2/lamc)/r2^2;
bz2(i,j,k) = bz2(i,j,k) + tmpz*xw*yw*zw/nspx/nspy/nspz*1E12; 
% gives B in units of picoT


r3=sqrt((x(i)-xs(is))^2+(y(j)-ys3(js))^2+(z(k)-zs(ks))^2);

tmpx = alpha*gA^2*Ngl*(vy*(z(k)-zs(ks))-vz*(y(j)-ys3(js)))*(1/lamc+1/r3)*exp(-r3/lamc)/r3^2;
bx3(i,j,k) = bx3(i,j,k) + tmpx*xw*yw*zw/nspx/nspy/nspz*1E12;  
%  gives B in units of picoT

tmpy = alpha*gA^2*Ngl*(vz*(x(i)-xs(is))-vx*(z(k)-zs(ks)))*(1/lamc+1/r3)*exp(-r3/lamc)/r3^2;
by3(i,j,k) = by3(i,j,k) + tmpy*xw*yw*zw/nspx/nspy/nspz*1E12;  
% gives B in units of picoT

tmpz = alpha*gA^2*Ngl*(vx*(y(j)-ys3(js))-vy*(x(i)-xs(is)))*(1/lamc+1/r3)*exp(-r3/lamc)/r3^2;
bz3(i,j,k) = bz3(i,j,k) + tmpz*xw*yw*zw/nspx/nspy/nspz*1E12; 
% gives B in units of picoT


r4=sqrt((x(i)-xs(is))^2+(y(j)-ys4(js))^2+(z(k)-zs(ks))^2);

tmpx = alpha*gA^2*Ncu*(vy*(z(k)-zs(ks))-vz*(y(j)-ys4(js)))*(1/lamc+1/r4)*exp(-r4/lamc)/r4^2;
bx4(i,j,k) = bx4(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/r4)*exp(-r4/lamc)/r4^2;
by4(i,j,k) = by4(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)-ys4(js))-vy*(x(i)-xs(is)))*(1/lamc+1/r4)*exp(-r4/lamc)/r4^2;
bz4(i,j,k) = bz4(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

phix1=0.0;
phiy1=0.0;
phiz1=0.0;
phiInfx1=0.0;
phiInfz1=0.0;
for i=1:nfpx;
for j=1:nfpy;
for k=1:nfpz;
bx(i,j,k) = bx1(i,j,k)+bx2(i,j,k)+bx3(i,j,k)+bx4(i,j,k);
by(i,j,k) = by1(i,j,k)+by2(i,j,k)+by3(i,j,k)+by4(i,j,k);
bz(i,j,k) = bz1(i,j,k)+bz2(i,j,k)+bz3(i,j,k)+bz4(i,j,k);

Binfx(i,j,k) = -0.0287*gA^2*lamc*vz*(Ncu*exp(-y(j)/lamc) + Ngl*exp(-(y(j)+yw)/lamc) ...
    - Ngl*exp(-(2.0-y(j))/lamc) - Ncu*exp(-(3.0-y(j))/lamc) )*1E12;  % pT

%for infinite slab, By=0

Binfz(i,j,k) = 0.0287*gA^2*lamc*vx*(Ncu*exp(-y(j)/lamc) + Ngl*exp(-(y(j)+yw)/lamc) ...
    - Ngl*exp(-(2.0*yw-y(j))/lamc) - Ncu*exp(-(3.0*yw-y(j))/lamc) )*1E12;  %pT


if isnan(bx(i,j,k))
     fprintf(' bx Nan : %i %i %i %d \n', i,j,k)
end
if isnan(Binfx(i,j,k))
     fprintf(' Binfx Nan : %i %i %i %d \n', i,j,k)
end


phiInfx1 = phiInfx1 - 1.832E8*Binfx(i,j,k)/1E12/660.0/nspx/nspy/nspz;   % average rotation rad/m
phiInfz1 = phiInfz1 - 1.832E8*Binfz(i,j,k)/1E12/660.0/nspx/nspy/nspz;
phix1 = phix1 - 1.832E8*bx(i,j,k)/1E12/660.0/nspx/nspy/nspz;
phiy1 = phiy1 - 1.832E8*by(i,j,k)/1E12/660.0/nspx/nspy/nspz;
phiz1 = phiz1 - 1.832E8*bz(i,j,k)/1E12/660.0/nspx/nspy/nspz;

end
end
end


% put in the rest of the slab for the phi calculation
xslab = 50.0-2*xw;
zslab = 500.0-2*zw;
phix2=0.0;
phiy2=0.0;
phiz2=0.0;
phiInfx2=0.0;
phiInfz2=0.0;

for j=1:nfpy;

phiInfx2 = phiInfx2 - 1.832E8*Binfx(nfpx,j,nfpz)/1E12/660.0/nspy;   % average rotation rad/m
phiInfz2 = phiInfz2 - 1.832E8*Binfz(nfpx,j,nfpz)/1E12/660.0/nspy;
phix2 = phix2 - 1.832E8*bx(nfpx,j,nfpz)/1E12/660.0/nspy;
phiy2 = phiy2 - 1.832E8*by(nfpx,j,nfpz)/1E12/660.0/nspy;
phiz2 = phiz2 - 1.832E8*bz(nfpx,j,nfpz)/1E12/660.0/nspy;


end
% now add up the pieces

phiInfx = (4.0*phiInfx1*xw*zw + 2.0*phiInfx1*xw*zslab + phiInfx2*xslab*zslab)/5.0;   % rotation rad, wgt average by width
phiInfz = (4.0*phiz1*xw*zw + 2.0*phiInfz1*xw*zslab + phiInfz2*xslab*zslab)/5.0;
phix = (4.0*phix1*xw*zw + 2.0*phix1*xw*zslab + phix2*xslab*zslab)/5.0;
phiy = (4.0*phiy1*xw*zw + 2.0*phiy1*xw*zslab + phiy2*xslab*zslab)/5.0;
phiz = (4.0*phiz1*xw*zw + 2.0*phiz1*xw*zslab + phiz2*xslab*zslab)/5.0;

bxx22 = bx(:,2,2);
byx22 = by(:,2,2);
bxx42 = bx(:,4,2);
byx42 = by(:,4,2);
%bxx1 = squeeze(bx(:,1,:));
%bxx2 = squeeze(bx(:,2,:));
%bxx3 = squeeze(bx(:,3,:));
%bxx4 = squeeze(bx(:,4,:));
%bxx5 = 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)')



fprintf(' dPhi/dz (rad) bx: %d   by: %d   bz: %d\n', phix,phiy,phiz)
fprintf(' dPhi/dz (rad) binfx: %d   binfz: %d\n', phiInfx,phiInfz)

save fielddata1.mat
%save xfieldout bx -ASCII