L=100;
N=250;
T=zeros(N+1,N+1);
xx=0:L/N:L;yy=0:L/N:L;
high=1;
int=10;
h=0;

%T(:,1)=high;
%T(1,:)=high;
T(:,end)=high;
%T(end,:)=high;
dx=1;dy=1;
vxraw=30; 
vyraw=0;
k=1;

factor=0.5;
locations=35;
dimfact=0.1;

rate=1*rand(N+1,N+1);
rate=rate>=max(max(rate))*factor;
ind=find(rate<max(max(rate)));
sel=randperm(length(ind));
sel=sel(1:locations);
ind=ind(sel);
rate=max(max(rate))*ones(size(rate,1),size(rate,2));
rate(ind)=0;

for j=1:length(ind)
    [row,col]=ind2sub(size(rate),ind(j));
    rad=rand()*dimfact*mean([size(rate,1),size(rate,2)]);
    for a=1:size(rate,1)
        for b=1:size(rate,2)
            if (a-row)^2+(b-col)^2<rad^2
                rate(a,b)=0;
            end
        end
    end
end
zer=find(rate==0);

rows=zeros(1,length(zer));
cols=zeros(1,length(zer));
for i=1:length(zer)
    [rows(i),cols(i)]=ind2sub(size(rate),zer(i));
end

etamin=0.2;
dt=etamin*mean([dx,dy])^2/(40*k);

for iter=1:1000 
 Tup=[T(2:end,:);zeros(1,N+1)];
 Tdown=[zeros(1,N+1);T(1:end-1,:)];
 Tleft=[T(:,2:end),zeros(N+1,1)];
 Tright=[zeros(N+1,1),T(:,1:end-1)];
 if sign(vxraw)==1 && sign(vyraw)==1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vxraw.*dt/dx*(T-Tleft)-vyraw.*dt/dy*(T-Tdown);
 elseif sign(vxraw)==1 && sign(vyraw)==-1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vxraw.*dt/dx*(T-Tleft)-vyraw.*dt/dy*(Tup-T);
 elseif sign(vxraw)==-1 && sign(vyraw)==1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vxraw.*dt/dx*(Tright-T)-vyraw.*dt/dy*(T-Tdown);
 elseif sign(vxraw)==-1 && sign(vyraw)==-1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vxraw.*dt/dx*(Tright-T)-vyraw.*dt/dy*(Tup-T);
 elseif sign(vxraw)==1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vxraw.*dt/dx*(T-Tleft);
 elseif sign(vxraw)==-1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vxraw.*dt/dx*(Tright-T);
 elseif sign(vyraw)==1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vyraw.*dt/dy*(T-Tdown);
 elseif sign(vyraw)==-1
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T)-vyraw.*dt/dy*(Tup-T);
 else
 T=T+dt*k/dx^2.*(Tright+Tleft-2*T)+dt*k/dy^2.*(Tup+Tdown-2*T);
 end
 
 %T(:,1)=high;
 %T(1,:)=high;
 T(:,end)=high;
 %T(end,:)=high;
 T(zer)=0;
 
 %[X,Y]=meshgrid(yy,xx);
 %contourf(T);
 %hold all
 %plot(rows,cols,'ok','MarkerSize',3);
 
 if mod(iter,int)==0
 imagesc(T,[0,high]);
 set(gca,'YDir','normal');
 axis equal;
 axis tight;
 hold all
 
 
  plot(cols,rows,'ok','MarkerSize',5);
  saveas(gcf,['projectim',sprintf('%04d',h),'.png']);
  h=h+1;
  pause(0.0001) 
  disp(iter);
 end
 end