%Hunt 1978 3-D dispersion solution Eqn.11.

clear all;
close('all')
[ X Y Z] = meshgrid( -1:0.05:3, -1:0.05:1, -1:0.05:1);
M3=100;
Dxx=.001;
Dyy=.0005;
Dzz=.0005;
sigma=.5;
V=0.04;

data = zeros( size(X,1), size(X,2), size(X,3));

disp(sprintf('size(X)    = [ %d %d %d]',size(X,1),    size(X,2),    size(X,3)));
disp(sprintf('size(Y)    = [ %d %d %d]',size(Y,1),    size(Y,2),    size(Y,3)));
disp(sprintf('size(Z)    = [ %d %d %d]',size(Z,1),    size(Z,2),    size(Z,3)));
disp(sprintf('size(data) = [ %d %d %d]',size(data,1), size(data,2), size(data,3)));

for t=1:10:501
disp(sprintf('t = %d',t));

yindex = 0;
for y=-1:0.05:1

    yindex = 1+yindex;

    xindex = 0;
    for x=-1:0.05:3

        xindex = xindex+1;

        zindex = 0;
        for z=-1:0.05:1

            zindex = zindex+1;

            R = sqrt(x^2+y^2*Dxx/Dyy+z^2*Dxx/Dzz); 

            if R==0
                data(yindex,xindex,zindex)=0;
            else
                data(yindex,xindex,zindex) = (M3*exp((x*V)/(2*Dxx))/(8*pi*sigma*R*sqrt(Dyy*Dzz)))*( exp(-R*V/(2*Dxx))*erfc((R-V*t)/(2*sqrt(Dxx*t)))+exp(R*V/(2*Dxx))*erfc((R+V*t)/(2*sqrt(Dxx*t))) ); 
            end

        end
    end
end
    
data(1,1,1)=1;%color scaling trick?

figure;
colormap(colorcube)
SLICE(X,Y,Z,data,[0],[0],[0])
title(sprintf('t = %d', t));

box on

clear data
view(150,30); daspect([1 1 1]);axis([-1,3,-1,1,-1,1])
a=gcf;%get current figure
saveas(a,sprintf('%04d%s',t,'.bmp'),'bmp');
%close

end
