lon=[-350 -10]; lat=[-65 65];%lon=[-180 -160]; lat=[20 40];  % small sample of open ocean
time=[datenum(1950,1,15) datenum(2008,12, 15)];
forcd=rnc_Extract_NOAA_SST(lon, lat, time, 0);
forcd=rnc_forcdStats(forcd, 'SST');
[npgo_time, npgo] = rout_npgo;
npgo=npgo(1:708);
npgo_time=npgo_time(1:708);
% 31 is the timestep of the data in days
o=rnc_corr2(forcd.datenum, forcd.ano, npgo_time, npgo,31);
mfig; rnc_map(forcd.mean, forcd);
mfig; clf
subplot(2,1,1); rnc_map(o.corr, forcd);
% r1 is the regression, that gives you the actual size of the anomalies
% associated with variations in the index (e.g. NPGO)
subplot(2,1,2); rnc_map(o.r1, forcd);

% to make an eps file
lpr npgo_maps.eps
% print a bitmap PNG versino of the figure
print -dpng  npgo_maps.png




%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

lon=[-245 -110]; lat=[25 62];  % North Pacific
time=[datenum(1950,1,15) datenum(2008,12, 15)];
forcd=rnc_Extract_NOAA_SST(lon, lat, time, 0);
slp=rnc_Extract_NCEP_SLP(lon, lat, time, 0);
forcd=rnc_forcdStats(forcd,'SST');
SSTa=forcd.ano;

% step 1)  SSTa(x,y,t) --> F(x,t)

F=ConvertXYT_into_ZT(SSTa, forcd.mask);
%SSTa_back=ConvertZT_into_XYT(F, forcd.mask);

% step 2) F(x,t) = E(x,k)*P(k,t)

% E(x,k) these are the eigenvector of the covariance matrix C(x,x)
% where C(x,x) = F(x,t)*F'(t,x)
%C=F*F';
C=F*F';
[E,L]=eig(C);
l=diag(L);
l2=l.^2;
l2=l2/sum(l2)*100;

% step 3) E'(k,x)*F(x,t) = E'(k,x)*E(x,k)*P(k,t), where E'(k,x)*E(x,k)=1
P=E'*F;

% keep the first 4
l2=l2(end:-1:end-3);
E=E(:, end:-1:end-3);
P=P(end:-1:end-3,:)';
E=ConvertZT_into_XYT(E, forcd.mask);

%fourier
len=length(P);
NFFT = 2^nextpow2(len); % Next power of 2 from length of y
Y = fft(P,NFFT)/len;
%f = Fs/2*linspace(0,1,NFFT/2+1);
plot(2*abs(Y(1:NFFT/2+1)))
hold on
Y = fft(npgo,NFFT)/len;
plot(2*abs(Y(1:NFFT/2+1)),'y')
hold off


x=P(:,2);
n = length(x);
N1=n*5;
X1 = abs(fft(x,N1));
F1 = [0 : N1 - 1]/N1;
plot(F1,X1,'-x'),title('Freq'),axis([0 1 0 20])

x=npgo;
n = length(x);
N1=n;
X1 = abs(fft(x,N1));
F1 = [0 : N1 - 1]/N1;
plot(F1,X1,'-x'),title('Freq'),axis([0 1 0 20])

%compare P to NPGO
[npgo_time, npgo] = rout_npgo;
npgo=npgo(1:708);
npgo_time=npgo_time(1:708);
mfig;
plot(npgo*10,'r');
hold on
plot(P(:,2)*-1,'g'); %2nd EOF
c=corrcoef(P(:,2)*-1,npgo);hold off

%comparing E1 and E2
pcolor(E(:,:,1));
shading interp
colorbar
mfig;clf
pcolor(E(:,:,2));
shading interp
caxis([-.08 .06])
colorbar





%Computing my own F=ConvertXYT_into_ZT(SSTa, forcd.mask);
i=1;
j=1;
k=1;
ind=1;
A=zeros(31*22,708); % MANU: this should be a general statement, read from the matrix

% MANU: loops like this no good! Want to use matrix algebra instead.
for i=1:31
    for k=1:22
        for j=1:708                %need an if func for the NaN values
            A(ind,j)=SSTa(i,j,k);
            j=j+1;
        end
        j=[];
        k=k+1;
        ind=ind+1;
    end
    k=[]; i=i+1;
end
% [I,J,T]=size(SSTa);
% MANU: the statement above could be: A=reshape(SSTa, [ I*J T]);
% to go back reshape(A, [I J T]);

C=A*A';            
[E,v]=eig(C);



%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%Reconstructing PDO and NPGO from AR-1 model

lon=[-245 -110]; lat=[25 62];  %N Pacific
time=[datenum(1950,1,15) datenum(2008,12, 15)];
slp=rnc_Extract_NCEP_SLP(lon, lat, time, 0);
slp=rnc_forcdStats(slp,'SLP');
SLPa=slp.ano;
F=ConvertXYT_into_ZT(SLPa, slp.mask);
C=F*F';
[E,L]=eig(C);
l=diag(L);
l2=l.^2;
l2=l2/sum(l2)*100;
P=E'*F;
l2=l2(end:-1:end-3);
E=E(:, end:-1:end-3);
P=P(end:-1:end-3,:)';
E=ConvertZT_into_XYT(E, slp.mask);

PC1=P(:,1);  PC2=P(:,2);
%AR-1
%PDO
dt=.05;
gamma=.2;
alpha=1;
O1=0;
for i=2:length(PC1)
    O1(i)=(1-gamma*dt)*O1(i-1)+alpha*dt*PC1(i);
    i=i+1;
end

%NPGO
dt=1;
gamma=.1;
alpha=1;
O2=0;
for i=2:length(PC1)
    O2(i)=(1-gamma*dt)*O2(i-1)+alpha*dt*PC2(i);
    i=i+1;
end

[npgo_time, npgo] = rout_npgo;
pdo=rout_pdo;  mfig; plot(pdo.datenum, pdo.index)
npgo=npgo(1:708);
clf;
plot(nn(npgo))
hold on
plot(nn(O2),'c')
red_signi(npgo_time, npgo, slp.datenum, O2, 31, 1500);