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

clear
close('all')
[x y z] = meshgrid(-1:0.05:3,-1:0.05:1,-1:0.05:1);

M3=1
Dxx=.001
Dyy=.0001
Dzz=.0001
sigma=.5
V=0.04


tindex=0
for t=1:1:51
    tindex=tindex+1

    data = M3*exp(-(x-V*t).^2/(4*Dxx*t)-y.^2/(4*Dyy*t)-z.^2/(4*Dzz*t))/(8*sigma*sqrt(pi^3*t^3*Dxx*Dyy*Dzz));

    well(tindex)=mean(data(20,40,:));%Vertically averaged concentrations at well (location given in matrix index values)  

    Ciso=10/t^(3/2) %Concentration of isosurface
    p = patch(isosurface(x,y,z,data,Ciso));
    isonormals(x,y,z,data,p);
    box on

    clear data
    set(p,'FaceColor','red','EdgeColor','none');
    alpha(1)%Transparency
    view(150,30); daspect([1 1 1]);axis([-1,3,-1,1,-1,1])
    camlight; lighting phong;
    
    a=gcf;%get current figure
    saveas(a,sprintf('%04d%s',t,'.bmp'),'bmp');
    close
    hold off

end
close