clear all;
close all;

load precip_sno.txt;

  % first three values are year,month, and day
  % datenum in this instance converts time into a single 
  % numerical value so that we know
  % that consecutive days are separated by a value of 1
  dates = datenum(precip_sno(:,1:3));
  data(:,1:7) = precip_sno(:,4:10);
  %data(:,1)=ben lomond peak
  %data(:,2)=ben lomond trail
  %data(:,3)=tony grove lake
  %data(:,4)=snowbird
  %data(:,5)=thaynes canyon
  %data(:,6)=timpanogos
  %data(:,7)=payson
  
  label = ['BLP';'BLT';'TGL';'SNW';'THC';'TIM';'PAY'];
  
  [nd,nn] = size(dates);

  %sum data over each season
  dateold = 0;
  total = zeros(18,7);
  iy=0;
  for i=1:nd
   %increment season by looking for break between may and oct
    if(dates(i)-dateold > 100)
              iy=iy+1;
              yr(iy,1)=iy+1989;
    end
     total(iy,:)=total(iy,:)+data(i,1:7);  
     dateold = dates(i);
  end
 
  ny = iy;
 
  %for convenience and to match notes, define total as x
  % columns are locations and rows are years

x = total;


figure(1)
plotmatrix(x);

xm = mean(x);
xs = std(x,1);

% anomalies
xp = x - ones(iy,1) * xm ;

% standardized anomalies

xstar = xp ./ (ones(iy,1) * xs);

% plot hovmuller diagram of standardized temporal anomalies
figure(2)
hold on
xa = 1:1:7;
[C,hc]=contour(xa,yr,xstar);
clabel(C,hc);
set(gca,'YDir','reverse')

% linear algebra way to compute cross correlations

rho_1 = transpose(xstar)*xstar / iy;

% matlab direct calculation
rho_2 = corrcoef(x);

figure(3)
hold on
[C,hc]=contour(rho_2,[0.99;0.95;0.90;0.85;0.80;0.75;0.70;0.65]);
clabel(C,hc);

%now have columns be years and rows be locations

xh = transpose(total);

xmh = mean(xh);
xsh = std(xh,1);

% anomalies
xph = xh - ones(7,1) * xmh ;

% standardized anomalies

xstarh = xph ./ (ones(7,1) * xsh);
% linear algebra way to compute cross correlations

s_1 = transpose(xstarh)*xstarh / 7;

% matlab direct calculation
s_2 = corrcoef(xh);

% plot hovmuller diagram of standardized spatial anomalies
figure(4)
hold on
[C,hc]=contour(yr,xa,xstarh);
set(gca,'YDir','reverse')
clabel(C,hc);

figure(5)
hold on
[C,hc]=contour(yr,yr,s_2);
clabel(C,hc);