function [pntws,pthws] = unc_bnds(tmsrs_ens,prct)
%UNC_BNDS Computes both pointwise and pathwise confidence intervals at the 
% specified level, given an ensemble of equally likely time series.
%   [pntws,pthws] = unc_bnds(tmsrs_ens,prct) where 'tmsrs_ens(t,ens)' is an
%   ensemble of time series, 'prct' is the desired percentile (default is
%   90%), 'pntws(t,2)' contains the upper and lower bounds of the pointwise
%   uncertainty, and 'pthws(t,2)' contains the upper and lower bounds of
%   the pathwise uncertainty.
%
%   Nathan Steiger, June 2014

% DETERMINE LEVEL OF CONFIDENCE INTERVAL
if nargin == 1
  ptl=5;ptu=95;
  pct_o=0.9;
elseif nargin == 2
  ptl=(100-prct)/2;ptu=100-ptl;
  pct_o=prct/100;
end

% POINTWISE CONFIDENCE INTERVAL
pntws=prctile(tmsrs_ens,[ptl ptu],2);

% PATHWISE CONFIDENCE INTERVAL
% Find iteratively by inflating the pointwise bounds
pthws=pntws;
enst=size(tmsrs_ens,2);pct=0;

j=1;
while pct < pct_o
  % while pu < 0.9 && pl < 0.9
  
  % Increase upper and lower independently
  %if pl < 0.9
  %  pthws(:,1)=pthws(:,1)-0.005;
  %elseif pu < 0.9
  %  pthws(:,2)=pthws(:,2)+0.005;
  %end
  
  % Increase both upper and lower together
  pthws(:,1)=pthws(:,1)-0.005;
  pthws(:,2)=pthws(:,2)+0.005;
  
  upper=zeros(enst,1);
  lower=zeros(enst,1);
  for i=1:enst
    upper(i)=all(pthws(:,2)>tmsrs_ens(:,i));
    lower(i)=all(pthws(:,1)<tmsrs_ens(:,i));
  end
  
  pu=sum(upper)/enst;
  pl=sum(lower)/enst;
  
  pct=min(pu,pl);
  
  j=j+1;
  
  % Stop the iteration if no clear convergence
  if j>1000
    break
  end
  
end

end

