function W = vertical_load(h_f, rw, bw, Csig, n, c1, c2, s)

    theta_f = acos(1 - h_f/rw); % [rad]
    theta_r = -theta_f; % [rad]
    theta_m = (c1 + c2*s)*theta_f; % [rad]
    % theta_m = 0.5*theta_f; 

    h_m = rw*(1 - cos(theta_m));
    sigma_max = Csig * h_m^n;

    N = 200;
    theta = linspace(theta_r, theta_f, N);
    dtheta = theta(2) - theta(1);

    sigma = zeros(size(theta));

    mask_left  = theta <= theta_m;
    mask_right = theta >  theta_m;

    sigma(mask_left)  = sigma_max * ...
        (theta(mask_left) - theta_r) / (theta_m - theta_r);

    sigma(mask_right) = sigma_max * ...
        (theta_f - theta(mask_right)) / (theta_f - theta_m);

    W = bw * rw * sum(sigma) * dtheta;

end