# How to calculate arbitrary area on surface of sphere?

I am trying to calculate the solid angle subtended by arbitrary-shaped loops on a sphere's surface.

First, I parametrize circular loops by:

$$\theta(t,k_{x0},r) = k_{x_0} + r \cos(t);$$ $$\phi(t,k_{y0},r) = k_{y_0} + r \sin(t);$$

where $$0\leq t\leq2\pi$$, and $$k_{x_0}$$, $$k_{y_0}$$ define the loop's center. So, we can say that this step draws out a circular loop on the $$\theta/\phi$$ plane.

Then I project these onto the sphere's surface using spherical coordinates, as follows:

$$x(\theta,\phi)=r \cos{\theta}\sin{\phi}, y(\theta,\phi)=r \sin{\theta}\sin{\phi}, z(\theta,\phi)=r \cos{\phi}$$

How do I go about calculating the surface area within these $$(x,y,z)$$ loops on the surface? This will allow me to calculate the solid angle I need.

The solid angle is given by: $$\Omega= \iint_S \frac{\hat{r}\cdot\hat{n}}{r^2} \, \mathrm{d}\Sigma = \iint_{\mathcal{R}}\sin \theta \, \mathrm{d} \theta \, \mathrm{d} \phi=\frac{\textrm{spherical surface area}}{r^2}$$

I tried using various types of RegionMeasures to calculate this area (such as defining the area within the loop on the sphere as a Region, and by varying the radius from 0 to r, calculating the length of each loop in between and summing it all up), but I feel like I am missing a simple answer to my problem. Maybe what I am missing is a way to somehow map my arbitrary loops into a appropriate integration bounds for $$\Omega$$, but I tried to avoid this by resorting to Mathematica.

So far, I found the following post most useful: Integrate to calculate enclosed area

Note: I am parametrizing these loops in a peculiar way because I am trying to investigate a physics problem where the functions $$x(\theta,\phi),y(\theta,\phi),z(\theta,\phi)$$ will be different, and make my loops twist and turn. The ultimate goal is to find the solid angle in these cases, but I wanted to start with the sphere.

• @Bill thank you, I did see that. However, I was looking for an answer that did not rely entirely on r, theta, a or h the way they are defined, because I am trying to build up to the case where my base shape for computing the solid angle w.r.t. the loop is not the unit sphere but some arbitrary shape. Maybe I am just confusing myself unnecessarily. – TribalChief Feb 25 '19 at 20:18

v = {0, 0, 1};
r = Sin[Pi/4];
R = ImplicitRegion[
{x^2 + y^2 + z^2 == 1, v.{x, y, z} >= Sqrt[1 - r^2]},
{{x, -1.1, 1.1}, {y, -1.1, 1.1}, {z, -1.1, 1.1}}
];
Show[
DiscretizeRegion[R],
Graphics3D[Sphere[]]
]
Area[R] 1.8403

• Thank you for your answer. Just so I don’t confuse myself, would you mind clarifying your choices of numbers? Specifically, what are the roles of your chosen r and v? How do I choose the interval -1.1,1.1 based off my loop? – TribalChief Feb 25 '19 at 20:16
• Oh, v is the vector that point from {0,0,0} to the center of the disc and r is the radius. I chose r = Sin[Pi/4] just for having a 45-degree opening angle. The intervals -1.1 and 1.1 in the ImplicitRegion are just arbitrary. They have to be a bit larger than the sphere. If they are not extensively large, it helps DiscretizeRegion to localized the region to discretize (think of searching for the disk in a cube of edge length 1000 - you would need quite some luck to find it). – Henrik Schumacher Feb 25 '19 at 20:21
• Thank you for clarifying! However, how do I work with my definition for a loop? I am trying to apply this to a physics problem and so I have to start from a parameteization like the one I have in the question. Sorry for being slow, but I am missing the connection. – TribalChief Feb 25 '19 at 20:23
• Honesty, I do not understand your parameterization. It parameterizes a curve in the plane, not on the sphere. I also don't get your statement "Then I project these onto the sphere's surface using spherical coordinates." Maybe you want to share the code that you used? – Henrik Schumacher Feb 25 '19 at 20:26
• I did try to post it, but StackExchange kept giving me an error saying that it wasn’t formatted properly even though there was no apparent issue! I tried to type equations 4-6 on mathworld.wolfram.com/SphericalCoordinates.html. The kx and ky correspond to theta and phi. So, you are right in the sense that I choose a loop on a plane, but my axes here are the angles. Sorry that wasn’t clear. – TribalChief Feb 25 '19 at 20:30

Another approach, while slow, is to use RegionIntersection with geometry primitives. For example:

Area @ RegionIntersection[
Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1/Sqrt],
Sphere[{0,0,0}, 1]
]
% //N


-(-2 + Sqrt) 𝜋

1.8403

• Thank you for your answer. However, as with my comment to Henrik’s response, how do I directly relate my circular loops to the regions you’ve chosen? I am trying to look at a physics problem with this, and so my circular loops have to be parameteized that way. Sorry for being slow, but I’m missing the link. – TribalChief Feb 25 '19 at 20:21
• I just edited my question in an attempt to make it clearer. – TribalChief Feb 26 '19 at 1:09

A simple solution here is be to use the Geo functionality, in particular GeoArea, which can compute the area of any polygon on the surface of the sphere (or of an ellipsoid of revolution). There is also the primitive GeoCircle, consisting of the points on the spheroid at a given distance from a given center.

Note that in the Geo functionality we use degrees, not radians.

For example suppose your disk has a center with the following latitude and longitude, radius and is on a sphere of this radius:

centerLat = Quantity[30, "AngularDegrees"];
centerLon = Quantity[-40, "AngularDegrees"];


Then we can compute the area in stereoradians as

QuantityMagnitude[GeoArea[GeoDisk[{centerlat, centerlon}, geodiskRadius], GeoModel -> sphereRadius], "Meters"^2]
(* 8.89791 *)


You can also compute the area of any spherical polygon (with edges being great circles). For example this spherical triangle covers an octant of the full surface:

geoTriangle = Polygon[GeoPosition[{{0, 0}, {90, 0}, {0, 90}}]];


Therefore its area is 4Pi / 8:

QuantityMagnitude[GeoArea[geoTriangle, GeoModel -> sphereRadius], "Meters"^2]
(* 1.5708 *)


As I said, these functions can work on ellipsoids of revolution. Use GeoModel -> {a, b} where a and b are the equatorial and polar radiuses respectively.