function [ result ] = mu_rel_iron(r,A,Ar,Az)

    persistent data_B data_mu
    % Persistent cache of B-H data
    if( isempty(data_mu) )  % Load and define on first run
        % Load B-H curve data for iron (Tesla vs relative permeability)
        load BHcurve.mat;
        data_B = BHcurve(:,1);    % Magnetic flux density values [T]
        data_mu = BHcurve(:,2);  % Corresponding relative permeability values
    end
    
    % Calculate magnetic flux density magnitude
    if r < 0.01
        B_magnitude = ((2.*Ar).^2+(-Az).^2).^(1/2);
    else
        B_magnitude = ((A./r + Ar).^2+(-Az).^2).^(1/2);
    end

    array_check = isnan(B_magnitude) | isinf(B_magnitude);
    B_magnitude(array_check) = 0;
    B_magnitude = abs(B_magnitude);

    % Interpolate relative permeability
    result = interp1( data_B, data_mu, B_magnitude, 'linear', 'extrap' );
    
    % Apply bounds and smoothing if needed    
    % Handle saturation region
    index_limit = find( result < 1.01 );
    result(index_limit) = 1.01;  % Minimum value to avoid numerical issues
    % result = result.*0+100;  % Minimum value to avoid numerical issues
end
