%
% Program to Fit Exponential Flux Model to
%                    Settling Curves
%
% Bob Newell & Lisa Hopkins, 1996
%
%-------------------------------------------------------
% Insert your settling test data in arrays:
%    conc = row vector of concentrations (mg/l) one
%       for each test done
%    timeht = matrix of settling test data where
%       column 1 = time in hours of height readings
%       columns 2... = heights in metres at corresponding
%           times in column 1 (enter NaN if no reading)
% then adjust k1 and k2 in second section appropriately
%-------------------------------------------------------
% Data from Loganholme for Book Figures 8.6 and 8.7
%
conc = [ 4500 6400 1400 900 700 ];
timeht = [ 0 1 1 1 1 1; ...
       0.05 0.981964 0.998 0.89 0.84 0.85; ...
       0.083333 0.971944 0.99 0.79 0.69 0.69; ...
       0.116667 0.956914 0.985 0.695 0.5 0.49; ...
       0.15 0.936874 0.975 0.605 0.35 0.29; ...
       0.183333 0.926854 0.97 0.53 0.24 0.15; ...
       0.216667 0.906814 0.965 0.475 0.2 0.12; ...
       0.25 0.891784 0.96 0.425 0.17 0.1; ...
       0.283333 0.871743 0.955 0.385 0.15 NaN; ...
       0.316667 0.856713 0.95 0.35 0.135 NaN; ...
       0.35 0.831663 0.94 0.325 0.115 NaN; ...
       0.383333 0.811623 0.935 0.3 0.11 0.075; ...
       0.416667 0.796593 0.93 0.28 0.105 NaN; ...
       0.45 0.776553 0.923 0.26 NaN NaN; ...
       0.483333 0.756513 0.913 NaN NaN NaN; ...
       0.516667 0.736473 0.905 NaN NaN NaN; ...
       0.55 0.718437 NaN NaN NaN NaN; ...
       0.583333 0.701403 NaN 0.215 NaN NaN; ...
       0.616667 0.681363 NaN NaN NaN NaN; ...
       0.65 0.666333 0.87 NaN NaN NaN; ...
       0.75 0.624248 NaN NaN NaN NaN; ...
       0.85 0.584168 0.815 NaN NaN NaN; ...
       0.95 NaN NaN NaN NaN NaN; ...
       1.05 0.521042 0.76 NaN NaN NaN; ...
       1.15 0.490982 NaN NaN NaN NaN ];
%--------------------------------------------------------
%  plot the settling curves
%
[ m, n ] = size( timeht );
mark = [ 'ko', 'kx', 'k+', 'k*' ];
plot( timeht(:,1), timeht(:,2), mark(7:8) );
set( gca, 'fontsize', 16 );
xlabel( 'Time in hours' );
ylabel( 'Height in metres' );
hold on;
for i = 3:n,
  k =  2 * ( i - 3 - 4 * floor( ( i - 3 ) / 4 ) ) + 1;
  plot( timeht(:,1), timeht(:,i), mark(k:k+1) );
end;
hold off;
%
% add concentrations
ax = axis;
hold on;
for i = 2:n,
  for j = m:-1:1,
    if ~isnan( timeht(j,i) ),
      k = j; break;
    end;
  end;
  x(i) = timeht(k,1) + 0.02;
  y(i) = timeht(k,i);
  if x(i) + 0.3 > ax(2),
    x(i) = ax(2) - 0.3;
    y(i) = y(i) - 0.08; 
  end;
  if i > 2 & y(i-1)-y(i) < 0.08,
    y(i-1) = y(i-1) + 0.05;
  end;
end;
for i = 2:n,
  string = [ num2str(conc(i-1)) ' mg/l' ]
  h = text( x(i), y(i), string );
  set( h, 'fontsize', 16 );
end;
hold off;
%
print -f1 -deps2 setdata.eps
%-----------------------------------------------------------
%  estimate initial velocity from settling curves
%  and plot them with the regression line
%
%  k1, k2 define range of points to use to calculate
%         the velocity and then average
%
k1 = 4;
k2 = 6;
vel = zeros( k2-k1+1, n-1 );
for i = 1:k2-k1+1,
  vel(i,:) = ( ones(1,n-1) - timeht(k1+i-1,2:n) ) ...
            / timeht(k1+i-1,1);
end;
vbar = mean( vel );
flux = conc .* vbar / 1000;
%
figure;
plot( conc, log(vbar), '*' );
set( gca, 'fontsize', 16 );
xlabel( 'Concentration in mg/l' );
ylabel( 'Natural log of velocity' );
ax = axis
hold on;
c = [ ones(5,1) conc' ] \ log(vbar')
[ ax(1); ax(2) ]
[ c(1)+c(2)*ax(1); c(1)+c(2)*ax(2) ]
plot( [ ax(1); ax(2) ], [ c(1)+c(2)*ax(1); c(1)+c(2)*ax(2) ] );
hold off;
%
print -f2 -deps2 setfit.eps
%------------------------------------------------------------
% plot fitted flux curve with experimental points
%
figure;
plot( conc, flux, '*' );
set( gca, 'fontsize', 16 );
sylabel( '\16Flux in kg/m^2/h' );
xlabel( 'Concentration in mg/l' );
ax = axis;
%
hold on
c(1) = exp( c(1) );
c(2) = - c(2);
co = ax(1):100:ax(2);
fl = c(1) * co .* exp( - c(2) * co ) / 1000 ;
plot( co, fl, '-' );
hold off
%
printsto -f3 -deps2 setflux.eps
%------------------------------------------------------------
%  print out the flux model coefficients
%
disp( 'flux = conc * velocity' );
disp( '     = conc * C1 * exp( - C2 * conc )' );
disp( [ 'where C1 = ' num2str( c(1) ) ' m/h' ] );
disp( [ '      C2 = ' num2str( c(2) ) ' l/mg' ] );
%
%-------------------------------------------------------------
% That's All
