%%
% v5: 仅包含所有AHSI结果
% v6：添加了EMIT的计算
% v7：添加了PRISMA的计算
% v8：修改了tape5参数，改成slant path，增加TAPE5_SET函数 等距排列HEIGHT 0.1:0.2:2.5
% v9: 非等距排列HEIGHT [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0];
% v10
% v25.1.1 PARM_INPUT.albedo 0.2
% v25.1.2 PARM_INPUT.albedo 0.5
%% 清理变量
% clear; close all;
%% 加载路径
fun_path = genpath('functions'); addpath(fun_path)
%%
version = '25.1.2';
filename_xlsx = ['CO2_k_v', version, '.xlsx'];
fpath_xlsx = fullfile('data', filename_xlsx);
filename_mat = ['CO2_k_v', version, '.mat'];
fpath_mat = fullfile('data\CO2_k', filename_mat);
PARM_INPUT.albedo = 0.5; % 地表反射率 albedo
PARM_INPUT.H2OVP = 0.5; % 水汽压力 gram/cm2
PARM_INPUT.ANGLE = 180; % 俯仰角 (180 - SAZ,satellite zenith)
PARM_INPUT.MODEL = 6; % 1976美国标准大气模型
PARM_INPUT.HEIGHT = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0];
PARM_INPUT.n_iter = length(PARM_INPUT.HEIGHT); % 迭代次数

%% ZY1E
tic;
mydisp('开始计算ZY1E卫星的气体吸收系数');
load('data/FWHM_SWIR/ZY1E.mat')
k_ZY1E = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% ZY1F
tic;
mydisp('开始计算ZY1F卫星的气体吸收系数');
load('data/FWHM_SWIR/ZY1F.mat')
k_ZY1F = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% GF5B
tic;
mydisp('开始计算GF5B卫星的气体吸收系数');
load('data/FWHM_SWIR/GF5B.mat')
k_GF5B = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% GF5A
tic;
mydisp('开始计算GF5A卫星的气体吸收系数');
load('data/FWHM_SWIR/GF5A.mat')
k_GF5A = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% GF5
tic;
mydisp('开始计算GF5卫星的气体吸收系数');
load('data/FWHM_SWIR/GF5.mat')
k_GF5 = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% EMIT
tic;
mydisp('开始计算EMIT卫星的气体吸收系数');
load('data/FWHM_SWIR/EMIT.mat')
k_EMIT = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% PRISMA
tic;
mydisp('开始计算PRISMA卫星的气体吸收系数');
load('data/FWHM_SWIR/PRISMA.mat')
k_PRISMA = Get_data_sheet2(FWHM_SWIR, PARM_INPUT);
toc;
%% 将数据写入 Excel 文件
mydisp('开始写入 Excel 文件');
sheet_names = {'k_ZY1E', 'k_ZY1F', 'k_GF5B', 'k_GF5A', 'k_GF5', 'k_EMIT', 'k_PRISMA'};
data_sheets = {k_ZY1E, k_ZY1F, k_GF5B, k_GF5A, k_GF5, k_EMIT, k_PRISMA};
% 列名
column_names = {'Wavelength', 'FWHM', 'k_Summer_L', 'k_Summer_M' , 'k_Summer_H' ,...
    'k_Winter_L','k_Winter_M','k_Winter_H', 'k_1976_L', 'k_1976_M', 'k_1976_H'};
for i = 1:length(sheet_names)
    writecell(column_names, fpath_xlsx, 'Sheet', sheet_names{i}, 'Range', 'A1'); % 写入列名
    writematrix(data_sheets{i}, fpath_xlsx, 'Sheet', sheet_names{i}, 'Range', 'A2'); % 写入数据
end
mydisp('数据写入完成');

%% 将数据写入.mat文件
save(fpath_mat, 'k_ZY1E', 'k_ZY1F', 'k_GF5B', 'k_GF5A', 'k_GF5', 'k_EMIT', 'k_PRISMA', 'PARM_INPUT');

%%
function data_sheet = Get_data_sheet2(FWHM_SWIR, PARM_INPUT)
Gas.Band_wav1 = FWHM_SWIR(1, 1);
Gas.Band_wav2 = FWHM_SWIR(end, 1);
Gas.bandindex = GetWavIndexes(FWHM_SWIR, Gas);

HEIGHT = PARM_INPUT.HEIGHT;
k_1976 = zeros(length(Gas.bandindex), PARM_INPUT.n_iter);

for i = 1:PARM_INPUT.n_iter
    PARM_INPUT.HEIGHT = HEIGHT(i); % km
    [k_1976(:, i), ~] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT);
end
mydisp('模型计算完成')

data_sheet = [FWHM_SWIR, k_1976];
end

%% 考虑水汽和季节影响的气体吸收系数
function data_sheet = Get_data_sheet(FWHM_SWIR, PARM_INPUT, satellite)
Gas.Band_wav1 = FWHM_SWIR(1, 1);
Gas.Band_wav2 = FWHM_SWIR(end, 1);
Gas.bandindex = GetWavIndexes(FWHM_SWIR, Gas);
H2OVP = [0.15, 0.5, 1.3];

% 1976大气模型
PARM_INPUT.MODEL = 6;
PARM_INPUT.H2OVP = H2OVP(1);
[k_1976_L, radiance_79L] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
PARM_INPUT.H2OVP = H2OVP(2);
[k_1976_M, radiance_76M] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
PARM_INPUT.H2OVP = H2OVP(3);
[k_1976_H, radiance_76H] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
k_1976 = [k_1976_L k_1976_M k_1976_H];
mydisp('1976大气模型计算完成')

% 中纬度夏季模型
PARM_INPUT.MODEL = 2;
PARM_INPUT.H2OVP = H2OVP(1);
[k_Summer_L, radiance_SL] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
PARM_INPUT.H2OVP = H2OVP(2);
[k_Summer_M, radiance_SM] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
PARM_INPUT.H2OVP = H2OVP(3);
[k_Summer_H, radiance_SH] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
k_Summer = [k_Summer_L k_Summer_M k_Summer_H];
mydisp('中纬度夏季模型计算完成')

% 中纬度冬季模型
PARM_INPUT.H2OVP = H2OVP(1);
[k_Winter_L, radiance_WL] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
PARM_INPUT.H2OVP = H2OVP(2);
[k_Winter_M, radiance_WM] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
PARM_INPUT.H2OVP = H2OVP(3);
[k_Winter_H, radiance_WH] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT, satellite);
k_Winter = [k_Winter_L k_Winter_M k_Winter_H];
mydisp('中纬度冬季模型计算完成')
data_sheet = [FWHM_SWIR, k_Summer, k_Winter, k_1976];
% save('radiance_76M.mat', 'radiance_76M');
end