function output_grid = k_(sand_path, silt_path, clay_path, soc_path)
    % 读取Geotiff数据及其地理信息
    [sand_data, R_sand] = geotiffread(char(sand_path));
    [silt_data, R_silt] = geotiffread(char(silt_path));
    [clay_data, R_clay] = geotiffread(char(clay_path));
    [soc_data, R_soc] = geotiffread(char(soc_path));

    % 将数据归一化到范围LO，
    sand_data = double(sand_data) / max(sand_data(:));
    silt_data = double(silt_data) / max(silt_data(:));
    clay_data = double(clay_data) / max(clay_data(:));
    soc_data = double(soc_data) / max(soc_data(:));

    % 检查是否所有的纹理信息都相同
    assert(isequal(R_sand, R_silt) && isequal(R_sand, R_clay) && isequal(R_sand, R_soc), 'Geographic information mismatch');

    % 计算值
    [rows, cols] = size(sand_data);
    k_values = zeros(rows, cols);
    for i = 1:rows
        for j = 1:cols
            k_values(i, j) = k_calcu([sand_data(i, j), silt_data(i, j), clay_data(i, j), soc_data(i, j)]);
        end
    end

    % 将k值归一化到范围内 [0, 1]
    normalized_k_values = (k_values - min(k_values(:))) / (max(k_values(:)) - min(k_values(:)));

    % 创建归一化k值的Geotiff
    output_filename = 'normalized_k_values.tif';
    geotiffwrite(output_filename, normalized_k_values, R_sand);

    output_grid = normalized_k_values;
end


function [k] = k_calcu(mat)
 disp(mat)
%K_CALCU 计算k值
%   m   砂粒%
%   f   粉沙粒%
%   n   黏粒%
%   t   有机碳质量分数%
if isnan(max(mat))
    k=nan;
else
    m=mat(1);% g/kg
    f=mat(2);
    n=mat(3);
    t=mat(4)/10;%soc单位为 dg/kg
    sub=m+f+n;
    m=m/sub*100;
    f=f/sub*100;
    n=n/sub*100;

    theta=1-(m/100);
    q1=0.3*exp(0.0256*m*(1-(f/100)))+0.2;
    q2=(f/(f+n))^(0.3);
    q3=-0.25*t*(exp(-2.96*t+3.72)+t)+1;
    q4=-0.7*theta/(exp(22.9*theta-5.51)+theta)+1;

    k=q1*q2*q3*q4;
end

end

function output_grid = k_(sand_path, silt_path, clay_path, soc_path,output_path)
    % 读取Geotiff数据及其地理信息
    [sand_data, R_sand] = geotiffread(char(sand_path));
    [silt_data, R_silt] = geotiffread(char(silt_path));
    [clay_data, R_clay] = geotiffread(char(clay_path));
    [soc_data, R_soc] = geotiffread(char(soc_path));

    % 将数据归一化到范围LO，
    sand_data = double(sand_data) / max(sand_data(:));
    silt_data = double(silt_data) / max(silt_data(:));
    clay_data = double(clay_data) / max(clay_data(:));
    soc_data = double(soc_data) / max(soc_data(:));

    % 检查是否所有的纹理信息都相同
    assert(isequal(R_sand, R_silt) && isequal(R_sand, R_clay) && isequal(R_sand, R_soc), 'Geographic information mismatch');

    % 计算值
    [rows, cols] = size(sand_data);
    k_values = zeros(rows, cols);
    for i = 1:rows
        for j = 1:cols
            k_values(i, j) = k_calcu([sand_data(i, j), silt_data(i, j), clay_data(i, j), soc_data(i, j)]);
        end
    end

    % 将k值归一化到范围内 [0, 1]
    normalized_k_values = (k_values - min(k_values(:))) / (max(k_values(:)) - min(k_values(:)));

    % 创建归一化k值的Geotiff
%     output_filename = 'normalized_k_values.tif';
    geotiffwrite(char(output_path), normalized_k_values, R_sand);

    output_grid = normalized_k_values;
end


function [k] = k_calcu(mat)
 disp(mat)
%K_CALCU 计算k值
%   m   砂粒%
%   f   粉沙粒%
%   n   黏粒%
%   t   有机碳质量分数%
if isnan(max(mat))
    k=nan;
else
    m=mat(1);% g/kg
    f=mat(2);
    n=mat(3);
    t=mat(4)/10;%soc单位为 dg/kg
    sub=m+f+n;
    m=m/sub*100;
    f=f/sub*100;
    n=n/sub*100;

    theta=1-(m/100);
    q1=0.3*exp(0.0256*m*(1-(f/100)))+0.2;
    q2=(f/(f+n))^(0.3);
    q3=-0.25*t*(exp(-2.96*t+3.72)+t)+1;
    q4=-0.7*theta/(exp(22.9*theta-5.51)+theta)+1;

    k=q1*q2*q3*q4;
end

end