% to select all variables independence of r criteria we need to switch line 224/225

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   revised file from propick.m. This file is suitable for          %%
%% all climate and proxy matrix with first three rows are locations  %%
%% and ID(not important), and the first column is the year for the   %%
%% data. the climate sensitive proxies will be put in a matrix       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

X_1=load('/holocene/s1/zuz10/work1/temann/data/clidatah');
%-- change X_1 to the climate data format in the propick.m --
[i1,j1]=size(X_1);
X(1:2012,1:j1)=NaN;
X(1:3,1:j1)=X_1(1:3,1:j1); % The first three rows: locations and ID
for i=1:2009
    X(i+3,1)=i-1; % the first column: year
end
if X_1(4,1)<0
    X(4:i1-abs(X_1(4,1)),1:j1)=X_1(4+abs(X_1(4,1)):i1,1:j1);
else
    X(4+X_1(4,1):i1+X_1(4,1),1:j1)=X_1(4:i1,1:j1);
end
[im,in]=size(X)
%++++++++++ modify in 8/3/2004 for normalizing the data ++++++++++++++

%Xn=X(4:2012,2:j1);
%[Xn]=normf(Xn);
%X(4:2012,2:j1)=Xn;

%++++++++++ end of modify in 8/3/2004 for normalizing the data +++++++

%-------------------------------------------------------------

%load D:\scottzhrec1\data\itrdbmatrix % proxy dataset
load /holocene/s1/zuz10/work1/temann/data/itrdbmatrixh % proxy dataset

itrdbmatrix=itrdbmatrixh;
%-- change itrdbmatrix to the climate data format in the propick.m --
[i1,j1]=size(itrdbmatrix);
proxyfile(1:2012,1:j1)=NaN;
proxyfile(1:3,1:j1)=itrdbmatrix(1:3,1:j1); % The first three rows: locations and ID
for i=1:2009
    proxyfile(i+3,1)=i-1; % the first column: year
end
if itrdbmatrix(4,1)<0
    proxyfile(4:i1-abs(itrdbmatrix(4,1)),1:j1)=itrdbmatrix(4+abs(itrdbmatrix(4,1)):i1,1:j1);
else
    proxyfile(4+itrdbmatrix(4,1):i1+itrdbmatrix(4,1),1:j1)=itrdbmatrix(4:i1,1:j1);
end

Y=proxyfile;
[pm,pn]=size(Y)
  % ==============           end of options        ========================
byr = 1850; % beginning year to calculate r
eyr = 1995; % end year to calculate r
%byr = 1900; % 
%eyr = 1995; % late corr100
  
  
k1=find(X(:,1)==byr)
k2=find(X(:,1)==eyr)

igrid = 2; %numbers of gridpoints near proxy site used to cal. r
igr = 1; % numbers of gripoints need to pass the CL for selecting a proxy series


n4=0; d4=[]; % n4--> number of selected proxy series, d4-->column of selected proxy matrix

for i=1:pn-1
    n1=0; d1=[];
    for j=1:in-1
        n1=n1+1;
%----------------------------great circle distance------------------------------------------        
        eqrad = 6378; %km (3963 mi) Equatorial radius (surface to center distance)
%        b = 6357; %km (3950 mi) Polar radius (surface to center distance)
        ecc = 0.081082; %Eccentricity
        
%                 lon1=proxylonglat(i,1)*pi/180; % longitude of the proxy
%                 lat1=proxylonglat(i,2)*pi/180; % latitude of the proxy
%                 lon2=keeplonglat(:,1)*pi/180; % longitude of all instrumental gridboxes
%                 lat2=keeplonglat(:,2)*pi/180; % latitude of all instrumental gridboxes
        
                 lon1=X(1,j+1)*pi/180; % longitude of the proxy
                 lat1=X(2,j+1)*pi/180; % latitude of the proxy
                 lon2=Y(1,i+1)*pi/180; % longitude of all instrumental gridboxes
                 lat2=Y(2,i+1)*pi/180; % latitude of all instrumental gridboxes
                 rad = eqrad * sqrt(1 - ecc^2) / (1 - ecc^2 * (sin(lat1))^2);%radius of curvature
                 dlon = lon2 - lon1;
                 dlat = lat2 - lat1;
                 one = (sin(dlat/2)).^2 + cos(lat1).*cos(lat2).*(sin(dlon/2)).^2;
                 ccccc = 2 * atan2( sqrt(one),sqrt(1-one) );
%                 d = rad.*ccccc;
                 d1(n1) = rad.*ccccc;
%                 d(:,2)=1:2008; % just add row number for each d
%                 sorteddistance=sortrows(d,1); % sort the distance from the proxy to the center of each gridbox
%                 closest2columns=sorteddistance(1:ngridsearch,2); % that would be the  nearest gridboxes

