%AR-1 model
%forcings: dPDO/dt=SLP_PC1(t)-PDO/tao(PDO)
%          dNPGO/dt=...

%call data
slp_pc=read_SLPs;

%AR-1: PDO
pdo=rout_pdo;
tao_PDO=5;


for i=2:length(slp_pc.pc1)
    PDO(i)=slp_pc.pc1(i-1)-PDO(i-1)/tao_PDO;
end

%AR-1: NPGO
[npgo_time, npgo] = rout_npgo;
[year_npgo,month_npgo, day_npgo]=dates_datenum(npgo_time);

Start=find(slp_pc.year==year_npgo(1) & slp_pc.month==month_npgo(1));
PC2=slp_pc.pc2(Start:end);
NPGO=npgo(1);
tao_NPGO=10;
for i=2:length(PC2)
   NPGO(i)=PC2(i-1)-NPGO(i-1)/tao_NPGO;
end


plot(NPGO)
hold on
plot(npgo,'r')
cor=corrcoef(NPGO,npgo);
cor=cor(1,2);