% -----------------------------------------------------------------
data = load('cool5adcp');

% 0.01 to convert to SI units
u = 0.01*data.U.avg; 
v = 0.01*data.V.avg;
t =      data.dnum;
z =     -data.Z;

if 1
  % run a 3-point boxcar filter over the 20-min ADCP averages so that
  % we will have roughly comparable 1-hour averaged data
  uraw = u;
  vraw = v;
  tmpu = u';
  tmpv = v';
  for k=1:size(tmpu,2)
    tmpu2 = conv(tmpu(:,k),[1 1 1]/3);
    tmpv2 = conv(tmpv(:,k),[1 1 1]/3);
    tmpu(:,k) = tmpu2(2:end-1);
    tmpv(:,k) = tmpv2(2:end-1);
  end
  u = tmpu';
  v = tmpv';
end

if 1
  % work with small time range to limit cov matrix size
  daylim = [datenum([2001 7 13 0 0 0]); datenum([2001 8 6 0 0 0])];
  clip = find(t<daylim(1)|t>daylim(2));
  u(:,clip) = [];
  v(:,clip) = [];
  t(clip) = [];
end

if 1
  % don't trust the near surface bins
  clip = find(z>-6);
  u(clip,:) = [];
  v(clip,:) = [];
  z(clip) = [];
end

% -----------------------------------------------------------------
% Compute the EOFs for the raw ADCP data %%%%%%%%%%%%%%%%%%%%%%%%%%  

% detrend (transpose twice because detrend works down columns
up = detrend(u')';
vp = detrend(v')';

D = up + sqrt(-1)*vp;

[U,S,V] = svd(D,0);

percentvar = 100*diag(S)/trace(S)

E = U;
A = S*V';

% plot the mode structures
M = 1:2;
figure(1)
subplot(121)
plot(abs(E(:,M)),z)
title('Magnitude')
legend('mode 1','mode 2')
subplot(122)
plot(180/pi*angle(E(:,M)),z)
title('Angle')
legend('mode 1','mode 2')

% plot single mode reconstruction
figure(2)
K = 1:1;
Dhat = U(:,K)*A(K,:);
subplot(411)
pcolorjw(t,z,abs(Dhat))
title('abs mode 1')
subplot(412)
pcolorjw(t,z,angle(Dhat))
title('angle mode 1')
K = 2:2;
Dhat = U(:,K)*A(K,:);
subplot(413)
pcolorjw(t,z,abs(Dhat))
title('abs mode 2')
subplot(414)
pcolorjw(t,z,angle(Dhat))
title('angle mode 2')

% component time series
figure(3)
K = 1:1;
Dhat = U(:,K)*A(K,:);
uhat = real(Dhat);
vhat = imag(Dhat);
subplot(311)
pcolorjw(t,z,up)
title('u data')
subplot(312)
pcolorjw(t,z,uhat)
title('u mode 1 only')
K = 2;
Dhat = U(:,K)*A(K,:);
uhat = real(Dhat);
vhat = imag(Dhat);
subplot(313)
pcolorjw(t,z,uhat)
title('u mode 2 only')