% clc;
% clear;
% close all;
%%
% 读取xlsx文件
filename = fullfile("data", "Togtoh.xlsx"); sheetName = '2025.1.16';
T_ori = readtable(filename,'Sheet', sheetName, 'VariableNamingRule', 'preserve');
T = T_ori(T_ori.plume_quality>=0.7, :);
% 提取日期和排放量及不确定度
dates = datetime(T{:,'date'}, 'InputFormat', 'yyyyMMdd');
dates_num = datenum(dates);
[dates_num, index_date] = unique(dates_num);
emissions = T.emission;
emissions = emissions(index_date);
emission_uncertainty = T.emission_uncertainty;
emission_uncertainty = emission_uncertainty(index_date);
satellite = string(T.platform);
satellite = satellite(index_date);
%%
idx1 = strcmp(satellite, 'GF5A');
idx2 = strcmp(satellite, 'GF5B');
idx3 = strcmp(satellite, 'ZY1E');
idx4 = strcmp(satellite, 'ZY1F');
% bias1= -3.48*0.01;
% bias2 = 4.87*0.01;
% emissions(idx1 | idx2) = emissions(idx1 | idx2) / (1 + bias1);
% emissions(idx3 | idx4) = emissions(idx3 | idx4) / (1 + bias2);
%%
fig = figure('Position', [200 100 1500 400]);
hold on;
yyaxis right; % 激活右侧纵轴
%% 绘制检修时间
% % 自定义颜色映射
% colorMap = containers.Map({'A', 'B', 'C'}, {'#B9D4E7', '#D0E2EF', '#E9F0F8'});
% 
% % 读取检修时间表
% maintenanceSheet = readtable(filename, 'Sheet', '检修时间表', 'VariableNamingRule', 'preserve');
% 
% % 将起始时间和结束时间转换为datetime类型
% Time.start = datetime(num2str(maintenanceSheet{:,'起始时间'}), 'InputFormat', 'yyyyMMdd');
% Time.end = datetime(num2str(maintenanceSheet{:,'结束时间'}), 'InputFormat', 'yyyyMMdd');
% 
% % 绘制检修时间区间色块
% for i = 1:height(maintenanceSheet)
%     startDate = Time.start(i);
%     endDate = Time.end(i);
%     level = maintenanceSheet{i,"级别"};
%     color = colorMap(char(level));
% 
%     % 计算矩形的位置和宽度
%     x = startDate;
%     y = 8500; % 假设y坐标从20开始
%     width = days(endDate - startDate); % 将宽度转换为天数
%     height = 480;
% 
%     % 转换x和width为datenum以便rectangle函数使用
%     x_datenum = datenum(x);
%     width_datenum = datenum(startDate + width) - x_datenum;
% 
%     rectangle('Position', [x_datenum, y, width_datenum, height], 'FaceColor', color, 'EdgeColor', 'none');
% end

%% 绘制年排放量
% 计算年排放量
idx1_2 = dates_num >= datenum(datetime("2021-07-01")) & dates_num <= datenum(datetime("2021-12-31"));
idx2_1 = dates_num >= datenum(datetime("2022-01-01")) & dates_num <= datenum(datetime("2022-06-30"));
idx2_2 = dates_num >= datenum(datetime("2022-07-01")) & dates_num <= datenum(datetime("2022-12-31"));
idx3_1 = dates_num >= datenum(datetime("2023-01-01")) & dates_num <= datenum(datetime("2023-06-30"));
idx3_2 = dates_num >= datenum(datetime("2023-07-01")) & dates_num <= datenum(datetime("2023-12-31"));
idx4_1 = dates_num >= datenum(datetime("2024-01-01")) & dates_num <= datenum(datetime("2024-06-30"));
idx4_2 = dates_num >= datenum(datetime("2024-07-01")) & dates_num <= datenum(datetime("2024-12-31"));
emissions_1 = mean(emissions(idx1_2));
emissions_2 = mean(emissions(idx2_1));
emissions_3 = mean(emissions(idx2_2));
emissions_4 = mean(emissions(idx3_1));
emissions_5 = mean(emissions(idx3_2));
emissions_6 = mean(emissions(idx4_1));
emissions_7 = mean(emissions(idx4_2));

