clear all 

% define analytic response function f with 
%   t = time since eruption (years)
%   m = areasol mass (Tg)
f=inline('-0.5*(t).*exp(-0.6*(t)).*m/25','t','m');  

% construct linear superposition of SAT response to Sigl et al. 2015 volcanic forcing  
sigl15=load('Data/sigl15.txt');
years=sigl15(:,1);
colmass=sigl15(:,3); % unit Tg
Y1=-500; YN=2000; % first and last year of response time series
SATL=zeros(1,(YN-Y1+1)*12); % length is 2500 * 12 months
for n=1:numel(years)
  % start year of response to eruption n
  y1=years(n); 
  % last year of response to eruption n 
  %   set maximum eruption influence to 10 years
  %   truncate if out of range 
  y2=min(y1+9,YN); 
  % time index range for eruption n
  tind=(y1-Y1)*12+1:(y2-Y1)*12; 
  % superpose responses
  SATL(tind)=SATL(tind)+f([1:numel(tind)]/12,colmass(n)); 
end

% read GMST from 60 members of NoVolc and Volc 
SATR=load('Data/gmst_rcp45NoVolc_2006-2099_monthly.txt');
SATR=SATR(:,3:end);
SATV=load('Data/gmst_rcp45Volc_2006-2099_monthly.txt');
SATV=SATV(:,3:end);

% compute annual averages
month2year=[31 28 31 30 31 30 31 31 30 31 30 31]/365;
SATL=month2year*reshape(SATL,12,[]);
SATR=month2year*reshape(SATR,12,[]); SATR=reshape(SATR,[],60)';
SATV=month2year*reshape(SATV,12,[]); SATV=reshape(SATV,[],60)';

% remove anthropogenic trend by subtracting NoVolc ensemble mea
SATR_ensmean=mean(SATR);
for m=1:size(SATR,1);
  SATR(m,:)=SATR(m,:)-SATR_ensmean;
  SATV(m,:)=SATV(m,:)-SATR_ensmean;
end

% compute PDFs for annual GMST 
d=0.04;
bins=[(-4:d:4)-d/2;(-4:d:4)+d/2];
for n=1:size(bins,2)
   PDF_SATR(n)=length(find(SATR>bins(1,n) & SATR<=bins(2,n)))/numel(SATR(:));
   PDF_SATV(n)=length(find(SATV>bins(1,n) & SATV<=bins(2,n)))/numel(SATV(:));
   PDF_SATL(n)=length(find(SATL>bins(1,n) & SATL<=bins(2,n)))/numel(SATL(:)); 
end
x=mean(bins); 

% compute convolution NoVolc PDF and volcanic response PDF
PDF_CONV=zeros(size(PDF_SATR));  
NMAX=size(bins,2);
for n=1:NMAX
  for m=1:NMAX
    nm=min(NMAX,max(1,n-m+(NMAX-1)/2)+1);
    PDF_CONV(n)=PDF_CONV(n)+PDF_SATL(m)*PDF_SATR(nm);
  end
end 

% smooth curves
x_int=[-1:0.005:1];
PDF_SATR_int=interp1(x,PDF_SATR,x_int,'pchip',nan);
PDF_SATV_int=interp1(x,PDF_SATV,x_int,'pchip',nan);
PDF_SATL_int=interp1(x,PDF_SATL,x_int,'pchip',nan);
PDF_CONV_int=interp1(x,PDF_CONV,x_int,'pchip',nan);

% plot
figure(1)
clf
c2=[0.4 0.6 1];c1=[1 0.6 0.4];c12=(c1+c2)/2;
figure_height_scale=0.6;
set(gcf,'paperunits','inches','papersize',[8 6*figure_height_scale],'inverthardcopy','off','color',[1 1 1],...
        'paperposition',[0 0 8 6*figure_height_scale],'renderer','painters')
h(1)=bar(x,PDF_SATR); 
set(h(1),'facecolor',[0.8 0.8 0.8],'edgecolor','none')
ylabel('Probability density (K-1)')
xlabel('GMST anomaly after removing anthropogenic trend (K)')
hold on
h(3)=plot(x_int,PDF_SATL_int,'b','linewidth',1.5);
h(1)=plot(x_int,PDF_SATR_int,'color',c1,'linewidth',1.5);
h(2)=plot(x_int,PDF_SATV_int,'m','linewidth',1.5);
h(4)=plot(x_int,PDF_CONV_int,'g','linewidth',1.5);
hold off
legend(h,{'Model NoVolc','Model Volc','Sigl response','Convolution(NoVolc*Sigl)'},'location','northeast')
set(gca,'xlim',[-1 1],'ylim',[0 0.2],'ytick',[0 0.1 0.2],'fontsize',11,'fontweight','bold','linewidth',1);
print(gcf,'-dpng','-r300','convolution_vs_ensemble.png')

