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=21;   % 5 per lambda givse accruate value for lamc=1
nspy=21;
nspz=21;
nfpx=20;
nfpy=6;
nfpz=4;
xw=10.0*lamc ;
yw=1.0 ;
zw=10.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(ys1)
%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

avgbx1=0.0;
avgby1=0.0;
avgbz1=0.0;
avgbInfx1=0.0;
avgbInfz1=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*yw-y(j))/lamc) - Ncu*exp(-(3.0*yw-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


avgbInfx1 = avgbInfx1 + Binfx(i,j,k)/nfpx/nfpy/nfpz;   % average B-field (pT)
avgbInfz1 = avgbInfz1 + Binfz(i,j,k)/nfpx/nfpy/nfpz;
avgbx1 = avgbx1 + bx(i,j,k)/nfpx/nfpy/nfpz;
avgby1 = avgby1 + by(i,j,k)/nfpx/nfpy/nfpz;
avgbz1 = avgbz1 + bz(i,j,k)/nfpx/nfpy/nfpz;

end
end
end


% put in the rest of the slab for the phi calculation
xslab = 50.0-2*xw;
zslab = 500.0-2*zw;
avgbx2=0.0;
avgby2=0.0;
avgbz2=0.0;
avgbInfx2=0.0;
avgbInfz2=0.0;

for j=1:nfpy;

avgbInfx2 = avgbInfx2 + Binfx(nfpx,j,nfpz)/nfpy;   % average rotation rad/m
avgbInfz2 = avgbInfz2 + Binfz(nfpx,j,nfpz)/nfpy;   % here we're finding avgb in central region
avgbx2 = avgbx2 + bx(nfpx,j,nfpz)/nfpy;
avgby2 = avgby2 + by(nfpx,j,nfpz)/nfpy;
avgbz2 = avgbz2 + bz(nfpx,j,nfpz)/nfpy;

%fprintf('     %d %d %d ',nfpy, Binfx(nfpx,j,nfpz), avgbInfx2)

%fprintf(' y: %d   Bx(pT): %d   Bx_inf(pT): %d\n', y(j), bx(nfpx,j,nfpz),Binfx(nfpx,j,nfpz))


%fprintf('BxInf part 1 %d   ', -0.0287*gA^2*lamc*vz*Ncu*exp(-y(j)/lamc)*1E12)
%fprintf('BxInf part 2 %d   ', -0.0287*gA^2*lamc*vz*Ngl*exp(-(y(j)+yw)/lamc)*1E12)
%fprintf('BxInf part 3 %d   ', 0.0287*gA^2*lamc*vz*Ngl*exp(-(2.0*yw-y(j))/lamc)*1E12)
%fprintf('BxInf part 4 %d\n', 0.0287*gA^2*lamc*vz*Ncu*exp(-(3.0*yw-y(j))/lamc)*1E12)

end
% now add up the pieces

% three regions...corners, sides between corners, central slab, but we'll
% assume not much z-dependence on sides, so two regions, sides and center
phiInfx = -1.832E8*(2.0*avgbInfx1*xw*500.0 + avgbInfx2*xslab*500.0)/1E12/660.0/50.0;   % rotation (rad), wgt average by width
phiInfz = -1.832E8*(2.0*avgbInfz1*xw*500.0 + avgbInfz2*xslab*500.0)/1E12/660.0/50.0;
phix = -1.832E8*(2.0*avgbx1*xw*500.0 + avgbx2*xslab*500.0)/1E12/660.0/50.0;
phiy = -1.832E8*(2.0*avgby1*xw*500.0 + avgby2*xslab*500.0)/1E12/660.0/50.0;
phiz = -1.832E8*(2.0*avgbz1*xw*500.0 + avgbz2*xslab*500.0)/1E12/660.0/50.0;

phixcorner = -1.832E8*avgbx1*zw/1E12/660.0;
phiycorner = -1.832E8*avgby1*zw/1E12/660.0;
phizcorner = -1.832E8*avgbz1*zw/1E12/660.0;
phixside   = -1.832E8*avgbx1*zslab/1E12/660.0;
phixcenter = -1.832E8*avgbx2*500.0/1E12/660.0;

fprintf('  \n')
fprintf('  Average dPhi/dz (rad/m) \n')
fprintf('    corner                      central                  Semi-inf slab \n')
fprintf(' avgbx1 %d    avgbx2 %d    avgbInfx1 %d  avgbInfx2 %d\n ',avgbx1,avgbx2, avgbInfx1,avgbInfx2)
fprintf(' avgby1 %d    avgby2 %d    \n ',avgby1,avgby2)
fprintf(' avgbz1 %d    avgbz2 %d    avgbInfz1 %d  avgbInfz2 %d\n ',avgbz1,avgbz2, avgbInfz1,avgbInfz2)

fprintf('   \n')
fprintf(' Phi (rad) in corner phix: %d   phiy: %d   phiz: %d \n',phixcorner, phiycorner, phizcorner)
fprintf(' Phi (rad) in side   phix: %d \n',phixside)
fprintf(' Phi (rad) in center phix: %d \n',phixcenter)
fprintf(' Phi (rad) full calc phix: %d   phiy: %d   phiz: %d\n', phix,phiy,phiz)
fprintf(' Phi inf slab (rad) phiInfx: %d   phiInfz: %d\n', phiInfx,phiInfz)

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
subplot(2,2,1)
plot(x,bx(:,nfpy,nfpz),x,Binfx(:,nfpy,nfpz))
title(['Bx(x,',num2str(y(nfpy)),num2str(z(nfpz)),')(pT) lamc=',num2str(lamc),'mm'])
xlabel('y (mm)')
ylabel('Bx (pT)')

subplot(2,2,2)
plot(y,bx(nfpx,:,nfpz),y,Binfx(nfpx,:,nfpz))
title(['Bx(',num2str(x(nfpx)),',y,',num2str(z(nfpz)),')(pT) lamc=',num2str(lamc),'mm'])
xlabel('y (mm)')
ylabel('Bx (pT)')

dummy1=squeeze(bx(nfpx,nfpy,:))
dummy2=squeeze(Binfx(nfpx,nfpy,:))
subplot(2,2,3)
plot(z,dummy1,z,dummy2)
title(['Bx(',num2str(x(nfpx)),num2str(y(nfpy)),'z)(pT) lamc=',num2str(lamc),'mm'])
xlabel('y (mm)')
ylabel('Bx (pT)')


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 fielddata1.mat
%save xfieldout bx -ASCII