
% COMPUTE VS-LITE TREE ESTIMATES DURING OBSERVATIONAL PERIOD
addpath('/home/nsteiger/MATLAB/annual-da')

%==================
% LOAD PROXY DATA
%==================

load('/d1/nsteiger/proxy-data/PAGES2k_v2/v1_9_0/pages2kTSv1_9_0_unpack.mat',...
    'archive','p_lat','p_lon','proxy_ann','proxy_ama','resMin','year')
disp('Using PAGES2k data v 1.9.0')
prxy_yrs=year;
% Convert proxy lons into 0 to 360
nlns=find(p_lon<0);
p_lon(nlns)=p_lon(nlns)+360;

% standardize proxy data
proxy_ama=bsxfun(@minus,proxy_ama,nanmean(proxy_ama,1));
proxy_ama=bsxfun(@rdivide,proxy_ama,nanstd(proxy_ama,0,1));


% pull out proxy type of interest
% Just use high resolution ice cores
p_idx = find(ismember(archive,['tree'])); % 'ice', 'coral', 'documents', 'glacier ice'
lowres_id=find(resMin>1);
[~,p_ai,~]=setxor(p_idx,lowres_id);
p_idx=p_idx(p_ai);

prx_i=proxy_ama(:,p_idx);
p_lat=p_lat(p_idx);
p_lon=p_lon(p_idx);

pNum=length(p_idx);

calPeriod=[1901 2000];
%calPeriod=[1911 1995];

% Get parameter estimates and model error variance estimates
[T1,T2,M1,M2,eVar,eVarDist] = tree_vslite_calib(calPeriod,prx_i,p_lat,p_lon,...
    prxy_yrs,pNum);

%save('vs_lite_calib_params_T1upper75_calib_1911_1995.mat')

% Run VS-Lite
% trw=zeros(pNum,length(c_yrs));
% for q=1:pNum
%     [trw(q,:)] = VSLite_v2_3(1901,2000,p_lat(q),T1(q),T2(q),...
%         M1(q),M2(q),squeeze(t_cal(q,:,:)),squeeze(p_cal(q,:,:)),...
%         'lbparams',[0.76,0.01,0.093,4.886,5.80,1000,0.2,1]);
% end
% 
% save('vs_lite_trw_historical_2.mat')












