% We don't want to print on each timestep -> would be too time consuming
if(mod(t,500)==0)
    
        
    % draw 2d temperature
    figure(1), clf, box
    imagesc([min(min(X)) max(max(X))],[min(min(Y)) max(max(Y))], T(iny,inx)), axis equal
    colormap jet    
    axis tight
    ylabel('y [m]')
    xlabel('x [m]'), colorbar,  drawnow
    axis([0,Lx,0,Ly])
    print -dpng Temperature2dfins
    pause(0.1)
    
    
    % draw 2d velocity
    figure(2), clf, box
    imagesc([min(min(X)) max(max(X))],[min(min(Y)) max(max(Y))], U), axis equal, hold on
    qiny = iny(1:5:length(iny)); qinx = inx(1:5:length(inx));
    quiver( X(qiny,qinx), Y(qiny,qinx), U(qiny,qinx), V(qiny,qinx),4, 'w')
    colormap jet    
    axis tight
    ylabel('y [m]')
    xlabel('x [m]'), colorbar,  drawnow
    axis tight
    axis([0,Lx,0,Ly])
    
    print -dpng Velocity2dfins
    pause(0.1)
    


figure(3), clf, box
Tm = sum(T(iny,inx).*U)./(sum(U)); Um = mean(U);

            
end