clc;
clear;
close all;
fun_path = genpath('colorbar'); addpath(fun_path)
%%
foldername = 'E:\PRISMA\L0data';
%% 设置目标
facility = 'Matla';
files = Dirfiles(foldername, facility);
if(strcmp(facility, 'White Bluff'))
    Item.lat = 34.4234830;
    Item.lon = -92.1390973;
elseif(strcmp(facility, 'Labadie'))
    Item.lat = 38.5583;
    Item.lon = -90.8361;
elseif(strcmp(facility, 'Craig'))
    Item.lat = 40.4627;
    Item.lon = -107.5912;
elseif(strcmp(facility, 'Martin Lake'))
    Item.lat = 32.2597;
    Item.lon = -94.5703;
elseif(strcmp(facility, 'Matla'))
    Item.lat = -26.2815;
    Item.lon = 29.1415;
elseif(strcmp(facility, 'Matimba'))
    Item.lat = -23.7055;
    Item.lon = 27.5626;
elseif(strcmp(facility, 'PLU'))
    Item.lat = -7.71446;
    Item.lon = 113.58326;
end
%% 指定文件路径
len = length(files);
for j = 1:len
    close all;
    %% 读取数据
    filename = files(j).name;
    fname = fullfile(foldername, facility, filename);
    finfo = h5info(fname); % 获取文件信息
    disp(filename)
    %%
    swir = h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Data Fields/SWIR_Cube');
    vnir = h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Data Fields/VNIR_Cube');
    swir = single(permute(swir(:,1:171,:), [1, 3, 2]));
    swir = swir / 100; % 辐射定标
    swir = flip(swir, 3); % 翻转，短波长在前
    vnir = single(permute(vnir(:,4:66,:), [1, 3, 2]));
    vnir = vnir / 100; % 辐射定标
    vnir = flip(vnir, 3); % 翻转，短波长在前
    %% 经纬度坐标
    Latitude_SWIR = h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Geolocation Fields/Latitude_SWIR');
    Longitude_SWIR = h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Geolocation Fields/Longitude_SWIR');
    %% Data Fields
    % Cloud_Mask = single(h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Data Fields/Cloud_Mask'));
    % Cloud_Mask(Cloud_Mask == 255) = nan;
    %
    % LandCover_Mask = single(h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Data Fields/LandCover_Mask'));
    % LandCover_Mask(LandCover_Mask == 255) = nan;
    % SunGlint_Mask =  single(h5read(fname, '/HDFEOS/SWATHS/PRS_L1_HCO/Data Fields/SunGlint_Mask'));
    % SunGlint_Mask(SunGlint_Mask == 255) = nan;

    %%
    Cw_Swir = h5read(fname, '/KDP_AUX/Cw_Swir_Matrix');
    Fwhm_Swir = h5read(fname, '/KDP_AUX/Fwhm_Swir_Matrix');
    Cw_Vnir = h5read(fname, '/KDP_AUX/Cw_Vnir_Matrix');
    Fwhm_Vnir = h5read(fname, '/KDP_AUX/Fwhm_Vnir_Matrix');
    %%
    CW_VN = mean(Cw_Vnir(100:162,:), 2); % 4-66
    FWHM_VN = mean(Fwhm_Vnir(100:162,:), 2); % 4-66
    CW_SW = mean(Cw_Swir(82:252,:), 2); % 1-171
    FWHM_SW = mean(Fwhm_Swir(82:252,:), 2); % 1-171
    FWHM_VNIR = double([CW_VN, FWHM_VN]);
    FWHM_VNIR = flipud(FWHM_VNIR);
    FWHM_SWIR = double([CW_SW, FWHM_SW]);
    FWHM_SWIR = flipud(FWHM_SWIR);
    %% 加载气体吸收系数
    load("data/FWHM_SWIR/PRISMA.mat");
    load("data/CO2_k/CO2_k_v7.mat", "k_PRISMA");
    %% 读取反演波段
    datafile = 'data/Band_CO2_v1.xlsx';
    Tabledata = readtable(datafile, 'VariableNamingRule', 'preserve');
    bandindex = GetBandIndex(Tabledata, 'PRISMA', FWHM_SWIR);
    %% 匹配滤波算法
    bandlength = FWHM_SWIR(bandindex, 1);
    IsAC = false; % 不使用albedo factor矫正
    swir_ = swir(:, :, bandindex);
    swir_ = PreProcessSWIR(swir_);
    swir_2 = ConvertUnit(swir_, bandlength);
    %% 反演结果
    k_ = k_PRISMA(bandindex, 3);
    result = RMF(swir_2, k_, IsAC);
    %%
    facility_lat = Item.lat;
    facility_lon = Item.lon;
    % 计算每个点与给定经纬度的距离
    distances = sqrt((Latitude_SWIR - facility_lat).^2 + (Longitude_SWIR - facility_lon).^2);
    % 找到最小距离的索引
    [min_dist, linear_idx] = min(distances(:));
    % 转换为行列坐标
    [row, col] = ind2sub(size(distances), linear_idx);
    % 输出结果
    fprintf('最近的点坐标在 (%d, %d)，经纬度为 (%.4f, %.4f)\n', row, col, Latitude_SWIR(row, col), Longitude_SWIR(row, col));
    %%
    figure
    range = 70;
    imagesc(result(row-range:row+range,col-range:col+range));
    axis square;
    colorbar;
    toc;

end
%%
function files = Dirfiles(foldername, facility)
% 定义目录路径
directory = fullfile(foldername, facility);

% 获取该目录下的所有文件信息
allFiles = dir(fullfile(directory, '*.*'));

% 创建逻辑索引以排除 "." 和 ".."
validFilesIdx = ~ismember({allFiles.name}, {'.', '..'});

% 应用逻辑索引过滤文件列表
files = allFiles(validFilesIdx);
end