% load sst anomalies made using read-ersst.m (detrended)
% etc.

load ssta_short.mat ssta3 Lat Lon N M;  % ssta3 is detrended sst anomalies 1900-1996
                                   % data is 30S to 60N; missing data is NaN 
                                   % Lat is the latitudes (19) and Lon is
                                   % the Longitudes (36), so data is 
                                   % 10 deg long by 5 deg latitude

time = 1950:1/12:1996.99;

% get only locations with data and put into an 2-d matrix z(x,t)

% do north pacific only
min_lat = +20; 
max_lat =+60;
min_lon = 120;
max_lon = 280;

% make mask
mask=repmat(0,[length(Lat) length(Lon)]);

for j=1:length(Lat)
      if Lat(j) >= min_lat & Lat(j) <= max_lat
              for k=1:length(Lon)
                        if Lon(k) >= min_lon & Lon(k) <= max_lon
                          mask(j,k)=1;
                        else
                        end
              end
      else
      end
end

ssta=ssta3;
for i=1:N
    ssta(i,:,:)=squeeze(ssta3(i,:,:)).*mask;
end

    
[data index] = three_two(ssta);  % 'index' is the vector of index locations that have data


% get eofs and pcs

c=data*data'/N;
[v lamda]=eig(c);
v = fliplr(v);
lamda = fliplr(flipud(lamda));
p = v' * data;
Lamda=diag(lamda);
Lamda=abs(Lamda);

% Scale the eofs and pcs so that the amplitude is in the eofs, and the pcs
% have unit variance.

for i = 1:M
    v(:,i)=v(:,i)*sqrt(Lamda(i));
    if(std(p(i,:)')~=0)
    p(i,:)=p(i,:)/std(p(i,:)');
    end
end


dummy = repmat(NaN,[19 36]);   % an array used when converting  data stored in a vector (by location index) 
                               % back to the Lat vs. Lon matrix format

figure(1)
for i=1:4
    subplot(2,2,i)
    dummy(index) = v(:,i);  % map the location data in column back to a matrix 
                               % that is dimension Lat by Lon
    pcolor(Lon,Lat,dummy)
    caxis([-1 1]); colorbar
    title(['EOF ' int2str(i)])
end

figure(2)

for i=1:2
    subplot(2,1,i)
    mstart = (i-1)*2+1; mend =mstart+1; 
    plot(time,p(mend,:)',time,p(mstart,:)')
    title(['PC ' int2str( (i-1)*2 + 1) ' and ' int2str( (i-1)*2 + 2)])
end

figure(3)
subplot(2,1,1)
plot(1:M, cumsum(Lamda)/trace(lamda));grid on;
axis([1 40 0 1])
xlabel('Eigenvalue Number'); ylabel('Cumulative Fraction of total variance')
subplot(2,1,2)
plot(1:M,Lamda/trace(lamda),'x');grid on;
axis([1 20 0 1])
xlabel('Eigenvalue Number'); ylabel('Fraction of total variance')

% regress (project) onto normalized pc1 (called pc_n for the north pacific pc)

pc_n=p(1,:);

[data2 index2]=three_two(ssta3);
sstr=data2*pc_n'/N;
dummy(index2)=sstr;
sst_regress_n=dummy;

% Now plot how the entire SST field looks as a function of the strength of
% the first EOF
figure(4);
pcolor(Lon,Lat,sst_regress_n);caxis([-1 1]); colorbar
title('SST projected onto PC1 from the north Pacific data')

