clear
close all
clc
format compact

set_list;
%% GDP per capita in 2010
population_2010 = xlsread('plot.xlsx','Cost-benefit','B4:B33');
gdp_r_2010 = xlsread('plot.xlsx','Cost-benefit','C4:C33');

gdp_per_capita = (gdp_r_2010*1e9)./population_2010;
[~, id_sort] = sort(gdp_per_capita, 'descend'); % Sort province by GDP per capita
r2_sort = r2(id_sort);
r2_full_sort = r2_full(id_sort);

% Consumption
[~, ~, raw] = xlsread('plot.xlsx','4 pct','D3:E32');
data = reshape([raw{:}],size(raw));
c_r_2030_bau = data(:,1);
c_r_2030_p4 = data(:,2);

% CO2
[~, ~, raw] = xlsread('plot.xlsx','4 pct','H3:I32');
data = reshape([raw{:}],size(raw));
co2_r_2030_bau = data(:,1);
co2_r_2030_p4 = data(:,2);
clearvars data raw;

% Energy
col_consumption = xlsread('plot.xlsx','4 pct', 'AG3:AG32');
gas_consumption = xlsread('plot.xlsx','4 pct', 'AH3:AH32');
oil_consumption = xlsread('plot.xlsx','4 pct', 'AI3:AI32');
nhw_consumption = xlsread('plot.xlsx','4 pct', 'AJ3:AJ32');
egy_bau2030 = transpose([col_consumption, gas_consumption, oil_consumption, nhw_consumption]);

% ==========
col_consumption = xlsread('plot.xlsx','4 pct', 'AK3:AK32');
gas_consumption = xlsread('plot.xlsx','4 pct', 'AL3:AL32');
oil_consumption = xlsread('plot.xlsx','4 pct', 'AM3:AM32');
nhw_consumption = xlsread('plot.xlsx','4 pct', 'AN3:AN32');
egy_policy2030 = transpose([col_consumption, gas_consumption, oil_consumption, nhw_consumption]);

% Anthoprogenic PM2.5
[~, ~, raw] = xlsread('plot.xlsx','4 pct','T3:U32');
data = reshape([raw{:}],size(raw));
pm_anth_r_2030_bau = data(:,1);
pm_anth_r_2030_p4 = data(:,2);
clearvars data raw;


% Totoal PM2.5 (anth. PM + other sources)
[~, ~, raw] = xlsread('plot.xlsx','4 pct','V3:W32');
data = reshape([raw{:}],size(raw));
pm_total_r_2030_bau = data(:,1);
pm_total_r_2030_p4 = data(:,2);
clearvars data raw;

% Precursor emissions
[~, ~, raw] = xlsread('plot.xlsx','4 pct','J3:S32');
data = reshape([raw{:}],size(raw));
urban_r_2030bau = data(:,1:2:9);
urban_r_2030policy = data(:,2:2:10);
d_urb = urban_r_2030policy - urban_r_2030bau;
d_urb_pctg = d_urb./urban_r_2030bau;
clearvars data raw;

%% Figure S3: BAU vs policy in 2030
w = (c_r_2030_p4 - c_r_2030_bau);               % Consumption [bn USD]
x = egy_policy2030 - egy_bau2030;               % Energy consumption [mtce]
y = (co2_r_2030_p4 - co2_r_2030_bau);           % CO2 [mega ton]
z = (pm_total_r_2030_p4 - pm_total_r_2030_bau); % total PM [um/m3]
z2 = (pm_anth_r_2030_p4 - pm_anth_r_2030_bau);  % anth PM [um/m3]

w_sort = w(id_sort);
x_sort = x(:,id_sort);
y_sort = y(id_sort);
z_sort = z(id_sort);
z2_sort = z2(id_sort);

% Percentage
wp = (c_r_2030_p4 - c_r_2030_bau)./c_r_2030_bau*100;
xp = (sum(egy_policy2030) - sum(egy_bau2030))./sum(egy_bau2030)*100;
yp = (co2_r_2030_p4 - co2_r_2030_bau)./co2_r_2030_bau*100;
zp = (pm_total_r_2030_p4 - pm_total_r_2030_bau)./pm_total_r_2030_bau*100;
z2p = (pm_anth_r_2030_p4 - pm_anth_r_2030_bau)./pm_anth_r_2030_bau*100;