%----------------------------great circle distance------------------------------------------        
        
        
%        d1(n1)=(Y(1,i+1)-X(1,j+1))^2+(Y(2,i+1)-X(2,j+1))^2;
    end
    [B,isort] = sort(d1); % sort the min distant gridpoints
    
    n3=0; d3=[];
    for k=1:igrid % %numbers of gridpoints near proxy site used to cal. corrcoef
        n2=0; d2=[];
        for kk=1:k2-k1+1
            if Y(kk+k1-1,i+1)>-9999 & X(kk+k1-1,isort(k)+1)>-9999
                n2=n2+1; 
                d2(n2,1)=Y(kk+k1-1,i+1);
                d2(n2,2)=X(kk+k1-1,isort(k)+1);
            end
        end
        
%        if n2>3 % at least 4 pair of records
        if n2>20 % at least 20 pair of records % modified 08/08/04
            n3=n3+1;
            rr=corrcoef(d2);
            aa(k,1)=(rr(1,2));
%            k3=find(rtable1(:,1)==n2);
%            k3=find(rtable1(:,1)==n2-fix(n2/20)); %04//19/05 degrees of freedom
%            k4=find(rtable1(1,:)==CL);
%            if abs(rr(1,2))>rtable1(k3,k4)
            if abs(rr(1,2))>=0.0 % select entire proxies
%            if abs(rr(1,2))>0.14 % mike's others
%--------------------------------------------------------------
%            if (rr(1,2))>0.14 % mike's tree +, (one sided, N=140, 95%)
%            if (rr(1,2))>0.109 % mike's tree +, (one sided, N=140, 90%)
%--------------------------------------------------------------                
                d3(n3)=1;
            else
                d3(n3)=0;
            end
        end
    end
    
    if sum(d3)>=igr % at least 'igr' correlation of 'igrid' pass the CL
        n4=n4+1;
        
        bb=abs(aa);
        ra=max(bb);
        [iii,jjj]=find(bb==ra);
        if sum(iii)>2
            rra(n4)=aa(1,1); %07/20/06
        else
            rra(n4)=aa(iii,jjj);
        end

        d4(n4)=i+1;
    end
    
    
end

d4
n4

%==================================
% to form matrix for reconstruction
%==================================

%XX(1:im,1:in+n4) = NaN;
%XX(:,1:in) = X;

%for i=1:n4
%    XX(:,in+i) = Y(:,d4(i)); % this matrix includes first column of years
%end

for i=1:n4
%    PP(:,i) = Y(:,d4(i)); % this proxy network 07/20/06
    PPP(1:3,i)=Y(1:3,d4(i));
    PPP(4,i)=rra(i);
end


save('/holocene/s1/zuz10/work1/temann/data/infoh','PPP','-ascii')




% to select all variables independence of r criteria we need to switch line 224/225

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   revised file from propick.m. This file is suitable for          %%
%% all climate and proxy matrix with first three rows are locations  %%
%% and ID(not important), and the first column is the year for the   %%
%% data. the climate sensitive proxies will be put in a matrix       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

X_1=load('/holocene/s1/zuz10/work1/temann/data/clidatah');
%-- change X_1 to the climate data format in the propick.m --
[i1,j1]=size(X_1);
X(1:2012,1:j1)=NaN;
X(1:3,1:j1)=X_1(1:3,1:j1); % The first three rows: locations and ID
for i=1:2009
    X(i+3,1)=i-1; % the first column: year
end
if X_1(4,1)<0
    X(4:i1-abs(X_1(4,1)),1:j1)=X_1(4+abs(X_1(4,1)):i1,1:j1);
else
    X(4+X_1(4,1):i1+X_1(4,1),1:j1)=X_1(4:i1,1:j1);
end
[im,in]=size(X)
%++++++++++ modify in 8/3/2004 for normalizing the data ++++++++++++++

%Xn=X(4:2012,2:j1);
%[Xn]=normf(Xn);
%X(4:2012,2:j1)=Xn;

%++++++++++ end of modify in 8/3/2004 for normalizing the data +++++++

%-------------------------------------------------------------

%load D:\scottzhrec1\data\itrdbmatrix % proxy dataset
load /holocene/s1/zuz10/work1/temann/data/itrdbmatrixh % proxy dataset

itrdbmatrix=itrdbmatrixh;
%-- change itrdbmatrix to the climate data format in the propick.m --
[i1,j1]=size(itrdbmatrix);
proxyfile(1:2012,1:j1)=NaN;
proxyfile(1:3,1:j1)=itrdbmatrix(1:3,1:j1); % The first three rows: locations and ID
for i=1:2009
    proxyfile(i+3,1)=i-1; % the first column: year
end
if itrdbmatrix(4,1)<0
    proxyfile(4:i1-abs(itrdbmatrix(4,1)),1:j1)=itrdbmatrix(4+abs(itrdbmatrix(4,1)):i1,1:j1);
else
    proxyfile(4+itrdbmatrix(4,1):i1+itrdbmatrix(4,1),1:j1)=itrdbmatrix(4:i1,1:j1);
end

Y=proxyfile;
[pm,pn]=size(Y)
  % ==============           end of options        ========================
