% preparing input data of reconstruction based on all proxy network

%##################################################
%%%%%%%% for index-only reconstruction %%%%%%%%%%%%
%##################################################

%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% high-f proxy for full/subproxy network
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

clear all; close all;

z=load('/holocene/s1/zuz10/work1/temann/data/rtable1209');
%---------------------------------------------------
load /holocene/s1/zuz10/work1/temann/data/itrdbmatrixh
x=itrdbmatrixh;
[n1,m1]=size(x)

%-- select all these chronologies that located within a certain area --
%                   (e.g. -120W~-60W, 20N~60N)

%%-- global proxy network---
ilon1 = -180; ilon2 = 180; 
ilat1 = -90; ilat2 = 90;   
%%--------------------------

%%-- North hemisphere proxy network---
%ilon1 = -180; ilon2 = 180; 
%ilat1 = 0; ilat2 = 90;   
%%--------------------------

%%-- South hemisphere proxy network---
%ilon1 = -180; ilon2 = 180; 
%ilat1 = -90; ilat2 = 0;   
%%--------------------------
%------------------------------------------------------------------

%----------------------------------------------------------------
%                                                               %
%  If you want a non-dendro dataset, then you just need to      %
%  delete "z(3,i)==9000/1(9000/1 is tree-ring ID) %
%                                                               %
%----------------------------------------------------------------

for istep=1:19
%istep=1;

byr=1800-(istep-1)*100;
eyr=1995;

kk=find(x(:,1)==byr);
kkk=find(x(:,1)==eyr); 

x(kk,1:3)
x(kkk,1:3)

n=0;
for i=1:m1-1 % 
%    if (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | z(3,i)==7000 |  z(3,i)==6000 |  z(3,i)==5000 |  z(3,i)==4000 |  z(3,i)==3000 |  z(3,i)==2000)  &...
    if (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | z(3,i)==7000 |  z(3,i)==6000 |  z(3,i)==5000 |  z(3,i)==4000 |  z(3,i)==3000 |  z(3,i)==2000)  &...
         x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &...
            z(1,i)>=ilon1 & z(1,i)<=ilon2 & z(2,i)>=ilat1 & z(2,i)<=ilat2 & abs(z(4,i))>=0
        n=n+1;
        y(1:kkk-kk+1,n)=x(kk:kkk,i+1);
        zz(:,n)=z(:,i);
    end
end
num(istep,1)=n;

