clear
close all
clc

set_list;
%% National data
co2_n_pctg = xlsread('plot.xlsx','Fig2','C8:E8');
egy_n_pctg = xlsread('plot.xlsx','Fig2','C64:E64');
SO2_n_pctg = xlsread('plot.xlsx','Fig2','C17:E17');
NOx_n_pctg = xlsread('plot.xlsx','Fig2','C26:E26');
NH3_n_pctg = xlsread('plot.xlsx','Fig2','C35:E35');
PM25_n_pctg = xlsread('plot.xlsx','PM_total','D34:F34')/100;
mort_n_pctg = xlsread('plot.xlsx','Mortality_total','E34:G34');

color_code = [
0.5    0.5    0.5    %EGY
0      0.4470 0.7410 %SO2
0.8500 0.3250 0.0980 %NOx
0.9290 0.6940 0.1250 %NH3
0.4660 0.6740 0.1880 %PM
0.3010 0.7450 0.9330 %Mortality
];

marker_symbol = {'x', '^', 's', 'o'};
%% Regional data under 4 percent policy
[~, ~, raw] = xlsread('plot.xlsx','4 pct','C3:AD32');
data = reshape([raw{:}],size(raw));

population_2030 = data(:,1);
consumption_loss = (data(:,2) - data(:,3))./population_2030*1e9; % [USD/capita]
health_benefit = data(:,end)./population_2030*1e9; % [USD/capita]

d_CO2 = (data(:,6) - data(:,7))./data(:,6); % [percent]
d_egy = (data(:,4) - data(:,5))./data(:,4); % [percent]
d_SO2 = (data(:,8) - data(:,9))./data(:,8); % [percent]
d_NOx = (data(:,10) - data(:,11))./data(:,10); % [percent]
d_NH3 = (data(:,12) - data(:,13))./data(:,12); % [percent]
d_PM_anth = (data(:,18) - data(:,19))./data(:,18); % [percent] anth. PM
d_PM_total = (data(:,20) - data(:,21))./data(:,20); % [percent] total PM
d_mort = (data(:,22) - data(:,23))./data(:,22); % [percent]
clearvars data raw;


%% Figure 4
figure(1); clf;
set(gcf, 'unit', 'centimeter', 'pos', [1    4.5   17.8    5.75*2]);

%% Figure 4a: national data
ax1 = axes('pos', [0.015    0.6    0.36    0.36]); hold on; box on;
plot([0 60], [0 60], '--', 'color', [0.7 0.7 0.7]);
plot(co2_n_pctg*100, egy_n_pctg*100, 'color', color_code(1,:));
plot(co2_n_pctg*100, SO2_n_pctg*100, 'color', color_code(2,:));
plot(co2_n_pctg*100, NOx_n_pctg*100, 'color', color_code(3,:));
plot(co2_n_pctg*100, NH3_n_pctg*100, 'color', color_code(4,:));
plot(co2_n_pctg*100, PM25_n_pctg*100, 'color', color_code(5,:));
plot(co2_n_pctg*100, mort_n_pctg*100, 'color', color_code(6,:));

marker_size = 4;
for i = 1:3
    plot(co2_n_pctg(i)*100, egy_n_pctg(i)*100, 'marker', marker_symbol{i+1}, 'color', color_code(1,:), 'markersize', marker_size);
    plot(co2_n_pctg(i)*100, SO2_n_pctg(i)*100, 'marker', marker_symbol{i+1}, 'color', color_code(2,:), 'markersize', marker_size);
    plot(co2_n_pctg(i)*100, NOx_n_pctg(i)*100, 'marker', marker_symbol{i+1},'color', color_code(3,:), 'markersize', marker_size);
    plot(co2_n_pctg(i)*100, NH3_n_pctg(i)*100, 'marker', marker_symbol{i+1},'color', color_code(4,:), 'markersize', marker_size);
    plot(co2_n_pctg(i)*100, PM25_n_pctg(i)*100, 'marker', marker_symbol{i+1}, 'color', color_code(5,:), 'markersize', marker_size);
    plot(co2_n_pctg(i)*100, mort_n_pctg(i)*100, 'marker', marker_symbol{i+1}, 'color', color_code(6,:), 'markersize', marker_size);
