function  [I,J] = diffusion(radiance,diffcoeff,Nit,noise)
% Simulate the diffusion equation for Nit steps
%
% Copyright 2006 Paolo Favaro (p.favaro@hw.ac.uk)
% 
% School of Engineering and Physical Sciences
% Heriot-Watt University, Edinburgh, UK
% 
% Last revision: August 2006
%
% This program can be used only for research purposes.
% This program is distributed WITHOUT ANY WARRANTY; 
% without even the implied warranty of MERCHANTABILITY 
% or FITNESS FOR A PARTICULAR PURPOSE.

J = radiance;
I = radiance;
for i=1:Nit
    % diffuse (one step)
    if mod(i,2)==0
        div = divergence(J(:,:,i),diffcoeff/Nit);
    else 
        div = divergenceR(J(:,:,i),diffcoeff/Nit);
    end
    J(:,:,i+1) = J(:,:,i)+div;
end
I = J(:,:,Nit+1)+noise;
return