function k = calc_k(FWHM_SWIR, Gas, xmlStruct)
%% 参数：波数起止、步长
FWHM = FWHM_SWIR(Gas.bandindex, :); % 提取所选波段
WAVNUM1 = round(1e7 / Gas.Band_wav2 - 220);
WAVNUM2 = round(1e7 / Gas.Band_wav1 + 220);
SatelliteZenith = str2double(xmlStruct.SatelliteZenith);

if(SatelliteZenith >= 0 && SatelliteZenith < 90)
    ANGLE = 180 - SatelliteZenith;
elseif(SatelliteZenith >= -90 && SatelliteZenith < 0)
    ANGLE = 180 - (90 + SatelliteZenith);
else
    ANGLE = 180;
    fprintf("请检查卫星俯仰角参数！\n");
end
ANGLE = 179.64;
MON = month(str2date(xmlStruct.StartTime));
albedo = 0.05;
STEP = 1; % MODTRAN 运行步长，设置为1
%% 参数文件
if(MON >= 6 && MON <= 8)
    MODEL = 2; % 中纬度夏季模型
elseif((MON == 1) || (MON >= 11 && MON <= 12))
    MODEL = 3; % 中纬度冬季模型
else
    MODEL = 6; % 默认1976大气参数
end
%% 单位转换（波长转波数）
FWHM_(:,1) = 1e7 ./ FWHM(:,1);
FWHM_(:,2) = 1e7 ./ (FWHM(:,1) - FWHM(:,2)/2) - 1e7 ./ (FWHM(:,1)+FWHM(:,2)/2);
BandLen = size(FWHM_, 1); % 波段个数
%%
NUM = 10; % NUM 表示执行MODTRAN的次数，也表示不同浓度的个数
radiance_total = cell(NUM, 1);
c_CO2 = linspace(400, 500, NUM); % 单位：ppmv
delta_c = (c_CO2(2:end) - c_CO2(1))'; % 浓度差
for j = 1 : NUM
    %% 输入文件
    CO2MX = c_CO2(j);
    write_tape5(MODEL, CO2MX, ANGLE, WAVNUM1, WAVNUM2, albedo);
    %% 提取结果
    radiance_total{j} = execute_mod();
end

%% 获取卷积半径
s = zeros(BandLen, 1);
for i = 1:BandLen
    wavenumber = FWHM_(i, 1); % 取指定波数
    [~,index] = min(abs(radiance_total{1}(:, 1) - wavenumber)); % 找相应索引
    s(i)=index;
end
s1 = s(1:end-1);
s2 = s(2:end);
s3 = abs(s2 - s1); % 作差
AF_Wing = 20;
%% 截取波段
radiance_cut = cell(NUM, 1);
for i = 1 : BandLen
    wavenumber = FWHM_(i, 1); % 取指定波数
    [~,index] = min(abs(radiance_total{1}(:, 1) - wavenumber)); % 找相应索引
    sta_index = round(index - AF_Wing / STEP); % 起始索引
    stp_index = round(index + AF_Wing / STEP); % 终止索引
    for j = 1 : NUM
        radiance_cut{j}(i,:) = radiance_total{j}(sta_index:stp_index, 2);
    end
end

%% 狭缝函数
window = -AF_Wing : STEP : AF_Wing;
slit = zeros(BandLen, length(window)); % 初始化
for i = 1 : BandLen
    g = FWHM_(i, 2); % g为半宽 或者叫 resolution
    g = g / 2;
    slit(i,:) = sqrt(log(2)) / (sqrt(pi) * g) * exp(-log(2)*(window/g).^2); % 公式与HAPI.py 一致
    slit(i,:) = slit(i,:) / (sum(slit(i,:)) * STEP);
end

%% 高斯卷积
radiance_conv = cell(NUM, 1);
for j = 1 : NUM
    radiance_conv{j} = sum(radiance_cut{j} .* slit, 2) * STEP;
end

%% 计算k值
y = zeros(BandLen, NUM - 1);
for j = 2 : NUM
    y(:, j - 1) = log(radiance_conv{j}./radiance_conv{1}); % 分别除以conv_radiance1
end
x = delta_c;
k = zeros(BandLen, 1); % 单位每ppmv

for i = 1 : BandLen
    % 最小二乘拟合结果 斜率即为 k
    p = polyfit(x, y(i, :), 1); % p = [斜率 截距]
    % p = mean(y(i, :));
    k(i) = -1 * p(1);
    R = corrcoef(x,y(i, :));
    R2 = R(1,2)^2;

    %     f = polyval(p,x);
    %     plot(x,y(i,:),'o',x,f,'-')
    %     legend('data','linear fit')
end

if(sum(k) == 0)
    disp('k值计算错误，请检查代码！')
    return
end
% plot(FWHM(:,1),k)
end
