matlabcode

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.

%Checing of input paramters

%Checking of input parameter L
if L<=0
    error ('The length(L) must be positive number!');
end

%Checking of input parameter N
if (isinteger(N)==0)&&(N<=0)
    error ('The number of cells along each axis(N) must be intiger number greater than 0!');
end

%Checking of input parameter tmax
if tmax<0
    error('The simulation period(tmax) must be greater than 0')
end

%Checking of flow parameter
if (q<=0)&&(q>(1/8))
    error('The flow parameter (q) cannot be greater than 1/8!');
end

%Checking validity of Data file- if data file exisits

if (exist('Datafile')==0)
    error ('Sorry, Data files does not exist!');
end

%Checking if drain exists
if nargin==6
    if drain<=0
        error ('The number of drain must be positive values!');
    end
end


% 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
end

% Three subfunction called by simulate
% Author: Aleksandra Stempien
% Date: 08.05.2012


function rain= getrain (tnext,Data)
% Getrain get the amount of rain water, given in mm during the latest
% minute
%
% Call: rain= getrain (tnext,Data)
% Input Parameters:
% tnext             - the index to search for
% Data              - array containing two columns(  first with the number of minutes
%                      passed since simulation was started in increasing order, second rain fallen within the latest minute)
% Ouput Parameter:
% rain              - amount of rain in time period
%
% Author: Aleksandra Stempien
% Date: 08.05.2012


i = find(Data(:,1)==tnext,1);

if (isempty(i)==1)
    rain = 0 ;
else
    rain = Data (i,2);
end
end