end

text(co2_n_pctg(end)*100, egy_n_pctg(end)*100, '  Energy', 'color', color_code(1,:), 'fontsize', 7);
text(co2_n_pctg(end)*100, SO2_n_pctg(end)*100, '  SO_2', 'color', color_code(2,:), 'fontsize', 7);
text(co2_n_pctg(end)*100, NOx_n_pctg(end)*100-1, '  NO_x', 'color', color_code(3,:), 'fontsize', 7);
text(co2_n_pctg(end)*100, NH3_n_pctg(end)*100-0.5, '  NH_3', 'color', color_code(4,:), 'fontsize', 7);
text(co2_n_pctg(end)*100, PM25_n_pctg(end)*100, '  PM_{2.5}', 'color', color_code(5,:), 'fontsize', 7);
text(co2_n_pctg(end)*100, mort_n_pctg(end)*100+0.5, '  Mortality', 'color', color_code(6,:), 'fontsize', 7);

text(co2_n_pctg(1)*100-0.5, SO2_n_pctg(1)*100+0.5, ' \leftarrow 3% Policy', 'rotation', 90, 'fontsize', 7)
text(co2_n_pctg(2)*100-0.5, SO2_n_pctg(2)*100+0.5, ' \leftarrow 4% Policy', 'rotation', 90, 'fontsize', 7)
text(co2_n_pctg(3)*100-0.5, SO2_n_pctg(3)*100+0.5, ' \leftarrow 5%', 'rotation', 90, 'fontsize', 7)
text(co2_n_pctg(3)*100+1.5, SO2_n_pctg(3)*100-0.7, '     Policy', 'rotation', 90, 'fontsize', 7)

