function [Fzs, Rstar] = Fz_soil_subcircle(h, R, bw, soil, Csig, n, c1, c2, s, f)

    % --- Eq (20): substitute circle ---
    if f <= 0
        Rstar = R;
    else
        h_safe = max(h, 1e-6);
        Rstar = R * ( sqrt(1 + f/h_safe) + sqrt(f/h_safe) )^2;
    end

    % --- Contact angles (con Rstar) ---
    theta_f = acos( max(-1, min(1, 1 - h/Rstar)) );
    theta_r = -theta_f;
    theta_m = (c1 + c2*s) * theta_f;

    % --- Bekker sigma_max (come fai tu) ---
    h_m = Rstar*(1 - cos(theta_m));
    sigma_max = Csig * h_m^n;

    % --- discretizzazione ---
    N = 300;
    th = linspace(theta_r, theta_f, N);
    dth = th(2) - th(1);

    % --- sigma(th): triangolare come in vertical_load ---
    sigma = zeros(size(th));
    left  = th <= theta_m;
    right = th >  theta_m;

    sigma(left)  = sigma_max .* (th(left)  - theta_r) ./ (theta_m - theta_r);
    sigma(right) = sigma_max .* (theta_f - th(right)) ./ (theta_f - theta_m);

    % --- tau_x(th): Iwanosi-Hanamoto (stesso spirito del tuo) ---
    % j(th) = R* [ (theta_f - th) - (1-s)(sin(theta_f)-sin(th)) ]
    % j = Rstar * ( (theta_f - th) - (1 - s)*(sin(theta_f) - sin(th)) );
    % 
    % tau = (soil.c + sigma*tan(soil.phi)) .* (1 - exp(-abs(j)/soil.K));
    % 
    % % --- Eq (18): Fzs = b R* ∫ [ tau sin(th) + sigma cos(th) ] dth
    % Fzs = bw * Rstar * sum( tau.*sin(th) + sigma.*cos(th) ) * dth;

    Fzs = bw * Rstar * sum(sigma) * dth;
end
