# Find equation of the circular cross section of a unit sphere

I have a unit sphere in Cartesian coordinates:

$x^2 + y^2 + z^2 = 1$

or in spherical coordinates:

$x = \rho \sin(\phi) \cos(\theta)\\ y = \rho \sin(\phi) \sin(\theta)\\ z = \rho \cos(\phi)$

I select a point $P$ on the surface of the sphere, in the coordinate system that is more convenient. I want to identify the points, on the surface of the sphere, on the circle of radius $r$ centred in $P$, where $r$ is the distance of the shortest arc on the surface of the sphere from each point to $P$ (i.e. $r$ is the great arc between two points on a unit sphere).

I expect this problem to be pretty simple but I haven't found any clear resource on the web. Any idea? Essentially what I need is the equation in 3D of a "circular section" of a sphere, possibly as a function of one angle.

• I think you should define better the circular section centered in $P$. – Tony Piccolo Nov 4 '13 at 11:36
• As I understand it, he means the points that are at distance $r$ from $P$ along a great arc/geodesic through $P$. – HSN Nov 4 '13 at 11:38
• Do you mean a circle section of the sphere with a plane distant $\cos r\,$ ($r$ radians) from the center of the sphere ? – Tony Piccolo Nov 4 '13 at 11:53
• I meant the points on the surface of the sphere at shortest distance $r$ from $P$. The distance $r$ is measured along the surface of the sphere, i.e. it measures the great arc from each point to P. – randomatlabuser Nov 4 '13 at 11:58
• So I think I am right: please check. – Tony Piccolo Nov 4 '13 at 12:18

I'm not entirely sure if you're asking for a closed expression that gives the points, or just a description of these points and am going to assume the latter.

The first step in solving this problem is to change it into an easier one. We may rotate the sphere in such a way the point $P$ becomes the north pole $N=(0,0,1)$. If we call this rotation (which may be expressed as a $3\times3$ matrix) $R$, a point $Q$ is on the circle with radius $r$ around $P$ iff $R(Q)$ is on the circle with radius $r$ around $N$. This is due to the fact that a rotation of the sphere does not change distances on the sphere.

The next step is to consider these particular circles of radius $r$ around $N$. The best way to see this is to draw a sketch of the intersection of the sphere with the $yz$-plane, i.e. a circle of radius 1 around the origin. We can also see to points $P_1, P_2$ on the circle at distance $r$ from $N$. (I'll assume $r$ is small enough for this to make sense) These points should be such that $P_1$ is obtained from $P_2$ by reflecting around the $z$-axis. The spherical coordinate $\theta$ corresponding to these points can now easily be obtained: this angle is given as an angle wrt the $z$-axis, while $r$ is precisely the arc length corresponding to this angle on a circle with radius 1. This leads to the conclusion that $r=\theta$ in this special case.

The next thing to notice, is that all points at distance $r$ from $N$ have the same $z$-coordinate. By projecting either $P_1$ or $P_2$ at this $z$-coordinate (which again can be done insightful in the same picture), one finds this height to be $z_0=\cos(\theta)=\cos(r)$. Moreover, yet again from considering the same picture, the required circle is a circle parallel to the $xy$-plane with center $(0,0,z_0)$ and radius $\sin(\theta)=\sin(r)$.

We may thus conclude that the points at distance $r$ from $N$ are of the form $(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))$ with $\phi\in[0,2\pi)$. Using the rotation matrix $R$ from before, we thus get the set of points satisfying your original question: $$\{R^{-1}(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))\mid \phi\in[0,2\pi)\}.$$

• Yes, that's it! – randomatlabuser Nov 5 '13 at 1:39

This is what I have done so far, in Matlab:

% This is what I have done so far, in Matlab
% Spherical coordinates of a point on the surface of the unit sphere
theta_sph = 2 * rand * pi;    % first angle
phi_sph = 2 * rand * pi;    % second angle
% Cartesian coordinates of the same point
x_sph = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph = rho_sph * cos(phi_sph);
xyz_sph = [x_sph; y_sph; z_sph];

% Find Rotation matrix to rotate the sphere
% so that the point will end up at the North pole
a_x = 0; a_y = -phi_sph; a_z = -theta_sph;
R_x = [1 0 0;
0 cos(a_x) -sin(a_x);
0 sin(a_x) cos(a_x)];
R_y = [cos(a_y) 0 sin(a_y);
0 1 0;
-sin(a_y) 0 cos(a_y)];
R_z = [cos(a_z) -sin(a_z) 0;
sin(a_z) cos(a_z) 0;
0 0 1];
% Rotation matrix (to be applied to the cartesian coordinates)
R = R_x * R_y * R_z;
% Rotation matrix to move back the point to its original coordinates
R_back = R^(-1);
% Coordinates on rotated sphere
xyz_rot = R * xyz_sph;
% Coordinates rotated back
xyz_rot_back = R_back * xyz_rot;
disp('Error is infinitesimal:')
disp(xyz_sph - xyz_rot_back)
fprintf('\n\n\n')