set(gca, 'fontsize', 6);
set(gca, 'xcolor', [1 1 1]*0.3, 'ycolor', [1 1 1]*0.3);
xlabel('CO_2 Reduction in 2030', 'color', 'k');
ylabel({'Reduction in Non-Target', 'Emissions/Quantities in 2030'}', 'color', 'k');
grid on;
axis equal;
xlim([0 50]);
ylim([0 50]);
set(gca, 'xtick', 0:10:50, 'xticklabel', {'0%', '10%', '20%', '30%', '40%', '50%'});
set(gca, 'ytick', 0:10:50, 'yticklabel', {'0%', '10%', '20%', '30%', '40%', '50%'});

text(1.5, (50-2.5), '(a)', 'fontweight', 'bold');

%% Figure 4b: consumption loss (USD per capita) vs. avoided health damage (USD per capita)
ax2 = axes('pos', [0.40    0.6000    0.58   0.3600]); hold on; box on;
x = consumption_loss;
y = health_benefit;
plot([0 700], [0 700], '--', 'color', [1 1 1]*0.7);
%plot(x,y, 'o', 'markersize', 3, 'color', [1 1 1]*0.1, 'markerf', [1 1 1]*0.1);

set(gca, 'fontsize', 6);
xlabel('Consumption Loss (USD per Capita)', 'color', 'k');
ylabel({'Avoided Health Damages (USD per Capita)'}, 'color', 'k');
grid on;
for i = 1:length(x)
    if any(ismember({'SD','JS','TJ'}, r2{i}))
        text(x(i)-1, y(i), [' ', r2{i}, ' '], 'horizontalalignment', 'right', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i), y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'FJ','LN','JX','JL','YN','XJ','HE','QH','GD','ZJ','HI','HL','SN','NX'}, r2{i}))
        text(x(i)+1, y(i), [' ', r2{i}, ' '], 'horizontalalignment', 'left', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'HB'}, r2{i}))
        text(x(i)+1, y(i)+10, [' ', r2{i}, ' '], 'horizontalalignment', 'left', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'SC'}, r2{i}))
        text(x(i)+1, y(i)-20, [' ', r2{i}, ' '], 'horizontalalignment', 'left', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'HA','SH'}, r2{i}))
        text(x(i)-20, y(i), [' ', r2{i}, ' '], 'verticalalignment', 'top', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'AH'}, r2{i}))
        text(x(i)-10, y(i), [' ', r2{i}, ' '], 'verticalalignment', 'top', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'CQ'}, r2{i}))
        text(x(i)-10, y(i), [' ', r2{i}, ' '], 'verticalalignment', 'bottom', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'HN'}, r2{i}))
        text(x(i)-10, y(i), [' ', r2{i}, ' '], 'verticalalignment', 'bottom', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'GS'}, r2{i}))
        text(x(i)-20, y(i)+5, [' ', r2{i}, ' '], 'verticalalignment', 'bottom', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    if any(ismember({'GX'}, r2{i}))
        text(x(i)-15, y(i)+10, [' ', r2{i}, ' '], 'verticalalignment', 'bottom', 'fontsize', 6, 'color', [1 1 1]*0.7);
        plot(x(i) ,y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.7, 'markerf', [1 1 1]*0.7);
    end
    % highlight provinces with full name
    if any(ismember({'SX','BJ'}, r2{i})) 
        text(x(i)-1, y(i), [' ', r2_full{i}, ' '], 'horizontalalignment', 'right', 'fontsize', 6);
        plot(x(i), y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.1, 'markerf', [1 1 1]*0.1);
    end
    if any(ismember({'NM', 'GZ'}, r2{i})) 
        text(x(i)+1, y(i), [' ', r2_full{i}, ' '], 'horizontalalignment', 'left', 'fontsize', 6);
        plot(x(i), y(i), 'o', 'markersize', 3, 'color', [1 1 1]*0.1, 'markerf', [1 1 1]*0.1);
    end
end
text(320, 345, 'One-to-One Ratio Line', 'fontsize', 7, 'color', [1 1 1]*0.7, 'rotation', 19);
text(-90, (1000-50), '(b)', 'fontweight', 'bold');


%% Figure 4c: regional data
ax3 = axes('pos', [0.08, 0.08, 0.9, 0.42]); hold on; box on;

% Provincial vertical grouping
for i = 1:30
    plot(ones(1,6)*d_CO2(i), sort([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)]), 'color', [1 1 1]*0.6);
end

% All data points
h = plot(d_CO2, d_egy, 'o', ...
         d_CO2, d_SO2, 'o', ...
         d_CO2, d_NOx, 'o', ...
         d_CO2, d_NH3, 'o', ...
         d_CO2, d_PM_anth, 'o', ...
         d_CO2, d_mort, 'o');

marker_size = 4;
set(h, 'markersize', marker_size);
set(h(1), 'markerf', color_code(1,:), 'color', color_code(1,:), 'marker', 'o');
set(h(2), 'markerf', color_code(2,:), 'color', color_code(2,:), 'marker', 'v');
set(h(3), 'markerf', color_code(3,:), 'color', color_code(3,:), 'marker', 's');
set(h(4), 'markerf', color_code(4,:), 'color', color_code(4,:), 'marker', '^');
set(h(5), 'markerf', color_code(5,:), 'color', color_code(5,:), 'marker', 'x');
set(h(6), 'markerf', color_code(6,:), 'color', color_code(6,:), 'marker', 'd');

% Linear fit
x_range = [0.04 0.525];

p = polyfit(d_CO2, d_egy, 1);
y_range = x_range.*p(1) + p(2);
h1 = plot(x_range, y_range, ':', 'color', color_code(1,:));
text(x_range(end)+0.00035, y_range(end)-0.04, 'Energy', 'fontsize', 8, 'color', color_code(1,:));

p = polyfit(d_CO2, d_SO2, 1);
y_range = x_range.*p(1) + p(2);
h2 = plot(x_range, y_range, ':', 'color', color_code(2,:));
text(x_range(end)+0.00035, y_range(end)+0.01, 'SO_2', 'fontsize', 8, 'color', color_code(2,:));

