Code covered by the BSD License

# Environmental Modeling - using MATLAB

### Ekkehard Holzbecher (view profile)

Springer Book, 2nd Edition

GaussianPuff
```function GaussianPuff
% Gaussian puff model
%    using MATLAB analytical solutions
%
%   \$Ekkehard Holzbecher  \$Date: 2006/08/21 \$
%-----------------------------------------------------------------------
d = 3;                      % dimensionality
Dx = 0.004; Dy = 0.001; Dz = 0.001;  % diffusivities
v = 1;                      % velocity
M = 1;                      % mass
xmin = 0.85; xmax = 1.15;   % x-axis interval
ymin = -0.1; ymax = 0.1;    % y-axis interval (used only for d>1)
zmin = -0.1; zmax = 0.1;    % z-axis interval (used only for d>2)
t = [1:4:20];               % time
gplot = 1;                  % plot option (=1 yes; =0 no)
gcont = 2;                  % contour plot option (=2 filled; =1 yes; =0 none)
gsurf = 0;                  % surface plot option (=1 yes; =0 no)
%       (only for d<3)
%-------------------------------execution------------------------------
switch d
case 1
x = linspace (xmin,xmax,100);
c = [];
for i = 1:size(t,2)
xx = x - v*t(i);
c = [c; (M/sqrt(4*pi*Dx*t(i)))*ones(1,size(x)).*...
exp(-(xx.*xx)/(4*Dx*t(i)))];
end
case 2
if size(t,2) > 1 t = t(1); display('first time level is considered only');
end
[x,y] = meshgrid (linspace(xmin,xmax,100),linspace(ymin,ymax,100));
xx = x - v*t;
c = (M/4/pi/t/sqrt(Dx*Dy))*ones(100).*exp(-0.25*(xx.*xx/Dx + y.*y/Dy)/t);
case 3
if size(t,2) > 1 t = t(1); display('first time level is considered only');
end
[x,y,z] = meshgrid (linspace(xmin,xmax,45),...
linspace(ymin,ymax,45),linspace(zmin,zmax,45));
xx = x -v*t;
c = (M/(sqrt(4*pi*t)^3)/sqrt(Dx*Dy*Dz))*ones(45,45,45)...
.*exp(-0.25*(xx.*xx/Dx + y.*y/Dy + z.*z/Dz)/t);
end

%---------------------------------output-------------------------------
if gplot
switch d
case 1
plot (c'); hold on; xlabel ('space');
ylabel ('concentration'); title ('1D Gaussian puff');
case 2
for i = 10:10:100
plot (c(:,i)); hold on;
end
xlabel ('x'); ylabel ('concentration');
title ('2D Gaussian puff along slices y=const)');
case 3
for i = 5:5:45
plot (c(:,i,23)); hold on;
end
xlabel ('x'); ylabel ('concentration');
title ('3D Gaussian puff along line x=0');
end
colorbar
end

switch gcont
case 1               % contours
figure;
switch d
case 1
contour (x,t,c,30);
grid on; xlabel ('space'); ylabel ('time');
case 2
contour (x,y,c,30);
xlabel ('x'); ylabel ('y');
case 3
end
title ('concentration');
case 2               % filled contours
figure;
switch d
case 1
contourf(x,t,c,30);
colorbar; xlabel ('space'); ylabel ('time');
case 2
contourf (x,y,c,30);
colorbar; xlabel ('x'); ylabel ('y');
case 3
xslice = [xmin 0.5*(xmin+xmax) xmax];
yslice = [ymax];
zslice = [zmin 0.5*(zmin+zmax) zmax];
slice (x,y,z,c,xslice,yslice,zslice);
xlabel ('x'); ylabel ('y'); zlabel ('z');
end
title ('concentration');
end
if gsurf
figure;
switch d
case 1
surf(x,t,c);
colorbar; xlabel ('space'); ylabel ('time');
case 2
surf (x,y,c);
colorbar; xlabel ('x'); ylabel ('y');
case 3
end
title ('concentration');
end

```

Contact us