clc;
clear all;
L = 25;
N = 51;
x = 0:L/N:L;
y=x;
z=x;

T = zeros(N+1,N+1, N+1);
 
water=[0.5985 4.18 1]; % thermal_conductivity , capacity , density
wall=[85 3.57 0.535]; % Lithium
hydrogen_plasma=[0.1805 14.3 1.67E-9];
thermal_conductivity=[water(1) wall(1) hydrogen_plasma(1)]; % W/(meters*K)
capacity = [water(2) wall(2) hydrogen_plasma(2)]; %specific J/g*Kelvin
density = [water(3) wall(3) hydrogen_plasma(3)]; % gramm/cm^3 
k = thermal_conductivity./(capacity.*density); %%cm^2/seconds

iter=0;
deltax = L/N;
deltay = L/N;
deltaz = L/N;
deltat = 1E-8; %time step = 10ns
n = (deltat) / (6*(deltax^2));


inner_radius=6.5/2;
outer_radius=19.4/2;
height = 11.3/2;
avg_radius=(outer_radius+inner_radius)/2;
fusion_threshold=20;
thickness=1;
yc=L/2;
xc=L/2;
zc=L/2;
T_c_o=30;
for i=2:N
    for j=2:N
        for l=2:N
            d=sqrt((x(i)-xc)^2+(y(j)-yc)^2);
            h=abs(z(l)-zc);
            inside_torus=(d<=outer_radius-thickness && d>=inner_radius+thickness && h<=height-thickness);
            if inside_torus
                T(i,j,l)=T_c_o-2*abs(d-avg_radius)-2*h; % tens of million Kelvin E7
            end
        end
    end
end

% [X, Z] = meshgrid(x,z);
% [Y, Z] = meshgrid(y,z);
[X2, Y2] = meshgrid(x,y);
[X,Y,Z] = meshgrid(x,y,z);
T_y_z=reshape(T(round(N/2),:,:),[N+1 N+1]);
T_x_z=reshape(T(:,round(N/2),:),[N+1 N+1]);
T_x_y=reshape(T(:,:,round(N/2)),[N+1 N+1]);

%surf(X, Z, T_x_z);
%isosurface(X,Y,Z,T); % Use Temperature to set mark color
caxis([0 T_c_o])
axis equal
o=colorbar;
colorTitleHandle = get(o,'Title');
titleString = 'Temperature (ten million K)';
set(colorTitleHandle ,'String',titleString);



p=0; % the place, we start nowhere
deltaV=deltax*deltay*deltaz; % m^3
while iter <= 100
    Tnew = T;
    fusion_rate=0;
    for i=2:N %sets fusion_rate based on existing T
        for j =2:N
            for l=2:N
                h=abs(z(l)-zc);
                d=sqrt((x(i)-xc)^2+(y(j)-yc)^2); 
                inside_torus=(d<=outer_radius-thickness && d>=inner_radius+thickness && h<=height-thickness);
                if (T(i,j,l) >= fusion_threshold) && inside_torus
                    fusion_rate=fusion_rate+fusion_of_T(T(i,j,l),N)*deltaV; %temperature gain per volume * Volume
                end
            end
        end
    end
    for i=2:N
      for j =2:N
         for l=2:N
            h=abs(z(l)-zc);
            d=sqrt((x(i)-xc)^2+(y(j)-yc)^2);     
            inside_torus=(d<=outer_radius-thickness && d>=inner_radius+thickness && h<=height-thickness);
            inside_wall=(d>outer_radius-thickness && d<=outer_radius+thickness)||(d>inner_radius-thickness && d<=inner_radius+thickness)||(h<=height+thickness && h>height-thickness);
            inside_water=(d<inner_radius || d>outer_radius || h>height);

            if inside_wall
                p=2; 
                T(i,j,l)=T(i,j,l)+fusion_rate/(N^3*(x(i)^2+y(j)^2+z(l))^2); % heat the walls based on fusion(neutron flux) rate
            elseif inside_torus
                p=3;
            elseif inside_water
                p=1;
            end
         
            Tnew(i,j,l)=n*k(p)*(T(i-1, j,l) + T(i+1,j,l) + T(i, j+1,l) + T(i, j-1,l) + T(i, j,l-1) + T(i, j,l+1) - 6*T(i,j,l)) + T(i,j,l);
        
            
         end
       end
        
    end
    T = Tnew;
    iter = iter+1;
    T_y_z=reshape(T(round(N/2),:,:),[N+1 N+1]);
    T_x_z=reshape(T(:,round(N/2),:),[N+1 N+1]);
    T_x_y=reshape(T(:,:,round(N/2)),[N+1 N+1]);
    %surf(X, Z, T_x_z);
    %surf(X2, Y2, T_x_y);
    isosurface(X,Y,Z,T)
    pause(0.01)
    
end

 %%
pcolor(X, Y, T_x_y)
shading interp
hold on
contour (X, Y, T_x_y, 30, 'k');




