Cij

function f=Cij(geom,tau,i,j,m,nGauss,track,sigma)
xi=i-m+1;
xj=j-m+1;
if strcmp(geom,'cyl')
if i<=j && sigma(i)~=0 && sigma(j)~=0
f=1/(sigma(i)*sigma(j))*(akin(3,tau)-akin(3,tau+sigma(i) ...
*Lx(xi,m,nGauss,track)) ...
-akin(3,tau+sigma(j)*Lx(xj,m,nGauss,track)) ...
+akin(3,tau+sigma(i)*Lx(xi,m,nGauss,track) ...
+sigma(j)*Lx(xj,m,nGauss,track)));
elseif i<=j && sigma(i)==0 && sigma(j)~=0
f=Lx(xi,m,nGauss,track)/sigma(j) ...
*(akin(2,tau)-akin(2,tau+sigma(j)*Lx(xj,m,nGauss,track)));
elseif i<=j && sigma(j)==0 && sigma(i)~=0
f=Lx(xj,m,nGauss,track)/sigma(i) ...
*(akin(2,tau)-akin(2,tau+sigma(i)*Lx(xi,m,nGauss,track)));
elseif i<=j && sigma(j)==0 && sigma(i)==0
f=Lx(xi,m,nGauss,track)*Lx(xj,m,nGauss,track)*akin(1,tau);
end
elseif strcmp(geom,'sph')
if i<=j && sigma(i)~=0 && sigma(j)~=0
f=1/(sigma(i)*sigma(j))*(exp(-tau)-exp(-(tau+sigma(i) ...
*Lx(xi,m,nGauss,track)))- ...
exp(-(tau+sigma(j)*Lx(xj,m,nGauss,track))) ...
+exp(-(tau+sigma(i)*Lx(xi,m,nGauss,track) ...
+sigma(j)*Lx(xj,m,nGauss,track))));
elseif i<=j && sigma(i)==0 && sigma(j)~=0
f=Lx(xi,m,nGauss,track)/sigma(j) ...
*(exp(-tau)-exp(-(tau+sigma(j)*Lx(xj,m,nGauss,track))));
elseif i<=j && sigma(j)==0 && sigma(i)~=0
f=Lx(xj,m,nGauss,track)/sigma(i) ...
*(exp(-tau)-exp(-(tau+sigma(i)*Lx(xi,m,nGauss,track))));
elseif i<=j && sigma(j)==0 && sigma(i)==0
f=Lx(xi,m,nGauss,track)*Lx(xj,m,nGauss,track)*exp(-tau);
end
else
error('Undefined geometry!');
end
Di

function f=Di(geom,i,m,nGauss,track,sigma)
x=i-m+1;
if strcmp(geom,'cyl')
if sigma(i)~=0
f=Lx(x,m,nGauss,track)/sigma(i)-1/sigma(i)^2 ...
*(akin(3,0)-akin(3,sigma(i)*Lx(x,m,nGauss,track)));
else
f=pi*Lx(x,m,nGauss,track)^2/4;
end
elseif strcmp(geom,'sph')
if sigma(i)~=0
f=Lx(x,m,nGauss,track)/sigma(i)-1/sigma(i)^2 ...
*(1-exp(-(sigma(i)*Lx(x,m,nGauss,track))));
else
f=Lx(x,m,nGauss,track)^2/2;
end
else
error('Undefined geometry!');
end
Keff

function f=keff(pij,vol,nusigmaf,sigmas,ref)
n=length(vol);
for i=1:n
for j=1:n
SQ(i,j)=sigmas(j)*vol(j)*pij(j,i)/vol(i);
FQ(i,j)=nusigmaf(j)*vol(j)*pij(j,i)/vol(i);
end
end
SQ=eye(n)-SQ;
tmp=aleig(SQ,FQ,eps);
f=1/tmp{2};
function L

function f=L(x,y,sigma,vol)
f=0;
for i=1:length(sigma)
if (i-1/2>=x & i+1/2<=y)
f=f+sigma(i)*vol(i);
end
end