function I_values1 = calculate_I_values( diss_no, Rxy, Zxy, pol, varargin)
% 计算I_values1的函数

    [~, ~, ~, num_time_steps] = size(diss_no);
    
    % 解析可选参数
    p = inputParser;
    addParameter(p, 'method', 'diss', @ischar);
    addParameter(p, 'R_center', 1.86, @isnumeric);
    parse(p, varargin{:});
    
    
    R_center = p.Results.R_center;
    method = p.Results.method;
    % 初始化输出
    I_values1 = zeros(num_time_steps, 1);
    
    % 预计算r和theta的梯度（如果网格不随时间变化）
    r = sqrt((Rxy - R_center).^2 + Zxy.^2);
    theta = pol;
    
    [dr_dx, dr_dy] = gradient(r);
    [dtheta_dx, dtheta_dy] = gradient(theta);
    
    % 计算雅可比行列式
    J = dr_dx .* dtheta_dy - dr_dy .* dtheta_dx;
    J(abs(J) < 1e-10) = 1e-10;
    
    % 循环计算每个时间步
    for t = 1:num_time_steps
        % 计算累计量
        if strcmp(method, 'diss')
            
            delta_Ni = diss_no(:, :, 1, t) ;
            % 计算当前时刻总量
        else if strcmp(method, 'exp')
            delta_Ni = (diss_no(:, :, 1, t));    
            % 计算当前时刻相对初始变化量
        else
            delta_Ni = (diss_no(:, :, 1, t) - diss_no(:, :, 1, 1)) ;
        end
        
        % 计算被积函数f
        f = r .* Rxy * 2 * pi .* delta_Ni;
        
        % 计算积分
        I_values1(t) = sum(sum(f .* abs(J)));
    end
end

使用方法：calculate_I_values( 被积分变量, Rxy, Zxy, pol, 可选参数)

被积分量应当为(x,y,z,t)四维矩阵
Rxy, Zxy, pol 从netcdf文件中读取
可选参数包括积分方式和磁轴位置，默认积分方式为计算所有时刻累积量