%byr = 1850; % beginning year to calculate r
%eyr = 1995; % end year to calculate r
byr = 1896; % 
eyr = 1995; % late corr100
  
  
k1=find(X(:,1)==byr)
k2=find(X(:,1)==eyr)

igrid = 2; %numbers of gridpoints near proxy site used to cal. r
igr = 1; % numbers of gripoints need to pass the CL for selecting a proxy series


n4=0; d4=[]; % n4--> number of selected proxy series, d4-->column of selected proxy matrix

for i=1:pn-1
    n1=0; d1=[];
    for j=1:in-1
        n1=n1+1;
%----------------------------great circle distance------------------------------------------        
        eqrad = 6378; %km (3963 mi) Equatorial radius (surface to center distance)
%        b = 6357; %km (3950 mi) Polar radius (surface to center distance)
        ecc = 0.081082; %Eccentricity
        
%                 lon1=proxylonglat(i,1)*pi/180; % longitude of the proxy
%                 lat1=proxylonglat(i,2)*pi/180; % latitude of the proxy
%                 lon2=keeplonglat(:,1)*pi/180; % longitude of all instrumental gridboxes
%                 lat2=keeplonglat(:,2)*pi/180; % latitude of all instrumental gridboxes
        
                 lon1=X(1,j+1)*pi/180; % longitude of the proxy
                 lat1=X(2,j+1)*pi/180; % latitude of the proxy
                 lon2=Y(1,i+1)*pi/180; % longitude of all instrumental gridboxes
                 lat2=Y(2,i+1)*pi/180; % latitude of all instrumental gridboxes
                 rad = eqrad * sqrt(1 - ecc^2) / (1 - ecc^2 * (sin(lat1))^2);%radius of curvature
                 dlon = lon2 - lon1;
                 dlat = lat2 - lat1;
                 one = (sin(dlat/2)).^2 + cos(lat1).*cos(lat2).*(sin(dlon/2)).^2;
                 ccccc = 2 * atan2( sqrt(one),sqrt(1-one) );
%                 d = rad.*ccccc;
                 d1(n1) = rad.*ccccc;
%                 d(:,2)=1:2008; % just add row number for each d
%                 sorteddistance=sortrows(d,1); % sort the distance from the proxy to the center of each gridbox
%                 closest2columns=sorteddistance(1:ngridsearch,2); % that would be the  nearest gridboxes

%----------------------------great circle distance------------------------------------------        
        
        
%        d1(n1)=(Y(1,i+1)-X(1,j+1))^2+(Y(2,i+1)-X(2,j+1))^2;
    end
    [B,isort] = sort(d1); % sort the min distant gridpoints
    
    n3=0; d3=[];
    for k=1:igrid % %numbers of gridpoints near proxy site used to cal. corrcoef
        n2=0; d2=[];
        for kk=1:k2-k1+1
            if Y(kk+k1-1,i+1)>-9999 & X(kk+k1-1,isort(k)+1)>-9999
                n2=n2+1; 
                d2(n2,1)=Y(kk+k1-1,i+1);
                d2(n2,2)=X(kk+k1-1,isort(k)+1);
            end
        end
        
%        if n2>3 % at least 4 pair of records
        if n2>20 % at least 20 pair of records % modified 08/08/04
            n3=n3+1;
            rr=corrcoef(d2);
            aa(k,1)=(rr(1,2));
%            k3=find(rtable1(:,1)==n2);
%            k3=find(rtable1(:,1)==n2-fix(n2/20)); %04//19/05 degrees of freedom
%            k4=find(rtable1(1,:)==CL);
%            if abs(rr(1,2))>rtable1(k3,k4)
            if abs(rr(1,2))>=0.0 % select entire proxies
%            if abs(rr(1,2))>0.14 % mike's others
%--------------------------------------------------------------
%            if (rr(1,2))>0.14 % mike's tree +, (one sided, N=140, 95%)
%            if (rr(1,2))>0.109 % mike's tree +, (one sided, N=140, 90%)
%--------------------------------------------------------------                
                d3(n3)=1;
            else
                d3(n3)=0;
            end
        end
    end
    
    if sum(d3)>=igr % at least 'igr' correlation of 'igrid' pass the CL
        n4=n4+1;
        
        bb=abs(aa);
        ra=max(bb);
        [iii,jjj]=find(bb==ra);
        if sum(iii)>2
            rra(n4)=aa(1,1); %07/20/06
        else
            rra(n4)=aa(iii,jjj);
        end

        d4(n4)=i+1;
    end
    
    
end

d4
n4

%==================================
% to form matrix for reconstruction
%==================================

%XX(1:im,1:in+n4) = NaN;
%XX(:,1:in) = X;

%for i=1:n4
%    XX(:,in+i) = Y(:,d4(i)); % this matrix includes first column of years
%end

for i=1:n4
%    PP(:,i) = Y(:,d4(i)); % this proxy network 07/20/06
    PPP(1:3,i)=Y(1:3,d4(i));
    PPP(4,i)=rra(i);
end

save('/holocene/s1/zuz10/work1/temann/data/infohlater96','PPP','-ascii')
