function [k,radiance_total] = Calculate_k(FWHM_SWIR, Gas, PARM_INPUT)
%% 参数：波数起止、步长
FWHM = FWHM_SWIR(Gas.bandindex, :); % 提取所选波段
PARM_INPUT.WAVNUM1 = round(1e7 / Gas.Band_wav2 - 80);
PARM_INPUT.WAVNUM2 = round(1e7 / Gas.Band_wav1 + 80);
STEP = 1; % MODTRAN 运行步长，设置为1
%% 单位转换（波长转波数）
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
    %% 输入文件
    PARM_INPUT.CO2MX = c_CO2(j);
    TAPE5_SET(PARM_INPUT);
    % write_tape5(MODEL, CO2MX, H2OVP, ANGLE, WAVNUM1, WAVNUM2, albedo);
    %% 提取结果
    radiance_total{j} = execute_mod(); % W/(cm2 sr cm-1)
end
%% 截取波段
radiance_cut = cell(NUM, 1);
for j = 1 : NUM % 初始化
    radiance_cut{j} = zeros(BandLen, 200);
end
for i = 1 : BandLen
    wavenumber = FWHM_(i, 1); % 取指定波数
    g = FWHM_(i, 2); % 半峰全宽 或者叫 resolution
    g = g / 2; % 半宽
    [~,index] = min(abs(radiance_total{1}(:, 1) - wavenumber)); % 找相应索引
    sta_index = round(index - floor(g) / STEP); % 起始索引
    stp_index = round(index + floor(g) / STEP); % 终止索引
    len = length(sta_index:stp_index);
    for j = 1 : NUM
        radiance_cut{j}(i, 1:len) = radiance_total{j}(sta_index:stp_index, 2)';
    end
end

%% 狭缝函数
slit = zeros(BandLen, 200); % 初始化
for i = 1 : BandLen
    g = FWHM_(i, 2); % 半峰全宽 或者叫 resolution
    g = g / 2; % 半宽
    window = -floor(g) : STEP : floor(g);
    len = length(window);
    % if(strcmp(satellite,'GF5') || strcmp(satellite,'GF5A') || strcmp(satellite,'GF5B') || ...
    %         strcmp(satellite,'EMIT') || strcmp(satellite,'PRISMA')  || strcmp(satellite,'ENMAP'))
    %     slit(i, 1:len) = sqrt(log(2)) / (sqrt(pi) * g) * exp(-log(2)*(window/g).^2); % 公式与HAPI.py 一致 高斯狭缝函数
    % elseif(strcmp(satellite,'ZY1E') || strcmp(satellite,'ZY1F'))
    %     a = 0.5;
    %     slit(i, 1:len) = 0.5 * sqrt(log(2)) / (sqrt(pi) * g / 2) * ((exp(-log(2)*((window-a*g/2)/(g/2)).^2) + exp(-log(2)*((window+a*g/2)/(g/2)).^2))); % 高斯叠加，用于ZY1
    % end
    slit(i, 1:len) = 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; % W/(cm2 sr cm-1)
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 = [斜率 截距]
    k(i) = -1 * p(1);
end

if(sum(k) == 0)
    disp('k值计算错误，请检查代码！')
    return
end
end
