clf
clear

N=50;
L=20;
k=1;
T=zeros(N+1,N+1);
T(:,1)=0;
T(:,N+1)=0;
T(1,:)=0;
T(N+1,:)=0;
f=0;
solidsugar=50;
% Parameters for the glass
for i=1:N+1
    for j=1:N+1
        if j==1
            T(i,j)=-1;
        end
        if j==N+1
            T(i,j)=-1;
        end
        if i==N+1
            T(i,j)=-1;
        end
        if i==1
            T(i,j)=-1;
        end
        if i<N-4*j
            T(i,j)=-1;
        end
        if i<-3*N +4*j
            T(i,j)=-1;
        end
        if i<((j-N/2)/5)^2
            T(i,j)=-1;
        end
    end
end

% Building the sugar cube
Tcube=zeros(N+1,N+1);
for i=1:N
    for j=1:N
        if (N/2)-3<=j && (N/2)+3>=j && N-6<=i && N>=i
            Tcube(i,j)=1000+solidsugar;
        end
    end
end

x=zeros(1,N+1);
dx=L/N;
dy=L/N;
dt=1e-3;

fpd=2; %frames per one-pixel drop

% Dissolving the cube
v=dt/(dx^2)
Numit=20000;
Cubemin=N;
for k=1:Numit
    Tnew=T;
    if Tcube(N-6,N/2) >= 100
        fprintf('%d,%d:%d ',N-6,N/2,Tcube(N-6,N/2));
        fprintf('%d,%d:%d ',N-8,N/2,Tcube(N-8,N/2));
    end
    for i=2:N
        for j=2:N
            Ttemp=0;
            neighbors=0;
            if Tcube(i,j)<1000 && T(i,j)>=0 %Don't do anything if it's in the sugar cube or outside the glass.
                if T(i+1,j)>=0 && Tcube(i+1,j)<1000
                    Ttemp=Ttemp+T(i+1,j);
                    neighbors=neighbors+1;
                end
                if Tcube(i+1,j)>=1000
                    Ttemp=Ttemp+solidsugar;
                    Tcube(i+1,j)=Tcube(i+1,j)-v*solidsugar; %this introduces more error as v gets larger
                    neighbors=neighbors+1;
                end
                if T(i,j+1) >=0 && Tcube(i,j+1)<1000
                    Ttemp=Ttemp+T(i,j+1);
                    neighbors=neighbors+1;
                end
                if Tcube(i,j+1)>=1000
                    Ttemp=Ttemp+solidsugar;
                    Tcube(i,j+1)=Tcube(i,j+1)-v*solidsugar; %this introduces more error as v gets larger
                    neighbors=neighbors+1;
                end
                if T(i,j-1)>=0 && Tcube(i,j-1)<1000
                    Ttemp=Ttemp+T(i,j-1);
                    neighbors=neighbors+1;
                end
                if Tcube(i,j-1)>=1000
                    Ttemp=Ttemp+solidsugar;
                    Tcube(i,j-1)=Tcube(i,j-1)-v*solidsugar; %this introduces more error as v gets larger
                    neighbors=neighbors+1;
                end
                if  T(i-1,j) >=0 && Tcube(i-1,j)<1000
                    Ttemp=Ttemp+T(i-1,j);
                    neighbors=neighbors+1;
                end
                if Tcube(i-1,j)>=1000
                    Ttemp=Ttemp+solidsugar;
                    Tcube(i-1,j)=Tcube(i-1,j)-v*solidsugar; %this introduces more error as v gets larger
                    neighbors=neighbors+1;
                end
                
                Tnew(i,j)=Tnew(i,j)+v*(Ttemp-neighbors*T(i,j));
                
            end
        end
    end
    
    % Dropping the cube
    if mod(k,fpd)==0
        if Cubemin>2
            for i2=2:N
                for j2=2:N
                    if  Tcube(i2,j2)>=1000
                        Tcube(i2-1,j2)=Tcube(i2,j2);
                        Tcube(i2,j2)=0;
                        if i2-1<Cubemin
                            Cubemin=i2-1
                        end
                    end
                end
            end
        end
    end
    
    T=Tnew;
    
    % Plotting the mess
    if mod(k+9,10)==0
        f=f+1;
        xx=0:L/N:L;
        yy=0:L/N:L;
        [X,Y]=meshgrid(yy,xx);
        pcolor(X,Y,T) %2d
        shading interp %2d
        %surf(X,Y,T); %3d
        %axis([0 L 0 L 0 N/2]) %3d
        saveas(gcf, ['louchesugar_',sprintf('%03d',f),'.png'])
        pause(0.1)
    end
end

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