emission_uncertainty_2021 = sqrt(sum(emission_uncertainty(idx1_2).^2) / sum(idx1_2)) .* 365 * 24;
emission_uncertainty_2022 = sqrt(sum(emission_uncertainty(idx2_1|idx2_2).^2) / sum(idx2_1|idx2_2)) .* 365 * 24;
emission_uncertainty_2023 = sqrt(sum(emission_uncertainty(idx3_1|idx3_2).^2) / sum(idx3_1|idx3_2)) .* 365 * 24;
emission_uncertainty_2024 = sqrt(sum(emission_uncertainty(idx4_1|idx4_2).^2) / sum(idx4_1|idx4_2)) .* 365 * 24;
%%
yyaxis left; % 激活左侧纵轴
days = [184; 181; 184; 181; 184; 182; 184];
heights = [emissions_1; emissions_2; emissions_3; emissions_4; emissions_5; emissions_6; emissions_7] .* days * 24; % 每个柱子的高度
start_points = datenum(datetime({'2021-07-03'; '2022-01-03'; '2022-07-03'; '2023-01-03'; '2023-07-03'; '2024-01-03'; '2024-07-03'})); % 每个柱子的起始位置
end_points = datenum(datetime({'2021-12-29'; '2022-06-28'; '2022-12-29'; '2023-06-28'; '2023-12-29'; '2024-06-28'; '2024-12-29'})); % 每个柱子的起始位置
width = 10; % 设置柱子的宽度（天数）

% 遍历数据点并为每一个绘制一个矩形
hold on;
for i = 1:length(heights)
    % 计算每个矩形的位置（左下角）
    pos = [start_points(i), 0, end_points(i) - start_points(i), heights(i)];

    % 使用 rectangle 函数绘制矩形
    rectangle('Position', pos, 'FaceColor', '#E9F0F8');