save(strcat('/holocene/s1/zuz10/work1/temann/info/info',num2str(istep)),'zz','-ASCII');
save(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/highf/proxy',num2str(istep)),'y','-ASCII'); % no normalize again this time
clear y zz 
end
save('/holocene/s1/zuz10/work1/temann/info/numh','num','-ASCII');

%----- do the climate data and ttls parameters ----
clear all; close all;

qqq=load('/holocene/s1/zuz10/work1/temann/newtemp/nhhinfxxx'); % ---nh temp


climate=qqq; % do not need to normalize again this time
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/highf/climate','climate','-ascii')


%-- estimate ttls values -----
for istep=1:19
    
    z=load(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/highf/proxy',num2str(istep)));
    [nrows,ncols]=size(z)
    
    ppp=z(nrows-145:nrows,:);%over 1850-1995
    cli=qqq
    
    keepers=[cli ppp];
    [nrows,ncols]=size(keepers)
    means=nanmean(keepers); % 
    stdevs=nanstd(keepers);
    stdkeepers=keepers-repmat(means,nrows,1);
    stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
    
    [U,S,V] = svd(stdkeepers,0); % do the decomposition, economy style. U's are the PCs
    highfevalues=diag(S.^2);
    highfpctvar=diag(S).^2./sum(diag(S).^2); % that gives the pct for each pc
    mbest=eigenselect(highfpctvar);
    eofnumb(istep)=mbest;
end
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/highf/eofnumb','eofnumb','-ascii')
%------------------------------------------------------------------





%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% low-f proxy for full/subproxy network
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

clear all; close all;

z=load('/holocene/s1/zuz10/work1/temann/data/rtable1209');

%---------------------------------------------------
load /holocene/s1/zuz10/work1/temann/data/itrdbmatrixl
x=itrdbmatrixl;
[n1,m1]=size(x)

%-- select all these chronologies that located within a certain area --
%                   (e.g. -120W~-60W, 20N~60N)


%%-- global proxy network---
ilon1 = -180; ilon2 = 180; 
ilat1 = -90; ilat2 = 90;   
%%--------------------------

%%-- North hemisphere proxy network---
%ilon1 = -180; ilon2 = 180; 
%ilat1 = 0; ilat2 = 90;   
%%--------------------------

%%-- South hemisphere proxy network---
%ilon1 = -180; ilon2 = 180; 
%ilat1 = -90; ilat2 = 0;   
%%--------------------------
%------------------------------------------------------------------
%------------------------------------------------------------------


for istep=1:19
%istep=1;

byr=1800-(istep-1)*100;
eyr=1995;

kk=find(x(:,1)==byr);
kkk=find(x(:,1)==eyr); 

x(kk,1:3)
x(kkk,1:3)

nc=0;
for i=1:m1-1 % 
%    (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | z(3,i)==7000 |  z(3,i)==6000 |  z(3,i)==5000 |  z(3,i)==4000 |  z(3,i)==3000 |  z(3,i)==2000)  &...
    if (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | z(3,i)==7000 |  z(3,i)==6000 |  z(3,i)==5000 |  z(3,i)==4000 |  z(3,i)==3000 |  z(3,i)==2000)  &...
         x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &...
            z(1,i)>=ilon1 & z(1,i)<=ilon2 & z(2,i)>=ilat1 & z(2,i)<=ilat2 & abs(z(4,i))>=0
        nc=nc+1;
        yc(1:kkk-kk+1,nc)=x(kk:kkk,i+1);
        zzc(:,nc)=z(:,i);
    end
end
nd=0; yd=[]; zzd=[];
for i=1:m1-1 % 
%    if (z(3,i)==9001 |  z(3,i)==8001 |  z(3,i)==7501 | z(3,i)==7001 |  z(3,i)==6001 |  z(3,i)==5001 |  z(3,i)==4001 |  z(3,i)==3001)  &...%no Lut data
    if (z(3,i)==9001 |  z(3,i)==8001 |  z(3,i)==7501 | z(3,i)==7001 |  z(3,i)==6001 |  z(3,i)==5001 |  z(3,i)==4001 |  z(3,i)==3001 |  z(3,i)==2001)  &...
         x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &...
            z(1,i)>=ilon1 & z(1,i)<=ilon2 & z(2,i)>=ilat1 & z(2,i)<=ilat2 & abs(z(5,i))>=0
        nd=nd+1;
        yd(1:kkk-kk+1,nd)=x(kk:kkk,i+1);
        zzd(:,nd)=z(:,i);
    end
end


%num(istep,1)=nc+nd;
num(istep,1)=nc;
num(istep,2)=nd;
y=[yc yd];
zz=[zzc zzd];

save(strcat('/holocene/s1/zuz10/work1/temann/info/info',num2str(istep)),'zz','-ASCII');
save(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/lowf/proxy',num2str(istep)),'y','-ASCII'); % no normalize again this time
clear y zz yc yd zzc zzd
end
save('/holocene/s1/zuz10/work1/temann/info/numl','num','-ASCII');


%----- do the climate data and ttls parameters ----
clear all; close all;

qqq=load('/holocene/s1/zuz10/work1/temann/newtemp/nhlinfxxx');  % ---nh temp

climate=qqq; % do not need to normalize again this time
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/lowf/climate','climate','-ascii')


%-- estimate ttls values -----
for istep=1:19
    
    z=load(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/lowf/proxy',num2str(istep)));
    [nrows,ncols]=size(z)
    
    ppp=z(nrows-145:nrows,:);%over 1850-1995
    cli=qqq
    
    keepers=[cli ppp];
    [nrows,ncols]=size(keepers)
    means=nanmean(keepers); % 
    stdevs=nanstd(keepers);
    stdkeepers=keepers-repmat(means,nrows,1);
    stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
    
    [U,S,V] = svd(stdkeepers,0); % do the decomposition, economy style. U's are the PCs
    lowfevalues=diag(S.^2);
    lowfpctvar=diag(S).^2./sum(diag(S).^2); % that gives the pct for each pc
    cumlowfpctvar=cumsum(lowfpctvar);
    ttlslowfpar=length(find(cumlowfpctvar<0.5))+1;
    if ttlslowfpar==0
    	ttlslowfpar=1;
    end
    eofnumb(istep)=ttlslowfpar;
end
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhglfulihad/lowf/eofnumb','eofnumb','-ascii')
%------------------------------------------------------------------
%------------------------------------------------------------------


