%  eda08_10
% CAC Pacific Sea Surface Temperature (SST) dataset
% plot singular values

% load the data
D=load('cac_sst.txt');

% set up sizes of arrays
[I, J] = size(D);
NBLOCK=85;
NLON = NBLOCK-1;
NLAT = J-1;
IMAGES = floor(I/NBLOCK);
Dt = 1/12;

% arrays for SST, date, lats, lons
SST = zeros(NLON,NLAT,IMAGES);
thedate = zeros(IMAGES,1); % in month.year format
theyear = zeros(IMAGES,1); % year extracted from thedate
themonth = zeros(IMAGES,1); % month extracted from thedate
lats=zeros(NLAT,1);
lons=zeros(NLON,1);

% cut up data into SST arrays
lats = D(1,2:J)';
lons = D(2:NBLOCK,1);
for j=[1:IMAGES]
    k1 = 2+(j-1)*NBLOCK;
    k2 = k1+NLON-1;
    thedate(j)=D(k1-1,1);
    themonth(j)=floor(thedate(j));
    theyear(j)=floor(10000*(thedate(j)-themonth(j))+0.1);
    SST(:,:,j) = D(k1:k2,2:J);
end

% some definitions related to plotting
MONTHS=12;
YEARS=ceil(IMAGES/12);
YBLOCKSIZE = 6;
YBLOCKS = ceil(YEARS/YBLOCKSIZE);

% fold out into array
SSTV = zeros( IMAGES, NLAT*NLON );
for j=[1:IMAGES]
    for p = [1:NLAT]
    for q = [1:NLON]
        k = p + (q-1)*NLAT;
        SSTV(j,k) = SST( q, p, j );
    end
    end
end

% singular value decomposition
[U,S,V] = svd(SSTV,0);
s=diag(S);
Ns=length(s);
C=U*S;
F=V';

% plot singular values
figure(YBLOCKS+2);
clf;
plot( [1:Ns]', s, 'k-');
hold on;
xlabel('index, i');
ylabel('s(i)');
title('singular values, s(i)');