end
ylim([0 1.6e7])
hold on;
ylabel('CO_2 Emission (Mt/half year)');
text(datenum(datetime('2021-08-01')), heights(1) + 1e6, ['2021H2: ', num2str(round(heights(1)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
text(datenum(datetime('2022-02-01')), heights(2) + 1e6, ['2022H1: ', num2str(round(heights(2)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
text(datenum(datetime('2022-08-01')), heights(3) + 1e6, ['2022H2: ', num2str(round(heights(3)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
text(datenum(datetime('2023-02-01')), heights(4) + 1e6, ['2023H1: ', num2str(round(heights(4)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
text(datenum(datetime('2023-08-01')), heights(5) + 1e6, ['2023H2: ', num2str(round(heights(5)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
text(datenum(datetime('2024-02-01')), heights(6) + 1e6, ['2024H1: ', num2str(round(heights(6)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
text(datenum(datetime('2024-08-01')), heights(7) + 1e6, ['2024H2: ', num2str(round(heights(7)/1e6, 2)), 'Mt'], 'FontSize', 10, 'FontWeight', 'bold');
%%
yyaxis right; % 激活左侧纵轴
color1 = '#C0321A';
color2 = '#DD7C4F';
color3 = '#547BB4';
color4 = '#6C61AF';
errorbar(dates_num, emissions, emission_uncertainty, '--', 'MarkerFaceColor', 'k', ...
    'LineWidth', 0.5, 'Color', 'k'); % 误差条颜色为深灰，线条粗细为2
h1 = errorbar(dates_num(idx1), emissions(idx1), emission_uncertainty(idx1), 'o', 'MarkerFaceColor', color1, ...
    'LineWidth', 1.6, 'Color', color1); % 误差条颜色为深灰，线条粗细为2
h2 = errorbar(dates_num(idx2), emissions(idx2), emission_uncertainty(idx2), 'o', 'MarkerFaceColor', color2, ...
    'LineWidth', 1.6, 'Color', color2); % 误差条颜色为深灰，线条粗细为2
h3 = errorbar(dates_num(idx3), emissions(idx3), emission_uncertainty(idx3), 'o', 'MarkerFaceColor', color3, ...
    'LineWidth', 1.6, 'Color', color3); % 误差条颜色为深灰，线条粗细为2
h4 = errorbar(dates_num(idx4), emissions(idx4), emission_uncertainty(idx4), 'o', 'MarkerFaceColor', color4, ...
    'LineWidth', 1.6, 'Color', color4); % 误差条颜色为深灰，线条粗细为2
% 设置图形属性
ylabel('CO_2 Emission (t/h)');

% % 定义固定的横轴日期
fixedDates = datetime({'2021-07-01', '2021-10-01', ...
    '2022-01-01', '2022-04-01', '2022-07-01', '2022-10-01', ...
    '2023-01-01', '2023-04-01', '2023-07-01', '2023-10-01', ...
    '2024-01-01', '2024-04-01', '2024-07-01', '2024-10-01'});
fixedDates_label = {'21-07-01', '21-10-01', ...
    '22-01-01', '22-04-01', '22-07-01', '22-10-01', ...
    '23-01-01', '23-04-01', '23-07-01', '23-10-01', ...
    '24-01-01', '24-04-01', '24-07-01', '24-10-01'};
fixedDates_num = datenum(fixedDates);

%% 标签
lgd = legend({'';'GF5A';'GF5B';'ZY1E';'ZY1F'});
% 设置图例位置为右下角
lgd.Location = 'southoutside';
% 使图例外部适应绘图区域
lgd.Box = 'off';
% 横向排列图例项
lgd.Orientation = 'horizontal';

%% 坐标轴设置
% 设置x轴刻度位置和标签
xticks(fixedDates_num);
xticklabels(fixedDates_label);

xlim([datenum(datetime('2021-07-01')), datenum(datetime('2024-12-31'))])
% 调整y轴限制以显示色块
ylim([0, 6000]);

grid on;
box on;
% 修改字体为 Arial，大小为 16
set(gca, 'FontName', 'Arial', 'FontSize', 12);
set(lgd, 'FontName', 'Arial', 'FontSize', 10); % 确保图例也使用相同的字体

hAx = gca; % 获取当前坐标轴句柄
% 设置左轴颜色为黑色
hAx.YAxis(1).Color = 'black'; 
% 设置右轴颜色为黑色
hAx.YAxis(2).Color = 'black';
% 设置x轴颜色为黑色
hAx.XColor = 'black';

% 获取左侧y轴现有的刻度位置
leftYTicks = hAx.YAxis(1).TickValues;

% 获取右侧y轴现有的刻度位置
rightYTicks = hAx.YAxis(2).TickValues;
% 手动定义左侧y轴的刻度标签文本
leftLabels = compose('%d', leftYTicks/1e6); % 使用compose函数生成标签文本，保留一位小数

% 手动定义右侧y轴的刻度标签文本
rightLabels = compose('%d', rightYTicks); % 例如，将右侧y轴标签设置为整数格式
% 设置左侧y轴的刻度标签
yyaxis left
set(hAx, 'YTickLabel', leftLabels);

% 设置右侧y轴的刻度标签
yyaxis right
set(hAx, 'YTickLabel', rightLabels);
hold off;

%% 导出图片
figname = ['C:\Users\huayi\Desktop\NC论文\Figure2\大唐托克托全年发电量\', 'Togtoh统计结果.png'];
print(fig, figname, '-dpng', '-r600');

%% 绘制上网电量与年排放量的散点图
fig2 = figure('Position', [200 200 800 600]);
hold on 
x = 2021:2024;

y1 = [24.622; 25.448; 26.631; 27.703]; % 上网电量
y2 = [emissions_1; (emissions_2 + emissions_3)/2; (emissions_4 + emissions_5)/2; (emissions_6 + emissions_7)/2] * 365 *24 / 1e6;
% 计算排放比例
ratio = y1 ./ y2;
power_coal_ratio = mean(ratio(1:3));
thermal_power_2024 = power_coal_ratio * y2(4); % 10^9 kWh
green_power_2024 = y1(4) - thermal_power_2024; % 10^9 kWh
reduced_CO2 = green_power_2024 / power_coal_ratio; % Mt CO2
reduced_coal = 143 / 350 * reduced_CO2; % Mt Coal
fprintf("2024年火电：%.2f（10^9kWh）; \n绿电：%.2f（10^9kWh）\n", thermal_power_2024, green_power_2024);
fprintf("节约标煤：%.2f（Mt）; \n减少碳排放：%.2f（Mt）\n", reduced_coal, reduced_CO2);
y1_a = [y1(1) 0;
    y1(2) 0;
    y1(3) 0;
    thermal_power_2024 green_power_2024];
%
yyaxis left;
b = bar(x - 0.2, y1_a, 0.35, 'stacked'); 
set(b(1), 'FaceColor', '#DD7C4F', 'EdgeColor', 'k', 'LineWidth', 1.2); % 设置火电颜色
set(b(2), 'FaceColor', '#548B64', 'EdgeColor', 'k', 'LineWidth', 1.2); % 设置绿电颜色
%
yyaxis right;
% 创建柱状图
b = bar(x + 0.2, y2, 0.35, 'FaceColor', '#6C61AF', 'LineWidth', 1.2);
% 获取当前轴
ax = gca;
% 获取柱状图对象的XData，这是每个柱子的基底位置
xdata = b.XData;
% 计算每个柱子的中心位置，因为bar的宽度是0.35，所以需要调整x坐标以指向柱子中心
x_centers = xdata;
err = [emission_uncertainty_2021; emission_uncertainty_2022; emission_uncertainty_2023; emission_uncertainty_2024] / 1e6;
% 添加误差线，这里假设 err 是一个和 y2 相同大小的数组，表示每个柱子的误差
hold on;
errorbar(ax, x_centers, y2, err, 'k.', 'LineStyle', 'none', 'LineWidth', 1.2);

xlim([2020.5 2024.5])
xticks(2021:2024)

hAx = gca; % 获取当前坐标轴句柄
% 设置左轴颜色为黑色
hAx.YAxis(1).Color = 'black'; 
% 设置右轴颜色为黑色
hAx.YAxis(2).Color = 'black';
% 设置x轴颜色为黑色
hAx.XColor = 'black';

yyaxis left;
% 添加数值标注
for i = 1:length(x)
    % 对于上网电量数据
    text(x(i) - 0.18, sum(y1_a(i,:)) + 0.2, sprintf('%.1f', sum(y1_a(i,:))), ...
         'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', ...
         'FontName', 'Arial', 'FontSize', 12, 'Color', 'k', 'FontWeight','bold');
end

% 标注2024年堆叠柱上的百分比
year_2024_index = length(x); % 假设x是按年排序的，那么最后一个元素是2024年
thermal_power_2024 = y1_a(year_2024_index, 1);
green_power_2024 = y1_a(year_2024_index, 2);
total_power_2024 = thermal_power_2024 + green_power_2024;

% 计算百分比
thermal_percentage = (thermal_power_2024 / total_power_2024) * 100;
green_percentage = (green_power_2024 / total_power_2024) * 100;


% 在火电和绿电的中间位置标注百分比
text(x(year_2024_index) - 0.18, thermal_power_2024 / 2, sprintf('%.1f%%', thermal_percentage), ...
     'VerticalAlignment', 'middle', 'HorizontalAlignment', 'center', ...
     'FontName', 'Arial', 'FontSize', 10, 'Color', 'w', 'FontWeight','bold'); % 使用白色以确保在深色背景上可见

text(x(year_2024_index) - 0.18, thermal_power_2024 + green_power_2024 / 2, sprintf('%.1f%%', green_percentage), ...
     'VerticalAlignment', 'middle', 'HorizontalAlignment', 'center', ...
     'FontName', 'Arial', 'FontSize', 10, 'Color', 'w',  'FontWeight','bold'); % 使用白色以确保在深色背景上可见

yyaxis right;
% 添加数值标注
for i = 1:length(x)
    % 对于排放量数据
    text(x(i) + 0.18, y2(i) + 0.2, sprintf('%.1f', y2(i)), ...
         'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center', ...
         'FontName', 'Arial', 'FontSize', 12, 'Color', 'k', 'FontWeight','bold');
end

box on;
axis square;
lgd = legend('Thermal Power', 'Green Power', 'CO_2 Emission');
% 设置图例位置为右下角
lgd.Location = 'north';
% 使图例外部适应绘图区域
lgd.Box = 'off';
% 横向排列图例项
lgd.Orientation = 'horizontal';

% 修改字体为 Arial，大小为 16
set(gca, 'FontName', 'Arial', 'FontSize', 16);
set(lgd, 'FontName', 'Arial', 'FontSize', 12); % 确保图例也使用相同的字体

% 如果有标题或标签，你也需要对它们进行同样的设置：
% title('Your Title Here', 'FontName', 'Arial', 'FontSize', 16);
xlabel('Year', 'FontName', 'Arial', 'FontSize', 16);
yyaxis left;
ylabel('Internet Power (10^9 kWh)', 'FontName', 'Arial', 'FontSize', 16);
ylim([0, 35])
yyaxis right;
ylabel('Annual CO_2 Emission (Mt/yr)', 'FontName', 'Arial', 'FontSize', 16);
ylim([0, 35])

% 对于双 Y 轴的情况，你可能还需要分别设置每个 Y 轴的字体：
hAx.YAxis(1).FontName = 'Arial'; 
hAx.YAxis(1).FontSize = 16;
hAx.YAxis(2).FontName = 'Arial'; 
hAx.YAxis(2).FontSize = 16;

%% 导出图片
figname = ['C:\Users\huayi\Desktop\NC论文\Figure2\大唐托克托全年发电量\', 'Togtoh上网电量与年排放统计结果.png'];
print(fig2, figname, '-dpng', '-r600');