r = sqrt(force_x.^2 + force_y.^2 + force_z.^2);
theta = atan2(force_y, force_x);
phi = acos(force_z ./ r);
figure('Color', 'white');
theta_edges = linspace(-pi, pi, n_theta+1);
phi_edges = linspace(0, pi, n_phi+1);
N = histcounts2(theta, phi, theta_edges, phi_edges);
theta_center = (theta_edges(i) + theta_edges(i+1))/2;
phi_center = (phi_edges(j) + phi_edges(j+1))/2;
dtheta = (theta_edges(i+1) - theta_edges(i));
dphi = (phi_edges(j+1) - phi_edges(j));
theta_patch = [theta_center-dtheta/2, theta_center+dtheta/2, ...
theta_center+dtheta/2, theta_center-dtheta/2];
phi_patch = [phi_center-dphi/2, phi_center-dphi/2, ...
phi_center+dphi/2, phi_center+dphi/2];
[x, y, z] = sph2cart(theta_patch, pi/2-phi_patch, scale);
patch(x, y, z, scale, 'EdgeColor', 'none');
title('3D Force Distribution Rose Plot');