% WATER_BALANCE - script to predict volume change in  a "lake"
% using predicted flow out and prescribed flow in -  
% Greg Hancock - College of William and Mary

clear           % clear the existing workspace variables

% SCRIPT PARAMETERS

V0 = 20;            % initial carboy volume (L)
time_max = 500;     % total time of simulation, in seconds
dt = 1;             % time step for simulation, in seconds
imax = time_max/dt; % total number of iterations

% CONSTANTS

v_to_h = 16.68;     % slope of relationship between carboy volume and height (L/cm)

% INFLOW HISTORY

% read prescribed inflow history (inflows.xls)

Q_in = xlsread('inflows.xls');  

% INITIAL VALUES

% volume of water in carboy (L)
volume = V0; 

% prescribed inflow from inflows.xls (L/s)
inflow = Q_in(1,2);  

% water height in carboy (mm)
height = volume * 16.68;

% outflow calculated from calibration equation (L/s)
  if volume <= 0
        outflow = 0;
    elseif volume < 5.5
        outflow = 0.007 * volume;
    else
        outflow = 0.0278 + 0.00193 * volume;
  end
  
  % initial time (secs)
time = 1;

% CREATE RECORDING MATRIX

record = zeros(time_max,5);

% record initial values
record(1,:) = [time volume inflow outflow height];

% TIME LOOP

for i = 2:imax
   
    % time
    time = i*dt;
     
    % update outflow (L/s)
    inflow = Q_in(i,2);
    
    % determine volume change (L)
    del_volume = (inflow - outflow)*dt;    
    
    % update volume (L)
    volume = volume + del_volume;
    
    % overflow
    if(volume > 20)
        volume = 20;
        'Carboy is overflowing'
    end
    
    % update height (mm)
    height = volume * 16.68;
  
    % update discharge (L/s)
    if volume <= 0
        outflow = 0;
    elseif volume < 5.5
        outflow = 0.007 * volume;
    else
        outflow = 0.0278 + 0.00193 * volume;
    end
    
    % update records
    record(i,:) = [time volume inflow outflow height];
    
end
 

% FIGURES

clf % clear figures

    figure(1)
    
    % plot inflow and outflow through time
    figure(1)  
    subplot(3,1,1)
    plot(record(:,1),record(:,3),'-b','LineWidth',3) % inflow
    hold on             
    plot(record(:,1),record(:,4),'-r','LineWidth',3) % outflow 
    
    title('Inflow and outflow')
    ylabel('Discharge (L/s)')
    legend('Inflow','Outflow')

    % plot water balance through time
    subplot(3,1,2) 
    plot(record(:,1),(record(:,3)-record(:,4)),'-g','LineWidth',3)
    hold on  
    
    hline = refline([0 0]);
    set(hline,'Color','k','LineStyle',':','LineWidth',2)
    
    ylabel('Water balance (L/s)')
    title('Water Balance')
    
    % plot volume through time
    subplot(3,1,3) 
    plot(record(:,1),record(:,2),'-k','LineWidth',3)
    hold on  
    
    title('Water Volume in Carboy')
    ylabel('Volume (L)')
    xlabel('Time (secs)')
  
    
% save volume records as comma-delimited file
xlswrite('record',record)



