%
% Figures 15.5 to 15.8
%
load syvab4;
[ n, m ] = size( dat );
tim = ( 0:n-1 )' * 0.1;
%
% plot time series
%
figure;
plot( tim, dat, '-' );
set( gca, 'fontsize', 16 );
xlabel( 'Time in hours' );
ylabel( 'NH4 and COD in mg/l, Flow in Ml/d' );
h = text( 20, 600, 'Flow' );
set( h, 'fontsize', 16 );
h = text( 20, 260, 'COD' );
set( h, 'fontsize', 16 );
h = text( 20, 80, 'NH4' );
set( h, 'fontsize', 16 );
%
print -f1 -deps2 fig15-5.eps
%
% fill in missing data
%
for i = 1:m,
  for j = 1:n,
    if isnan( dat(j,i) ),
      dat(j,i) = dat(j-1,i);
	end;
  end;
end;
%
mean1 = mean( dat(:,1) );
mean2 = mean( dat(:,2) );
mean3 = mean( dat(1:1586,3) );
disp( [ num2str( mean1 ) '  +-  ' ...
            num2str( std( dat(:,1) ) ) ] );
disp( [ num2str( mean2 ) '  +-  ' ...
            num2str( std( dat(:,2) ) ) ] );
disp( [ num2str( mean3 ) '  +-  ' ...
            num2str( std( dat(1:1586,3) ) ) ] );
%
figure;
subplot( 1, 2, 1 );
resdist( dat(:,1)-mean1 );
set( gca, 'fontsize', 16 );
xlabel( 'Amplitude' );
ylabel( 'Flow in Ml/d' );
subplot( 1, 2, 2 );
[ h, f, sp ] = frspec( dat(1:5:2000,1)-mean1, 0.5, 0, 0 );
size(f)
plot( f(1:40), sp(1:40), '-' );
set( gca, 'fontsize', 16 );
axis( [ 0 2 0 400000 ] );
xlabel( 'Radians/hour' );
ylabel( 'Magnitude' );
%
print -f2 -deps2 fig15-6.eps
%
figure;
subplot( 1, 2, 1 );
resdist( dat(:,2)-mean2 );
set( gca, 'fontsize', 16 );
xlabel( 'Amplitude' );
ylabel( 'COD in mg/l' );
subplot( 1, 2, 2 );
[ h, f, sp ] = frspec( dat(1:5:2000,2)-mean2, 0.5, 0, 0 );
plot( f(1:40), sp(1:40), '-' );
set( gca, 'fontsize', 16 );
axis( [ 0 2 0 20000 ] );
xlabel( 'Radians/hour' );
ylabel( 'Magnitude' );
%
print -f3 -deps2 fig15-7.eps
%
figure;
subplot( 1, 2, 1 );
resdist( dat(1:1586,3)-mean3 );
set( gca, 'fontsize', 16 );
xlabel( 'Amplitude' );
ylabel( 'NH4 in mg/l' );
subplot( 1, 2, 2 );
[ h, f, sp ] = frspec( dat(1:5:1586,3)-mean3, 0.5, 0, 0 );
plot( f(1:40), sp(1:40), '-' );
set( gca, 'fontsize', 16 );
axis( [ 0 2 0 4000 ] );
xlabel( 'Radians/hour' );
ylabel( 'Magnitude' );
%
print -f4 -deps2 fig15-8.eps
%
