r = 50;
NumPoints = 500;
[UnitX,UnitY,UnitZ] = sphere (NumPoints);
RealX = r * UnitX;
RealY = r * UnitY;
RealZ = r * UnitZ;
x = RealX(:);
y = RealY(:);
z = RealZ(:);
Mask = 2.*y .* (sqrt(r.^2-x.^2-y.^2)./(r.^2)).* sind(15) ...
+ ((r.^2-2*x.^2-2*y.^2)./(r.^2)) .* cosd(15) ...
>= cosd(5.8);
% ----------------
figure ('Position', [62,182,1172,338]);
MaskPatchPositive = all ([Mask,z>0], 2);
MaskPatchNegative = all ([Mask,z<0], 2);
% ----------------
subplot (1,3, 1);
hold on
surf (RealX,RealY,RealZ, ...
'EdgeAlpha', 0.05, ...
'FaceAlpha', 0.05);
scatter3 (x(Mask),y(Mask),z(Mask),9, ...
'MarkerEdgeColor', 'none', ...
'MarkerFaceColor', 'r');
hold off
view ([-35,45]);
grid ('on');
box ('on');
daspect ('auto');
ax = gca;
ax.BoxStyle = 'full';
title ('All points');
% ----------------
subplot (1,3, 2);
hold on
surf (RealX,RealY,RealZ, ...
'EdgeAlpha', 0.1, ...
'FaceAlpha', 0.1);
scatter3 (x(MaskPatchPositive),y(MaskPatchPositive),z(MaskPatchPositive),9, ...
'MarkerEdgeColor', 'none', ...
'MarkerFaceColor', 'r');
hold off
xlim ([-5,5]);
ylim ([0,15]);
zlim ([45,50]);
grid ('on');
box ('on');
daspect ([1,1,1]);
view ([-60,30]);
ax = gca;
ax.BoxStyle = 'full';
title ('Upper part');
% ----------------
subplot (1,3, 3);
hold on
surf (RealX,RealY,RealZ, ...
'EdgeAlpha', 0.1, ...
'FaceAlpha', 0.1);
scatter3 (x(MaskPatchNegative),y(MaskPatchNegative),z(MaskPatchNegative),9, ...
'MarkerEdgeColor', 'none', ...
'MarkerFaceColor', 'r');
hold off
xlim ([-5,5]);
ylim ([0,15]);
zlim ([-50,-45]);
grid ('on');
box ('on');
daspect ([1,1,1]);
view ([-60,30]);
ax = gca;
ax.BoxStyle = 'full';
title ('Lower part');