% Plot the circles on the sphere in the new rotated coordinates
figure
sphere
axis equal
hold on
plot3(xyz_sph(1), xyz_sph(2), xyz_sph(3), 'ow', 'MarkerFaceColor', 'w', 'MarkerSize', 7)
angles_circ = (0:360)' / 180 * pi;
m = 360;
dst_pnts1 = zeros(m, 1);
dst_pnts2 = zeros(m, 1);
chck = zeros(m, 3);
for hh = 1:m
% Distance from point, angles and coordinates of the circle
r = hh / m * pi;
x_new = rho_sph * sin(r) * cos(angles_circ);
y_new = rho_sph * sin(r) * sin(angles_circ);
z_new = rho_sph * cos(r) * ones(numel(angles_circ), 1);
% Finally, rotate back the points
xyz_new_rot = [x_new y_new z_new] * R_back';
plot3(xyz_new_rot(:, 1), xyz_new_rot(:, 2), xyz_new_rot(:, 3), '-k', 'LineWidth', 4)
% pause(0.05)
drawnow
% Calculate angles (    angle = atan2(norm(cross(a, b)), dot(a, b))    )
a = xyz_new_rot';
b = repmat(xyz_rot_back, 1, numel(angles_circ));
tmp = cross(a, b);
d1 = rho_sph * atan2(arrayfun(@(n) norm(tmp(:, n)), 1:size(tmp,2)), dot(a, b));
if all(fix(1e12 * (d1 - d1(1))) == 0)
dst_pnts1(hh) = d1(1);
else
dst_pnts1(hh) = NaN;
end

% Finally, to retrieve the spherical coordinates:
theta_sph_new = atan2(xyz_new_rot(:, 2), xyz_new_rot(:, 1));
phi_sph_new = acos(xyz_new_rot(:, 3) / rho_sph);
% Check that both cartesian and spherical coordinates are correct:
tmp = rho_sph * sin(phi_sph_new) .* cos(theta_sph_new) - xyz_new_rot(:, 1);
chck(m, 1) = max(abs(tmp));
tmp = rho_sph * sin(phi_sph_new) .* sin(theta_sph_new) - xyz_new_rot(:, 2);
chck(m, 2) = max(abs(tmp));
tmp = rho_sph * cos(phi_sph_new) - xyz_new_rot(:, 3);
chck(m, 3) = max(abs(tmp));

% Evaluate the distance between points on the circle and their centre:
d2 = rho_sph * acos(cos(phi_sph_new) .* cos(repmat(phi_sph, numel(angles_circ), 1)) + ...
sin(phi_sph_new) .* sin(repmat(phi_sph, numel(angles_circ), 1)) .* ...
cos(theta_sph_new - repmat(theta_sph, numel(angles_circ), 1)));
if all(fix(1e12 * (d2 - d2(1))) == 0)
dst_pnts2(hh) = d2(1);
else
dst_pnts2(hh) = NaN;
end

end

% Check number of bad approximations in the conversion into spherical coordinates
disp('Check number of bad approximations in the conversion into spherical coordinates:')
sum(chck(:) > 1e-12)

% For each circle the distance is constant
disp('For each circle the distance is constant:')
~any(isnan(dst_pnts1))
~any(isnan(dst_pnts2))
fprintf('\n\n\n')

% Plot calculated distances vs to known distances
figure
plot((1:m) / m * pi, (dst_pnts1), '.')
set(gca, 'XLim', [0 pi])
set(gca, 'YLim', [0 pi])
set(gca, 'XTick', (0:8) * pi / 8)
set(gca, 'YTick', (0:8) * pi / 8)
set(gca, 'XTickLabel', {'0' 'pi/8' 'pi/4' '3pi/8' 'pi/2' '5pi/8' '3pi/4' '7pi/8' 'pi'})
set(gca, 'YTickLabel', {'0' 'pi/8' 'pi/4' '3pi/8' 'pi/2' '5pi/8' '3pi/4' '7pi/8' 'pi'})

% Plot calculated distances vs to known distances
figure
plot((1:m) / m * pi, (dst_pnts2), '.')
set(gca, 'XLim', [0 pi])
set(gca, 'YLim', [0 pi])
set(gca, 'XTick', (0:8) * pi / 8)
set(gca, 'YTick', (0:8) * pi / 8)
set(gca, 'XTickLabel', {'0' 'pi/8' 'pi/4' '3pi/8' 'pi/2' '5pi/8' '3pi/4' '7pi/8' 'pi'})
set(gca, 'YTickLabel', {'0' 'pi/8' 'pi/4' '3pi/8' 'pi/2' '5pi/8' '3pi/4' '7pi/8' 'pi'})

• Notice that $\rho$ in the spherical coordinates is just 1 in your case, since you take a sphere of radius 1. – HSN Nov 5 '13 at 14:27
• Yes, but the if you decide to change that value the code should still work because it takes $\rho$ as a parameter. – randomatlabuser Nov 5 '13 at 14:31