%% 3d Heat Equation Trial with Polyval
p = [-0.000000219377116  0.000019037643653 -0.000695680017162  0.013659012592243 -0.150372230012195  0.815648271244546 -0.696731668656325  ];
x = linspace(1.033625312069473, 23.633886177386074,20);
y = polyval(p,x); % plot(x,y,'r')
[X,Y,Z] = cylinder(y);
%Tcylinder = cylinder(y);%.*100
% surf(Z,Y,X)
% axis square, set(gca,'color',[1 1 1]*0.8), grid off
% colormap([1 1 1])
Xh = X;
X = Z.*8; %All are 20x21 matrices
Z = Xh.*8;
Y = Y.*8;
 
%Shifting by 12
shifter = 12.5*(ones(20,21));
mildershifter = 11.5*(ones(20,21));
 
Xshifted = mildershifter+X;
Yshifted = shifter+Y;
Zshifted = shifter+Z;
 
surf(Zshifted,Yshifted,Xshifted);
%surf(Z,Y,X)
%axis square, set(gca,'color',[1 1 1]*0.8), grid off
%colormap([1 1 1])
xlabel('X-axis_Zmatrix')
ylabel('Y-axis_Ymatrix')
zlabel('Z-axis_Xmatrix')
 
%Resized capsule
 
%Ploting ZYX faces the proper direction
 
%To match y needs to go plus 12
%Z needs they all need to go plus 12
 
clc;
clear all;
L = 25;
N = 51;
x = 0:L/N:L;
y=x;
z=x;
 
T = zeros(N+1,N+1, N+1);

k =[0.143181818181818,44.5037828215398,7558309.95351953]; 
 
iteration=0;
f = 0;
deltax = L/N;
deltay = L/N;
deltaz = L/N;
deltat = 1E-8; 
n = (deltat) / (6*(deltax^2));
 
 
innerlimitradius=1;
outerlimitradius=9.9;
diskheight = 11.3/2;
meanlimit=(outerlimitradius+innerlimitradius)/2;

width=1;


ycenter=L/2;
xcenter=L/2;
zcenter=L/2;
Initialshieldtemp=25;


for i=2:N
    for j=2:N
        for l=2:N
            d=sqrt((x(i)-xcenter)^2+(y(j)-ycenter)^2);
            h=abs(z(l)-zcenter);
            shieldinterior= (d<=outerlimitradius-width && h<=diskheight-4.25);
            if shieldinterior
                T(i,j,l)=Initialshieldtemp+abs(d-meanlimit)-2*h; 
            end
        end
    end
end
 

[X2, Y2] = meshgrid(x,y);
[X,Y,Z] = meshgrid(x,y,z);
Tyz=reshape(T(round(N/2),:,:),[N+1 N+1]);
Tyz=reshape(T(:,round(N/2),:),[N+1 N+1]);
Tyy=reshape(T(:,:,round(N/2)),[N+1 N+1]);
 
p=0;
deltaV=deltax*deltay*deltaz;


shieldheatingrate=1.2; 


caxis([8 22])
axis equal
o=colorbar;
colorTitleHandle = get(o,'Title');
titleString = 'Temperature *10^2 Kelvin';
% R = 0;

while iteration <= 100
    Tnew = T;
        
    if iteration <= 25
     for i=2:N
      for j =2:N
         for l=2:N
            h=abs(z(l)-zcenter);
            d=sqrt((x(i)-xcenter)^2+(y(j)-ycenter)^2);     

            shieldedge=(d>outerlimitradius-width && d<=outerlimitradius+width)||(d>innerlimitradius-width && d<=innerlimitradius+width)||(h<=diskheight+width && h>diskheight-width);
            

            if shieldedge %True
                p=2; 
                T(i,j,l)=T(i,j,l)+shieldheatingrate/(N^3*(x(i)^2+y(j)^2+z(l))^2); % heat the walls based on fusion(neutron flux) rate
            elseif shieldinterior
                T(i,j,l)=T(i,j,l)+(1.2)/(N^3*(x(i)^2+y(j)^2+z(l))^2);
            end
            
         end
       end
     end
     end
    
    
    for i=2:N
      for j =2:N
         for l=2:N
            h=abs(z(l)-zcenter);
            d=sqrt((x(i)-xcenter)^2+(y(j)-ycenter)^2);     
            shieldinterior=(d<=outerlimitradius-width && d>=innerlimitradius+width && h<=diskheight-width);
            shieldedge=(d>outerlimitradius-width && d<=outerlimitradius+width)||(d>innerlimitradius-width && d<=innerlimitradius+width)||(h<=diskheight+width && h>diskheight-width);
            airsurrounding=(d<innerlimitradius || d>outerlimitradius || h>diskheight);
 
            if shieldedge %True false statement
                p=2; 
                T(i,j,l)=T(i,j,l)+shieldheatingrate/(N^3*(x(i)^2+y(j)^2+z(l))^2); % heat the walls based on fusion(neutron flux) rate
            elseif shieldinterior %True
                p=3;
            elseif airsurrounding %True 
                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
    
    
    if iteration >= 75
        
     for i=2:N
      for j =2:N
         for l=2:N
            h=abs(z(l)-zcenter);
            d=sqrt((x(i)-xcenter)^2+(y(j)-ycenter)^2);     
            shieldedge=(d>outerlimitradius-width && d<=outerlimitradius+width)||(d>innerlimitradius-width && d<=innerlimitradius+width)||(h<=diskheight+width && h>diskheight-width);          
 
            if shieldedge
                p=2; 
                T(i,j,l)=(9/10).*T(i,j,l); % Cooling heat shield
            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
    end
    
    T = Tnew;
    iteration = iteration+1;
    Tyz=reshape(T(round(N/2),:,:),[N+1 N+1]);
    Txz=reshape(T(:,round(N/2),:),[N+1 N+1]);
    Txy=reshape(T(:,:,round(N/2)),[N+1 N+1]);

    %view([-44 -4])
    isosurface(X,Y,Z,T) %-44 -6 azimuth for plotting
    title('3 Dimensional Heating of the Capsule on Reentry - Degrees Kelvin')
    f=f+1;
     saveas(gcf, ['3Dheatflow_',sprintf('%03d',f),'.png']) %Saving an image of each iteration in order to 
    pause(0.01)
    
    disp(iteration)
    
end