p = polyfit(d_CO2, d_NOx, 1);
y_range = x_range.*p(1) + p(2);
h3 = plot(x_range, y_range, ':', 'color', color_code(3,:));
text(x_range(end)+0.00035, y_range(end), 'NO_x', 'fontsize', 8, 'color', color_code(3,:));

p = polyfit(d_CO2, d_NH3, 1);
y_range = x_range.*p(1) + p(2);
h4 = plot(x_range, y_range-0.01, ':', 'color', color_code(4,:));
text(x_range(end)+0.00035, y_range(end), 'NH_3', 'fontsize', 8, 'color', color_code(4,:));

p = polyfit(d_CO2, d_PM_anth, 1);
y_range = x_range.*p(1) + p(2);
h5 = plot(x_range, y_range, ':', 'color', color_code(5,:));
text(x_range(end)+0.00035, y_range(end)+0.022, {'Anth.'}, 'fontsize', 8, 'color', color_code(5,:));
text(x_range(end)+0.00035, y_range(end)-0.022, {'PM_{2.5}'}, 'fontsize', 8, 'color', color_code(5,:));

p = polyfit(d_CO2, d_mort, 1);
y_range = x_range.*p(1) + p(2);
h6 = plot(x_range, y_range, ':', 'color', color_code(6,:));
text(x_range(end)+0.00035, y_range(end)+0.01, 'Mortality', 'fontsize', 8, 'color', color_code(6,:));

grid on;
set(gca, 'fontsize', 6);
xlim([0.05 0.57]);
ylim([-0.1 0.6]);
set(gca, 'xcolor', [1 1 1]*0.3);
set(gca, 'ycolor', [1 1 1]*0.3);
xlabel('CO_2 Reduction', 'color', 'k');
ylabel({'Reduction in Non-Target Emissions/Quantities'}, 'color', 'k');
set(gca, 'xtick', 0:0.05:0.55, 'xticklabel', {'0%', '5%', '10%', '15%', '20%', '25%', '30%', '35%', '40%', '45%', '50%', '55%'});
set(gca, 'ytick', -0.1:0.1:0.6, 'yticklabel', {'-10%', '0%', '10%', '20%', '30%', '40%', '50%', '60%'});

for i = 1:30
    if any(ismember({'ZJ','BJ','SN','YN','HB','NM','TJ','GZ'}, r2{i}))
        text(d_CO2(i), max([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)])+0.012, r2{i}, 'rotation', 90, 'horizontalalignment', 'left', 'fontsize', 7);
    elseif any(ismember({'CQ'}, r2{i}))
        text(d_CO2(i)+0.002, max([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)])+0.012, r2{i}, 'rotation', 90, 'horizontalalignment', 'left', 'fontsize', 7);
    elseif any(ismember({'GX'}, r2{i}))
        text(d_CO2(i)+0.005, max([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)])+0.012, r2{i}, 'rotation', 90, 'horizontalalignment', 'left', 'fontsize', 7);
    elseif any(ismember({'SC','HL'}, r2{i}))
        text(d_CO2(i)-0.003, max([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)])+0.012, r2{i}, 'rotation', 90, 'horizontalalignment', 'left', 'fontsize', 7);
    %elseif any(ismember({'HE'}, r2{i}))
    %    text(d_CO2(i)+0.001, min([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)])-0.012, r2{i}, 'rotation', 90, 'horizontalalignment', 'right', 'fontsize', 7);
    else
        text(d_CO2(i), min([d_SO2(i), d_NOx(i), d_NH3(i), d_egy(i), d_PM_anth(i), d_mort(i)])-0.012, r2{i}, 'rotation', 90, 'horizontalalignment', 'right', 'fontsize', 7);
    end
end

text(0.055, (0.6-0.04), '(c)', 'fontweight', 'bold');

%export_fig fig4_cobenefits_4pctg_tall -png -eps -r150 -painters -crop;
%export_fig fig4_cobenefits.pdf -native -crop