function V = simulate(L,N,tmax,q,Datafile,drain)
%SIMULATE performs a hydrological simulation
%
% It is called this way:
% V = simulate(L,N,tmax,q,Datafile), if there is no drain,
% V = simulate(L,N,tmax,q,Datafile,drain), otherwise.
%
% A call of this function will produce a hydrological simulation, and
% return the final water distribution in V.
%
% Input parameters:
% L = Length of domain along each axis (in meters).
% N = The number of cells along each axis.
% tmax = The simulation period is [0, tmax] (in minutes).
% q = Flow parameter (must be <= 1/8).
% When the heights of water surface in two neighbouring cells
% differ, the heights will change q times this difference
% during the next 3 seconds, unless this more than empties a cell.
% Cells further away will not be affected during these 3 seconds.
% Datafile = A MAT-file with rain data in a kx2 array called Data.
% drain = The rate of water flow through the drain cell (in m/minute).
% Author: Jørgen B. Sand, March, 2012.
% Check of the input parameters starts here
% nargin checks how many inputarguments there are:
switch nargin
case 5
numberOfArgumentsInput=5;
case 6
numberOfArgumentsInput=6;
otherwise
end
% Checks if any of the input are negative or zero and prints
% an error message in the case
if (L <= 0 || N <= 0 || tmax <= 0)
error('N,L or tmax is negative')
end
% Checks if tmax is an integer by checking the modulus
% to 1, which on integers always are 0.
if ~(mod(tmax,1) == 0)
error('tmax is not an integer ')
end
if (0> q || q >= (1/8))
error('q is negative or more than 1/8')
end
%The test of the input parameters starts here.
% We obtain the topology (in meters) as the NxN array H
[H,I,J] = topography(N);
% (I,J) are the array indexes of the drain cell, if there is an unblocked one.
% The initial water level (water surface minus topology)
% is stored (in meters) in the NxN array V
V = zeros(N);
% All the rain data is loaded into the kx2 array Data
load(Datafile);
if Data(end,1) < tmax
error('Data file does not contain data for the whole simulation period')
end
% The simulation is performed and visualized
figure()
if nargin==6 % There is an unblocked drain cell.
for tnext=1:tmax
% Get the amount of rain in the time period (tnext-1, tnext]
rain = getrain(tnext,Data);
% Compute the water distribution at time tnext
V = update(N,q,H,V,rain,I,J,drain);
% Show the topology H and the water surface H+V
show(L,N,q,H,V,tnext,rain,I,J,drain); drawnow;
end
else % There is no unblocked drain cell.
for tnext=1:tmax
% Get the amount of rain in the time period (tnext-1, tnext]
rain = getrain(tnext,Data);
% Compute the water distribution at time tnext
V = update(N,q,H,V,rain);
% Show the topology H and the water surface H+V
show(L,N,q,H,V,tnext,rain); drawnow;
end
end
%This is where the code I've written starts.
function r = getrain(tnext,Data)
%GETRAIN fetchecs the rain data for the next minute
%
% It is called this way:
% r = getrain(tnext,data)
%
% A call of this function will return a floating value
% with some amount of rain that is to fall in the next minute.
%
% Input parameters:
% tnext = The next minute of the simulation
% Data = The array containing the data with the rain.
%
% Author: Anders Nannerup, may 2012.
for i=1:length(Data(:,1))
index = find(Data(i,1)==tnext);
if isempty(index)
r = 0;
else
r = Data(i,2);
return
end
end
end
function V = update(N,q,H,V,rain,I,J,drain)
%UPDATE Returns the current level of water in the simulation
%
% It is called this way:
% V = update(N,q,H,V,rain,I,J,drain)
%
% A call of this function will return a NxN array containing
% the levels of water in every indice (in metres)
%
% Input parameters:
% N = The size of the array
% q = The flow level of the water.
% H = The array containing the heighs of the earth in metres.
% V = The current water level in the simulation.
% rain = a floating number indicating the rain level
% I = The location on the x-axis of the drain if a such exist
% J = the location on the y-axis of the drain if a such exist
% drain = a floating number indicating how much water the drain can drain.
% (in metres / minute)
% Author: Anders Nannerup, may 2012.
%two new arrays are initialized
HTemp = zeros(N+2,N+2);
VTemp = zeros(N+2,N+2);
column = 1;
%Finds the max of H+V in the current column.
for row=1:N+2
if (row>=1 && row <=N) && (column>=1 && column <=N)
h2=max(H(row,column)+V(row,column));
end
%sets the value of V into VTemp.
for column=1:N+2
if ~((row>=N+2 || row<=1) || (column >= N+2 || column<=1))
HTemp(row,column)=H(row-1,column-1);
VTemp(row,column)=V(row-1,column-1);
end
%if a larger maximum of heights is found the value
% becomes this value.
if (row>=1 && row <=N) && (column>=1 && column <=N)
h1=max(H(row,column)+V(row,column));
end
if h1>h2
h2=h1;
end
h=h2;
%h is set along the edges.
HTemp(1,column)=h;
HTemp(row,1)=h;
HTemp(row,N+2)=h;
HTemp(N+2,column)=h;
VTemp(row,1)=0;
VTemp(row,N+2)=0;
VTemp(1,column)=0;
VTemp(N+2,column)=0;
end
end
%V and H are set to the temporary matrixes with edges.
H=HTemp;
V=VTemp;
%Calculates how much rain falls in every 3-sec period
calcRain = rain/(20*1000); %rain is meters/minute, but we need mm/3-sec
%If the is a drain also calulates how much
%water is drain every 3-sec period.
if numberOfArgumentsInput ==6
drainCell = drain/(20*1000);% Same as above.
end
%For every 3-sec period.
for period=1:20
%Increases the level of water with calcRain every iteraion.
V = V + calcRain;
% an array dV is initialized with the changes in V due to do it's neightbour cells
dV = zeros(27);
delta = zeros(3,3);
%For all cells
for i=2:N+1
for j=2:N+1
%for each cell, delta is the neighbour cells.
for r=-1:1
for c=-1:1
delta(r+2,c+2) = (H(i,j)+V(i,j)-(H(i+r,j+c)+V(i+r,j+c)));
if delta(r+2,c+2) < 0
delta(r+2,c+2)=0;
end
end
end
S = q*sum(sum(delta));
if S>V(i,j)
for r=-1:1
for c=-1:1
dV(i+r,j+c) = dV(i+r,j+c) + (q*delta(r+2,c+2)*(V(i,j))/S);
end
end
dV(i,j) = dV(i,j) - V(i,j);
else
dV(i+r,j+c) = dV(i+r,j+c) + q*delta(r+2,c+2);
dV(i,j) = dV(i,j) - S;
end
end
end
% V is increased/decreased with dV (dV might be negative)
V = V + dV;
% If there is one or more non-blocked drain, the water level is
% reduced in these cells, by the amount the drain can drain in a
% 3-sec period.
if numberOfArgumentsInput ==6
if ~(H(I,J) < 0)
if (V(I,J)-drainCell) > 0
V(I,J)= V(I,J)-drainCell;
else
V(I,J)=0;
end
end
end
end
% The edges is removed from V. Back to NxN
V = V(2:26,2:26);
return
end
function [] = show(L,N,q,H,V,tnext,rain,I,J,drain)
%SHOW shows the simulation in a window.
%
% It is called this way:
% show(L,N,q,H,V,tnext,rain,I,J,drain
%
% A call of this function will show and update
% the simulation.
%
% Input parameters:
% L = Length of domain along each axis (in meters).
% N = The number of cells along each axis.
% q = The flow level of the water.
% H = The array containing the heighs of the earth in metres.
% V = The current water level in the simulation.
% tnext = The next minute there is to be simulated
% rain = a floating number indicating the rain level
% I = The location on the x-axis of the drain if a such exist
% J = the location on the y-axis of the drain if a such exist
% drain = a floating number indicating how much water the drain can drain.
% (in metres / minute)
% Author: Anders Nannerup, may 2012.
H = H*10; %From meter to decimeter.
V = V*10; %From mm til decimeter.
%Two arrays c and Z are initialized.
Z = zeros(N,N);
C = zeros(N,N);
%The values of C is set.
for i=1:N
for j=1:N
C(i,j)=ceil(3*V(i,j)/(V(i,j)+0.25));
end
end
xvalues = linspace(0,L,N);
yvalues = linspace(0,L,N);
%The soil is plotted.
soil = surf(xvalues,yvalues,H,Z);
hold on
%the soil and water is plotted.
water = surf(xvalues,yvalues,H+V,C);
hold off
%The colormaps is defined.
brown = [ 0.5 0.25 0; 0 0.8 1.0];
blue = [ 0.5 0.25 0; 0 0.8 1.0; 0 0.5 0.65];
darkblue = [ 0.5 0.25 0; 0 0.8 1.0; 0 0.5 0.65; 0 0.2 0.3];
%Colormap is set
if max(C(:))<2
colormap(brown);
elseif 2<=max(C(:)) && max(C(:))<3
colormap(blue);
elseif 3<=max(C(:))
colormap(darkblue);
end
%The maximum topographical height is calculated and the axes are set.
maxH = max(H(:));
axis([0 L 0 L 0 (maxH+2)]);
%Text along the axes, rain, current minute etc. is set.
text(0.1*L,L,maxH +1.8,['time (minutes)= ',num2str(tnext), ' / ', num2str(tmax)]);
text(0.1*L,L,maxH +0.9,['rain (mm/min)= ', num2str(rain)]);
xlabel('m')
ylabel('m')
zlabel('dm')
%The title depends on whether or not there is a drain.
if numberOfArgumentsInput ==5
title(['Simulation with Blocked or No Drain Cells with flow parameter q = ', num2str(q)]);
else
title(['Simulation with One Drain Cellwith flow parameter q = ', num2str(q)]);
drain = text(I,J,0,'*','Color','b');
end
%The legend depends on wheter or not there is a drain.
if numberOfArgumentsInput ==5
legend1 = legend([soil,water],{'soil','water'});
else
legend1 = legend([soil,water,drain],{'soil','water','drain'});
end
set(legend1,'Location','NorthEastOutside');
end
end