clear all
clf

roomSize=49;

R=zeros(roomSize,roomSize,roomSize);
R2=zeros(roomSize,roomSize,roomSize);
tmax=300;

set(gca,'nextplot','replacechildren');
aviobj=avifile('PollDiffusion','quality',100,'compression','Cinepak')


figure('Position',[15 15 1000 1000 ])
hold on
 
for t=1:tmax
    if t<100;
    R(15,25,45:49)=40;
    R(35,25,46:49)=40;
    end
    
    for i=2:roomSize-1
        for j=2:roomSize-1
            for k=2:roomSize-1
                 R(i,j,k)=R(i-1,j,k)*.12+R(i+1,j,k)*.12+R(i,j+1,k)*.12+R(i,j-1,k)*.12+R(i,j,k+1)*.40+R(i,j,k-1)*.12;             
            end
        end
    end

    subplot(2,2,1)
    view(0,45)
    hold on
    surf(squeeze(R(:,15,:)))
    axis([0 50 0 50 0 50])    
    xlabel('altitude')
    ylabel('depth')
    zlabel('concentration')
    title('Concentration for cross section, y=25')
    pause(.01)	 

    subplot(2,2,2)
    hold on
    view(0,45)
    surf(squeeze(R(:,:,40)))
    axis([0 50 0 50 0 5])  
    xlabel('length')
    ylabel('depth')
    zlabel('concentration')
    title('constant elevation =40') 
    pause(.01)    

    subplot(2,2,3)
    hold on
    view(0,45)
    surf(squeeze(R(:,:,25)))   
    xlabel('length')
    ylabel('depth')
    zlabel('concentration')
    title('constant elevation= 25')
    axis([0 50 0 50 0 1])
    pause(.01)
    
    subplot(2,2,4)
    hold on
    view(0,45)
    surf(squeeze(R(:,:,10)))    
    xlabel('length')
    ylabel('depth')
    zlabel('concentration')
    title('constant elevation= 10')
    axis([0 50 0 50 0 1])
    pause(.1)

    F(j)=getframe;
    aviobj=addframe(aviobj,F(j));    
end
