clear all
close all
clc

lamda = 0.6; ro = 1600; c = 0.92;
dx=1; dy=1; dz=1;
L1 = 10; L2 = 20; H = 20; K=3;
Th = 1200; Tc = 300; Te = 0; T0 = 0; t_end = 2500;
tau = t_end / 10;
x = 0:dx:(L2+5); y=0:dy:(L2+5); z=0:dz:21;
lx = length(x); ly=length(y); lz=length(z);

u = zeros(lx,ly,lz);
for k = 2:H+2
    for i = 1:lx
        for j = 1:ly
            u(i,j,k) = 300;
        end;
    end;
end;
for k = 2:H
    for i = 1:lx
        for j = 1:ly
            u(i,j,k) = 300;
            if (((i - lx / 2) * (i - lx / 2) + (j - ly / 2) * (j - ly / 2)) >= 100)
                u(i,j,k) = 1200;
            elseif (i>(lx/2-K) && i<=(lx/2+K) && j>(ly/2-K) && j<=(ly/2+K) && k >(lz/2-K) && k<=(lz/2+K))
                u(i,j,k) = 300;
            else u(i,j,k) = 0;
            end;
        end;
    end;
end;

for i = 1:lx
    for j = 1:ly
        u(i,j,1) = 0;
    end;
end;

un = u;

for t = 1:50:t_end
   for kk=2:H;
    for ii=2:(lx-1);
      for jj=2:(ly-1);
        un(ii,jj,kk) = lamda/ro/c*tau*(u(ii+1,jj,kk) -  2*u(ii,jj,kk)+u(ii-1,jj,kk)+u(ii,jj+1,kk)-2*u(ii,jj,kk)+u(ii,jj-1,kk)+u(ii,jj,kk+1)-2*u(ii,jj,kk)+u(ii,jj,kk-1))+u(ii,jj,kk);
           for i = 1:(lx-1);
             for j = 1:(ly-1);
               for k = 2:H;
                if (((i - lx / 2) * (i - lx / 2) + (j - ly / 2) * (j - ly / 2)) >= 100)
                    un(i,j,k) = 1200;
                elseif (i>(lx/2-K) && i<=(lx/2+K) && j>(ly/2-K) && j<=(ly/2+K))
                    un(i,j,k) = 300;
                end;
               end;
             end;
           end;
         end;
       end;
     end;
     for i = 1:(lx-1);
        for j = 1:(ly-1);
            for k = 2:H;
                if (((i - lx / 2) * (i - lx / 2) + (j - ly / 2) * (j - ly / 2)) >= 100)
                    un(i,j,k) = 1200;
                elseif (i>(lx/2-K) && i<=(lx/2+K) && j>(ly/2-K) && j<=(ly/2+K))
                    un(i,j,k) = 300;
             end;
          end;
        end;
    end;
   u=un;
    for i = 1:lx
        for j = 1:ly
            u(i,j,1) = u(i,j,2);
            u(i,j,21)=1200;
        end;
    end;
    colormap jet
    clf();
    xlabel('Координаты оси Ox');
    ylabel('Координаты оси Oy');
    zlabel('Координаты оси Oz');
    title([ ' t = ' num2str(t/100) ]);
    xslice = lx/2;
    yslice = ly/2;
    zslice = lz/2;
    slice(x,y,z,u,xslice,yslice,zslice);
    view(50,5);
    colorbar
    axis equal
    drawnow
end;