% MatLab script for implementing data analysis described in 
% "Teaching Students to Look at Their Data - a Tutorial" by
% William Menke, June 2012.

clear all;

% load merged discharge and precipitation data
% column 1: time in inter count of days starting 1/1/2002
% column 2: discharge in cfs
% column 3: percipitation in inches
D = load('days_disc_precip.txt');

% break out into separate variables, metric conversion
t = D(:,1);
d = D(:,2)/35.3146;
p = D(:,3)*25.4;
N = length(t);
tmin = min(t); tmax=max(t);

% plot 100 day portion of discharge data
% the portion was selected to have several well-separated storms
tmin2 = 650;
tmax2 = 750;
figure(1);
clf;
hold on;
set(gca,'LineWidth',3);
axis( [tmin2, tmax2, 0, max(d) ] );
plot( t, d, 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('discharge (cfs)');

% plot portion of discharge and precipitation data
figure(2);
clf;
hold on;
subplot(2,1,1);
set(gca,'LineWidth',3);
hold on;
axis( [tmin2, tmax2, 0, max(d) ] );
plot( t, d, 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('discharge (cfs)');
subplot(2,1,2);
hold on;
set(gca,'LineWidth',3);
axis( [tmin2, tmax2, 0, max(p) ] );
plot( t, p, 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('precipitation (in)');

% output some statistics, for the year 2002
fprintf('2002 annual rainfall %f\n', sum(p(1:365))/1000 );
Sd = sum(d(1:365))*86400;
Sp = sum(p(1:365))*1.8e10/1000;
fprintf('2002 annual volume in m3 discharge=%.4e precip=%.4e\n', Sd, Sp );

% plot discharge data
figure(3);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin2, tmax2, 0, max(d) ] );
plot( t, d, 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('discharge (m^3/s)');

% create weighted running average of precipitation data
% using a weighted running average with weights that
% monotonically decline with index.  I use
% here a simple formula based on exponentials that has both
% a short and a long time scale.  I tuned the coefficients
% by trial-and-error to give a reasonable match to the discharge data.
Ns=3; Ns2=30; a=0.1;
f = exp( -[1:N]'/Ns ) + a*exp( -[1:N]'/Ns2 );
fn = sum(f);
f = f/fn;
ps = conv( p, f ); % implement via a convolution
ps = ps(1:N);

% plot weighted running average of precipitation data
subplot(2,1,2);
hold on;
set(gca,'LineWidth',3);
axis( [tmin2, tmax2, 0, max(ps) ] );
plot( t, ps, 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('smoothed precipitation (mm)');

% plot weights of weighted running average
figure(4);
clf;
subplot(2,1,1)
hold on;
set(gca,'LineWidth',3);
axis( [0, Ns2, 0, 1.5 ] );
z = zeros(1,Ns2+1); z(2)=1;
plot( [0:Ns2]', z, 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('precipitation');
subplot(2,1,2)
hold on;
set(gca,'LineWidth',3);
axis( [0, Ns2, 0, max(f) ] );
c=200000;
plot( [0:Ns2]', [0, f(1:Ns2)']', 'k-', 'LineWidth', 3 );
xlabel('time (days)');
ylabel('precipitation');


