// ESERCITAZIONE BENCHMARK DI LAUGENBUCH

clc, clear, xdel(winsid())

disp('Opzioni di plotting:');
disp('1) Plot flusso veloce');
disp('2) Plot flusso veloce sul piano x');
disp('3) Plot flusso veloce sul piano y');
disp('4) Plot flusso termico');
disp('5) Plot flusso termico sul piano x');
disp('6) Plot flusso termico sul piano y');
disp('7) Plot potenza del core ');
disp('8) Termina il file');

I=fscanfMat('input.txt')

nx=I(1,1); // nodi lungo x
ny=I(1,2); // nodi lungo y
t=I(1,3); // durata del transitorio
p=I(1,4); // passo temporale

A=fscanfMat('flux.txt');
B=fscanfMat('powermap.txt');

x=(1:nx);
y=(1:ny);

for u=1:20
r=input('Cosa si vuole plottare? ');

if r==1
for i=1:t*2/p+2
if mod(i,2)==1;
i1=1+(i-1)*nx;
i2=i1+nx-1;    
flux1=A(i1:i2,1:ny);
plot3d(X,Y,flux1);
xlabel('X Coordinate [cm]')
ylabel('Y Coordinate [cm]')
zlabel('Mean Fast Flux [n/cm^2*s]')
//axis([0,nx,0,ny,0,1.0e+15])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(i-1)/2*p ) ); 
xpause(1000*0.5);
end
end
end

if r==2
disp('valore min = 1, valore max = nx');
nxv=input('Piano x da plottare = ');    
for i=1:t*2/p+2
if mod(i,2)==1;    
iv=nxv+(i-1)*nx;
flux2=A(iv,1:ny);
plot2d1:ny,flux2','b');
xlabel('Y Coordinate [cm]')
ylabel('Mean Fast Flux [n/cm^2*s]')
//axis([0,ny,0,2.0e+15])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(i-1)/2*p ) ); 
xpause(1000*0.5);
end 
end
end

if r==3
disp('valore min = 1, valore max = ny');
nyv=input('Piano y da plottare = ');
for i=1:t*2/p+2
if mod(i,2)==1;
i1=1+(i-1)*nx;
i2=i1+nx-1;
flux3=A(i1:i2,nyv);
plot2d(1:nx,flux3','b');
xlabel('X Coordinate [cm]')
ylabel('Mean Fast Flux [n/cm^2*s]')
//axis([0,nx,0,2.0e+15])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(i-1)/2*p ) ); 
xpause(1000*0.5);
end
end
end

if r==4
for j=1:t*2/p+2
if mod(j,2)==0;
j1=nx+1+(j-2)*nx;
j2=j1+nx-1;       
flux4=A(j1:j2,1:ny);
plot3d(X,Y,flux4);
xlabel('X Coordinate [cm]')
ylabel('Y Coordinate [cm]')
zlabel('Mean Thermal Flux [n/cm^2*s]')
//axis([0,nx,0,ny,0,1.8e+14])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(j-2)/2*p ) ); 
pause(1000*0.5);
end
end
end

if r==5
disp('valore min = 1, valore max = nx');
nxt=input('Piano x da plottare = ');
for j=1:t*2/p+2
if mod(j,2)==0;
jt=nxt+(j-2)*nx;
flux5=A(jt,1:ny);
plot2d(1:ny,flux5,'r');
xlabel('Y Coordinate [cm]')
ylabel('Mean Thermal Flux [n/cm^2*s]')
//axis([0,ny,0,1.8e+14])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(j-2)/2*p ) ); 
xpause(1000*0.5);
end
end
end

if r==6
disp('valore min = 1, valore max = ny');
nyt=input('Piano y da plottare = ');
for j=1:t*2/p+2
if mod(j,2)==0;
j1=nx+1+(j-2)*nx;
j2=j1+nx-1;  
flux6=A(j1:j2,nyt);
plot2d(1:nx,flux6,'r');
xlabel('X Coordinate [cm]')
ylabel('Mean Thermal Flux [n/cm^2*s]')
//axis([0,nx,0,1.8e+14])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(j-2)/2*p ) ); 
pause(1000*0.5);
end
end
end

if r==7
for z=1:t/p+1
z1=1+(z-1)*nx;
z2=z1+nx-1;    
flux7=B(z1:z2,1:ny);
plot3d(X,Y,flux7);
xlabel('X Coordinate [cm]')
ylabel('Y Coordinate [cm]')
zlabel('Core power [-]')
//axis([0,nx,0,ny,0,20])
xtitle(sprintf('                                                Transient Diffusion = %10.5e s ',(z-1)*p ) );
xpause(1000*0.5);
end
end

if r==8
break
end
end