D=D39+D40+DO+DNa+DFe; %coefficienti di diffusione
F=F39+F40; %prodotti di fissione
CAPTURE=[Pu39capture+Pu40capture+Nacapture+Fecapture+Ocapture]; %sigma di cattura
FISSION=[Pu39fission+Pu40fission]; %sigma di fissione
ELREM=[Pu39elrem+Pu40elrem+Naelrem+Feelrem+Oelrem]; %sigma di scattering elastico(direclty coupled)
%si passa alla creazione delle matrici che serviranno per i sucessivi
%calcoli dei flussi normalizzati e non, di B, di K. inoltre si assegna un valore di partenza a
%B, e si calcola un flusso iniziale
Buckling=zeros(1,nmax);
K=zeros(1,nmax);
FLUSSIN=zeros(26,nmax)
FLUSSI=zeros(26,nmax)
FLUX0=zeros(26,1);
B1=0.00006;
nmax=1000;
for ii=1:26
FLUX0(ii)=1/F(ii);
end
%definisco la matrice dei coefficienti andando prima a riempire la diagonale di una matrice vuota
A=zeros(26,26);
for ii=1:26
for jj=1:26
for ii=jj
A(ii,jj)=D(ii).*B1+REMOVAL(ii);
end
end
end
% e successivamente aggiungendo i termini di scattering
% ottenendo una matrice triangolare inferiore
for jj=1:26
if ii~=jj
M(ii,jj)=A(ii,jj)-ELREM(ii)-INELASTIC(ii,jj);
end
end
for ii=1:26
for jj=1:26
COEFFICIENTI(ii,jj)=M(ii,jj)+A(ii,jj);
end
end
%risolvo il sistema:
FLUSSI=COEFFICIENTI^-1.*MATFIS; %flusso non normalizzato
K1=F.*FLUSSI; %calcolo del K
FLUSSIN=FLUSSI/K1; %calcolo del flusso normalizzato
%devo cambiare B in base a K1:
if K1>1
B2=B1+0.5*B1;
elseif K1<1
B2=B1-0.5*B1;
end
%definisco nuovamente la matrice dei coefficienti
A=zeros(26,26);
for ii=1:26
for jj=1:26
for ii=jj
A(ii,jj)=D(ii).*B2+REMOVAL(ii);
end
end
end
for jj=1:26
if ii~=jj
M(ii,jj)=A(ii,jj)-ELREM(ii)-INELASTIC(ii,jj);
end
end
for ii=1:26
for jj=1:26
COEFFICIENTI(ii,jj)=M(ii,jj)+A(ii,jj);
end
end
%risolvo nuovamente il sistema:
FLUSSI(:,2)=COEFFICIENTI\X; %flusso non normalizzato
K2=F'*FLUSSI(:,1); %calcolo del K
FLUSSIN(:,2)=FLUSSI(:,2)/K2; %calcolo del flusso normalizzato
epsilon=10^-6; %definisco la tolleranza
er=abs((K2-1)); %definisco l'errore
while er>epsilon & n<nmax
n=3;
B(n)=B(n-1)+((B(n-1)-B(n-2))*(1-K(n-1)))/(K(n-1)-K(n-2));
%ridefinisco la matrice dei coefficienti
A=zeros(26,26);
for ii=1:26
for jj=1:26
for ii=jj
A(ii,jj)=D(ii).*B(n)+REMOVAL(ii);
end
end
end
for jj=1:26
if ii~=jj
M(ii,jj)=A(ii,jj)-ELREM(ii)-INELASTIC(ii,jj);
end
end
for ii=1:26
for jj=1:26
COEFFICIENTI(ii,jj)=M(ii,jj)+A(ii,jj);
end
end
%risolvo nuovamente il sistema:
FLUSSI(:,n)=COEFFICIENTI\X; %flusso non normalizzato
Kn=F*FLUSSI(:,n); %calcolo del K
FLUSSIN(:,n)=FLUSSI(:,n)/Kn; %calcolo del flusso normalizzato
%ridefinisco l'errore
err=abs(Kn-1);
iterazioni=n;
n=iterazioni+1;
end
BFin=B(n); %buckling finale
R=pi/sqrt(BFin)
