%

clear;
close all;

status=load('eof3_1030.txt;status=load('eof2_1030.txt');');

eof3_moy=mean(eof3_1030);
eof2_moy=mean(eof2_1030);

eof3_std=std(eof3_1030);
eof2_std=std(eof2_1030);
eof3_1030=zeros(3416);
eof2_1030=zeros(3416);
for a=1:3416;
   eof3_1030(a)=(eof3_1030(a)-eof3_moy)/eof3_std;
   eof2_1030(a)=(eof2_1030(a)-eof2_moy)/eof2_std;
end;

correlations2=zeros(21,1);
correlations3=zeros(21,1);

for a=1:21;
   if a<12;
      fin=3416-(11-a);
      debut=(11-a)+1;
      x=eof2_1030(1:fin);
      y2=eof2_1030(debut:3416);
      y3=eof3_1030(debut:3416);
   else
      fin=3416-(a-11);
      debut=(a-11)+1;
      y2=eof2_1030(1:fin);
      y3=eof3_1030(1:fin);
      x=eof2_1030(debut:3416);
   end;
   clear vad;
   vad=corrcoef(x,y2);
   correlations2(a)=vad(2,1);
   clear vad;
   vad=corrcoef(x,y3);
   correlations3(a)=vad(2,1);
end;
jour=[-10:10]';
plot(jour,correlations2,'k');
hold on;
plot(jour,correlations3,'k--');
axis([-10 10 -1 1]);
y=zeros(21,1);
plot(jour,y,'k');

figure(2);
orient('landscape')

plot(eof2_1030,eof3_1030,'b.');
axis([-5 5 -5 5]);
y=zeros(11,1);
x=[-5:5]';
hold on;
plot(x,y,'k')

XCentre=0;
YCentre=0;
Rayon=1.5;
VTheta=[0:1:360]';
XCercle=zeros(360);
YCercle=zeros(360);
for a=1:360;
   XCercle(a)=XCentre+Rayon*cos((2*pi*VTheta(a))/360);
   YCercle(a)=YCentre+Rayon*sin((2*pi*VTheta(a))/360);
end
plot(XCercle, YCercle,'k');
x=[-0.00001 0]';
y=[-5 5]';
plot(x,y,'k');
x=[-5 5]';
y=[-5 5]';
plot(x,y,'k');
x=[5 -5]';
y=[-5 5]';
plot(x,y,'k');
xlabel('EOF1');
ylabel('EOF2');
title('Phase spae points for JJAS days from 1979 to 2006');
text(-4,-1,'PHASE 1');
text(-2,-4,'PHASE 2');
text(1,-4,'PHASE 3');
text(3,-1,'PHASE 4');
text(3,1,'PHASE 5');
text(1,4,'PHASE 6');
text(-2,4,'PHASE 7');
text(-4,1,'PHASE 8');
set (gca,'tick','out');

figure(1);
status=print('-depsc2','lag_correlations_eof2_eof3_1030.eps');
figure(2);
status=print('-depsc2','phase_space_diagram_eof2_eof3_1030.eps');

phase1=zeros(3416,3);
phase2=zeros(3416,3);
phase3=zeros(3416,3);
phase4=zeros(3416,3);

phase5=zeros(3416,3);
phase6=zeros(3416,3);
phase7=zeros(3416,3);
phase8=zeros(3416,3);

for a=1:3416;
   % selection hors du cercle : il faut que la distance soit >1.5
   if sqrt((eof2_1030(a))^2+(eof3_1030(a))^2) >1.5;
      if eof3_1030(a)<0  % selection phase 1 2 3 4
          % selection phase 1
          if eof3_1030(a)>eof2_1030(a);
              phase1(a,1)=1;
              phase1(a,2)=eof2_1030(a);
              phase1(a,3)=eof3_1030(a);
          % selection phase 2
          elseif eof2_1030(a)<0;
              phase2(a,1)=1;
              phase2(a,2)=eof2_1030(a);
              phase2(a,3)=eof3_1030(a);
          % selection phase 3
          elseif abs(eof2_1030(a))<abs(eof3_1030(a));
              phase3(a,1)=1;
              phase3(a,2)=eof2_1030(a);
              phase3(a,3)=eof3_1030(a);
          else
              phase4(a,1)=1;
              phase4(a,2)=eof2_1030(a);
              phase4(a,3)=eof3_1030(a);
         end;
      % selection phase 5 6 7 8
      else
         % selection phase 5
         if eof3_1030(a)<eof2_1030(a);
             phase5(a,1)=1;
             phase5(a,2)=eof2_1030(a);
             phase5(a,3)=eof3_1030(a);
         elseif eof2_1030(a)>0;
             phase6(a,1)=1;
             phase6(a,2)=eof2_1030(a);
             phase6(a,3)=eof3_1030(a);
         % selection phase 7
         elseif abs(eof2_1030(a))<abs(eof3_1030(a));
             phase7(a,1)=1;
             phase7(a,2)=eof2_1030(a);
             phase7(a,3)=eof3_1030(a);
         else
             phase8(a,1)=1;
             phase8(a,2)=eof2_1030(a);
             phase8(a,3)=eof3_1030(a);
         end;
      end;
  end;
end;

status=save('phase1.txt','phase1','-ASCII');
status=save('phase2.txt','phase2','-ASCII');
status=save('phase3.txt','phase3','-ASCII');
status=save('phase4.txt','phase4','-ASCII');

status=save('phase5.txt','phase5','-ASCII');
status=save('phase6.txt','phase6','-ASCII');
status=save('phase7.txt','phase7','-ASCII');
status=save('phase8.txt','phase8','-ASCII');