function V = update (N,q,H,V,rain,I,J,drain)
% update calculates the depths of the water (from topography up to the
% water surface , after the latest minute
%
% Call: V = update (N,q,H,V,rain,I,J,drain)
% Input Parameters:
% N             - The number of cells along each axis.
% 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.
% H             - topology
% V             - initial depths of the water for during the latest minute
% rain          - amount of rain in time period
% I             - first index in the array drain
% J             - second index in the array drain
% drain         -The rate of water flow through the drain cell (in
% m/minute).
%
% Ouput Parameter:
% V             - updated depths of the water for during the latest minute
%
% Author: Aleksandra Stempien
% Date: 08.05.2012

% creating extended arrays of H an V, and put the initial values to 0
updateH=zeros(N+2);
updateV=zeros(N+2);

%boarding the highest point
top= max(H(:));
updateH(:)=top;
 %extending H and V arrays in the was that all the orginal cells have now 8
 %neighboring cells and coping the data from H and V matrixes to the
 %updated one
for i = 2:N+1
    for j = 2:N+1
       updateV(i,j) = V(i-1,j-1);
       updateH(i,j) = H(i-1,j-1);
    end;
end;


%Calculateing the avarage amount of rain, that has been falling, and the avarage
%decrease of the water depth  during 20 periods; the results are given in meters
AverageRain=((rain/20)/1000);
if nargin==8
    AverageDrain=drain/20;
end

%looping in each of twenty 3-seconds periods
for tp=1:20

    dV = zeros(N+2); % initiaily putting elements of array dV to 0

%calculating the chages of V due to the differences in heights of water
%surface in neighboring cells

% Adding avarage rain to the amount of rain
    for j=2:N+1
        for k=2:N+1
            updateV(j,k)=updateV(j,k)+AverageRain;
        end
    end


    %Calculating the amount of water in each cell
    for i=2:N+1
        for j=2:N+1
            %Defining of the parameter S storing the sum of differences
            %Defining the parameter SurfaceDif sotring the differences
            %between the water surface of its neighborning cells; initial
            %puting paramters to 0
            S=0;
            SurfaceDif=zeros(3);
            %Calculating how much water is above the water surface of its
            %neigborning cells ( difference between the cell)
            for r=-1:1
                for c=-1:1

                   if((updateH(i,j)+updateV(i,j)-(updateH(i+r,j+c)+updateV(i+r,j+c))) < 0)
                        SurfaceDif(r+2,c+2) = 0; %if difference is negative it is put to 0
                    else
                        SurfaceDif(r+2,c+2) = (updateH(i,j)+updateV(i,j)-(updateH(i+r,j+c)+updateV(i+r,j+c)));
                    end;
                    S=S+SurfaceDif(r+2,c+2).*q; % Sum of differences
                end
            end

            % Calculating new dV value dependently on the depth of the water
            if (S>updateV(i,j))
            for r=-1:1
                for c=-1:1
                    if (r==0&&c==0)
                        continue;
                    else
                        %incresing array dV by adding surface difference
                        dV(i+r,j+c)=dV(i+r,j+c)+ (((q*SurfaceDif(r+2,c+2)).*updateV(i,j))./S);
                    end;
                end;
            end
            dV(i,j)=dV(i,j)-updateV(i,j);

            else
                for r=-1:1
                    for c=-1:1
                        if ((r==0)&&(c==0))
                            continue;
                        else
                            dV(i+r,j+c)=dV(i+r,j+c)+ (q*SurfaceDif(r+2,c+2));
                        end
                    end
                end
                dV(i,j)=dV(i,j)-S;
            end
        end
    end
            % Updating V value
            updateV=updateV+dV;

            %If the there is non blocked drain cell, the water depth is
            %decreased
            if nargin==8
               if ((updateV(I,J)-AverageDrain)>=0)
                   updateV(I,J)=updateV(I,J)-AverageDrain;
               else
                   updateV(I,J)=0;
               end
            end
end


% Saving calulated values into matrixes
for i=2:N+1
    for j = 2:N+1
        V(i-1,j-1) = updateV(i,j);
        H(i-1,j-1) = updateH(i,j);
    end;
end
end

function show(L,N,q,H,V,tnext,rain,I,J,drain)
% show the topography and the water surface
%
% Call: show(L,N,q,H,V,tnext,rain,I,J,drain)
% Input Parameters:
% L             - Length of domain along each axis (in meters).
% N             - The number of cells along each axis.
% 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.
% H             - topology
% V             - initial depths of the water for during the latest minute
% tnext         - current minute of simulation
% rain          - amount of rain in time period
% I             - first index in the array drain
% J             - second index in the array drain
% drain         -The rate of water flow through the drain cell (in  m/minute).
%
%
% Author: Aleksandra Stempien
% Date: 08.05.2012

%Converting the values H and V into decimeters
H = H * 10;
V = V * 10;

clf reset;
%Creating Z anc C indexes to the colors that are to be used for i andj
%indexes cell on the plotted surface
C=zeros(N);
Z=zeros(N);

for i=1:N
    for j=1:N
        C(i,j)=ceil((3*V(i,j))/(V(i,j)+0.25));
    end
end

%Creating colormaps according to the given arrays
if (max(C(:))<2)
    colormap ([0.5 0.25 0; 0 0.8 1.0]);
elseif (2<=max(C(:)) && max(C(:))<3)
    colormap ([0.5 0.25 0; 0 0.8 1.0; 0 0.5 0.65]);
elseif (3<=max(C(:)))
    colormap([0.5 0.25 0; 0 0.8 1.0; 0 0.5 0.65; 0 0.2 0.3]);
end

%Calculating the maximal topographical height
HV=H+V;
top=max(HV(:));

xvalues=1:N;
yvalues=1:N;

step=L/N;

xvalues=xvalues*step;
yvalues=yvalues*step;

% limiting the ranges of axis

axis([0 L 0 L 0 top+2]);
hold on

layer1=surf(xvalues,yvalues,H,Z);


layer2=surf(xvalues,yvalues,H+V,C);

%mark the drain
if nargin==10
    hold on
    plot3(I,J,0,'b*');
end

% Writing required text in given positions, labeling the plot, adding the
% title, creating the grid

if(nargin == 10)
    hold on
    plot3(I,J,1,'marker','*')
    legend('soil', 'water', 'drain');
    titleText = sprintf('Simulation with One Drain Cell (drain = %1.1f m/min). Flow parameter q = %1.2f', drain, q);
    else
    legend('soil', 'water');
    titleText = sprintf('Simulation with Blocked or No Drain Cells. Flow parameter q = %1.2f', q);
end;


title on;
title(titleText);

grid on

xlabel('m');
ylabel('m');
zlabel('dm');

timeText = sprintf('time (minutes) = %1.0f', tnext);
rainText = sprintf('rain (mm/min) = %f', rain);
text(0.1*L,L, top+1.8,timeText);
text(0.1*L,L, top+0.9,rainText);


end