wp_sort = wp(id_sort);
xp_sort = xp(id_sort);
yp_sort = yp(id_sort);
zp_sort = zp(id_sort);
z2p_sort = z2p(id_sort);

%% ====================
figure(1); clf;
set(gcf, 'units', 'centimeters', 'pos', [2.5400    3.5454   17.5   9.5]);
ht = 0.17;
x_pos = 0.08;
y_pos = 1.01;
dy = 0.2;

% ==========
ax1 = axes;
set(ax1, 'pos', [x_pos    y_pos-dy    0.90    ht]);
data = w_sort;
bar(data, 'facec', [0 0.7 0], 'edgecolor', 'none');
set(gca, 'fontsize', 6);
ylabel({'\DeltaConsumption', '(bn USD)'}, 'fontweight', 'normal');
xlim([0.3 30.7]);
set(gca, 'xtick', 1:30, 'xticklabel', []);
ylim([-30 10]);
set(gca, 'ytick', -30:10:20);
set(gca, 'ticklength', [0.0100    0.0250]/2);
my_gridline('y');

pctg = wp_sort;
for i = 1:30
    if ismember(r2_sort(i), {'GD','LN','HE','JX','SX'})
        text(i, 0, ['', num2str(pctg(i), '%2.0f'), '% '], 'color', 'k', 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
    else
        text(i, data(i), ['', num2str(pctg(i), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
    end
end
text(0.75, -24, 'National Change = -2.1%', 'fontsize', 6);
%% ==========
ax2 = axes;
set(ax2, 'pos', [x_pos    y_pos-dy*2    0.90    ht]);
data = x_sort;
hb = bar(data', 'stacked', 'edgecolor', 'none');
set(hb(1), 'facec', [0.7412         0    0.1490]);
set(hb(2), 'facec', [0.9556    0.2379    0.1464]);
set(hb(3), 'facec', [0.9935    0.6013    0.2562]);
set(hb(4), 'facec', [0.9961    0.8510    0.4627]);
set(gca, 'fontsize', 6);
ylabel({'\DeltaEnergy Use', '(mtce)'}, 'fontweight', 'normal');
xlim([0.3 30.7]);
set(gca, 'xtick', 1:30, 'xticklabel', []);
ylim([-220 0]);
set(gca, 'ytick', -200:50:0);
set(gca, 'ticklength', [0.0100    0.0250]/2);
data = sum(x_sort);
pctg = xp_sort;
for i = 1:30
    if ismember(r2_sort(i), {'SD','SX','GZ'})
        text(i, 0, ['', num2str(pctg(i), '%2.0f'), '%  '], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
    else
        text(i, data(i), ['', num2str(pctg(i), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
    end
end
text(0.75, -170, 'National Change = -21%', 'fontsize', 6);

legend('Coal', 'Gas', 'Oil', 'Non-Fossil');
set(legend, 'fontsize', 6);
set(legend, 'orientation', 'horizontal');
set(legend, 'edge', 'none');
set(legend, 'color', 'none');
set(legend, 'pos', [0.525    0.61    0.45    0.03]);
my_gridline('y');

%% ==========
ax3 = axes;
set(ax3, 'pos', [x_pos    y_pos-dy*3    0.90    ht]);
data = y_sort;
bar(data, 'facec', [0 0.7 0], 'edgecolor', 'none');
set(gca, 'fontsize', 6);
ylabel({'\DeltaCO_2 (Mt)'}, 'fontweight', 'normal');
xlim([0.3 30.7]);
set(gca, 'xtick', 1:30, 'xticklabel', []);
ylim([-500 0]);
set(gca, 'ytick', -400:100:0);
set(gca, 'ticklength', [0.0100    0.0250]/2);
my_gridline('y');
pctg = yp_sort;
for i = 1:30
    if ismember(r2_sort(i), {'NM','SX','GZ'})
        text(i, 0, ['', num2str(pctg(i), '%2.0f'), '%  '], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
    else
        text(i, data(i), ['', num2str(pctg(i), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
    end
end
text(0.75, -400, 'National Change = -24%', 'fontsize', 6);

%% ==========
ax4 = axes;
set(ax4, 'pos', [x_pos    y_pos-dy*4-0.015    0.90    ht+0.015]);
data = z_sort;
bar(data, 'facec', [0 0.7 0], 'edgecolor', 'none');
set(gca, 'fontsize', 6);
xlim([0.3 30.7]);
ylim([-30 0])
set(gca, 'ytick', -25:5:0);
ylabel({'\DeltaPM_{2.5} (\mug/m^3)'}, 'fontweight', 'normal');
set(gca, 'xtick', 1:30, 'xticklabel', r2_full_sort, 'XTickLabelRotation', 90);
set(gca, 'ticklength', [0.0100    0.0250]/2);
my_gridline('y');
pctg = zp_sort; % total PM
pctg2 = z2p_sort; % anth PM
[~, id] = sort(data, 'ascend');
for i = 1:30
    text(i, 0, ['', num2str(pctg(i), '%2.0f'), '% (', num2str(pctg2(i), '%2.f'), '%)  '], 'color', 'k', 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
end
text(0.75, -26, 'National Change = -12% (-14%)', 'fontsize', 6);

%% Figure S4: precursor emission change
figure(2); clf;
set(gcf, 'units', 'centimeters', 'pos', [2.5400    3.5454   17.5   11]);
ht = 0.14;
x_pos = 0.08;
y_pos = 0.845;
dy = 0.17;

for i = 1:5
    so2_r_2030bau = urban_r_2030bau(:,i);
    so2_r_2030policy = urban_r_2030policy(:,i);
       
    x = d_urb(:,i);
    x_sort = x(id_sort);
    
    xp = d_urb_pctg(:,i)*100;
    xp_sort = xp(id_sort);
    
    ax1 = axes;
    set(ax1, 'pos', [x_pos    y_pos-dy*(i-1)    0.90    ht]);
    bar(x_sort, 'facec', [0 0.7 0], 'edgecolor', 'none');
    set(gca, 'fontsize', 6);
    xlim([0 31]);
    set(gca, 'ticklength', [0.0100    0.0250]/2);
    
    if i==1
        ylabel('\DeltaSO_2 (Mt)', 'fontweight', 'normal');
        ylim([-1.5 0]);
    elseif i==2
        ylabel('\DeltaNO_x (Mt)', 'fontweight', 'normal');
        ylim([-1 0]);
        set(gca, 'ytick', -1:0.25:0)
    elseif i == 3
        ylabel('\DeltaNH_3 (Mt)', 'fontweight', 'normal');
        ylim([-0.07 0.02]);
    elseif i==4
        ylabel(['\DeltaBC (Mt)'], 'fontweight', 'normal');
        ylim([-0.08 0]);
        set(gca, 'ytick', -0.06:0.02:0)
    else
        ylabel(['\DeltaOC (Mt)'], 'fontweight', 'normal');
        ylim([-0.08 0]);
        set(gca, 'ytick', -0.06:0.02:0)
    end
    
    if i==5
        set(gca, 'xtick', 1:30, 'xticklabel', r2_full_sort, 'XTickLabelRotation', 90);
    else
        set(gca, 'xtick', 1:30, 'xticklabel', []);
    end
    
    [~, id] = sort(x_sort, 'ascend');
    for ir = 1:30
        if i == 1
            if ismember(ir, id([1 2])), text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'color', 'k', 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'left');
            else text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
            end
        elseif i == 2
            if ismember(ir, id([1])), text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'color', 'k', 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'left');
            else text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
            end
        elseif i == 3
            if ismember(ir, id([1 2])), text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'color', 'k', 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'left');
            else text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
            end
        else
            if ismember(ir, id(1)), text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'color', 'k', 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'left');
            else text(ir, x_sort(ir), ['', num2str(xp_sort(ir), '%2.0f'), '%'], 'fontsize', 6, 'rotation', 90, 'horizontalalignment', 'right');
            end
        end
    end
    my_gridline('y');
    
    pctg = sum(d_urb(:,i))/sum(urban_r_2030bau(:,i));
    annotation('textbox', [x_pos+0.002, y_pos-dy*(i-1)+0.025, 0.3, 0.025], 'string', ['National Change = ', num2str(pctg*100, '%2.0f'), '%'], 'fontsize', 6, 'edgecolor', 'none');
end

