""" Sphere, a geometry
"""
import numpy as _np
[docs]class Sphere:
def __init__(self, radius):
self._r = radius
[docs] def cut(self, normal: _np.ndarray, distance: float, resolution: int=361) -> list[_np.ndarray]:
""" Get a list of vectors produced by a cut by a plane
Get a set of vectors produced by a cut by a plane, defined by the normal vector and distance from the center.
This method is useful to get vectors that forms
- the equator (distance zero with respective to the pole as normal)
- a latitutde (distance corresponding to the latitude, ie. r sin(lat), with respctive to the pole)
- the terminator (distance 0, w.r.t the subsolar point)
- a solar zenith angle (distance corresponding to the sza, i.e., r cos(sza), w.r.t the subsolar point)
- a longitude (distance 0 w.r.t. a point on the equator 90 degrees off from the latitude)
:param normal: The normal vector of the cut plane. (Not necessarily normalized, length does not care).
:param distance: The distance from the center. Negative allowed.
:keyword resolution: The resolution.
>>> import numpy as np
>>> sphere = Sphere(1) # A unit sphere
>>> parallel = sphere.cut([0, 1, 0], 0.5) # Cut by a xz-plane at y=0.5.
>>> print(len(parallel)) # A list of (3,) ndarray
361
>>> parallel = np.array(parallel)
>>> y = parallel[:, 1] # y shall be constant
>>> print(f'{y.min():.2f} {y.max():.2f}')
0.50 0.50
>>> x = parallel[:, 2] # x and z should range sq(0.75)
>>> print(f'{x.min():.3f} {x.max():.3f}')
-0.866 0.866
>>> z = parallel[:, 2] # x and z
>>> print(f'{z.min():.3f} {z.max():.3f}')
-0.866 0.866
"""
r = self._r
d = distance
nvec = _np.array(normal)
nvec = nvec / _np.linalg.norm(nvec)
theta = _np.arccos(_np.clip(nvec[2], -1, 1))
phi = _np.arctan2(nvec[1], nvec[0])
ct = _np.cos(theta)
st = _np.sin(theta)
cp = _np.cos(phi)
sp = _np.sin(phi)
if d > r:
return [] # No intersection
elif d == r:
return [nvec * d] # A single point
# Vector set to define
t = _np.linspace(0, 2 * _np.pi, resolution)
xdd = _np.sqrt(r * r - d * d) * _np.cos(t)
ydd = _np.sqrt(r * r - d * d) * _np.sin(t)
zdd = _np.broadcast_to(d, xdd.shape)
# Rotating theta around y.
xd = ct * xdd + st * zdd
yd = ydd
zd = -st * xdd + ct * zdd
# Rotating phi around z.
x = cp * xd - sp * yd
y = sp * xd + cp * yd
z = zd
xyz = _np.array([x, y, z]) # (3, resol) shape
xyz = xyz.T
xyz = xyz.tolist()
return xyz
[docs] def equator(self, northpole, resolution=361):
"""
:param northpole: The north pole direction. (Length doesn't care)
>>> import numpy as np
>>> sphere = Sphere(1000)
>>> vec = sphere.equator([0, 1, 0]) # If y-axis is the north pole, the the equator is in x-z plane (y=0)
>>> vec = np.array(vec)
>>> vec.shape
(361, 3)
>>> print(f'{vec[:, 1].min():.2f} {vec[:, 1].max():.2f}')
-0.00 0.00
>>> print(f'{vec[:, 0].min():.2f} {vec[:, 0].max():.2f}')
-1000.00 1000.00
"""
return self.cut(northpole, 0, resolution=resolution)
[docs] def latitude(self, northpole, latitude_deg, resolution=361):
""" Return the set of vectors, represnting the latitude.
:param latitude_deg: Latitude in degrees.
:param northpole: North pole coordinates. (Length doesn't care)
>>> import numpy as np
>>> sphere = Sphere(1000)
>>> vec = np.array(sphere.latitude([0, 1, 0], 30)) # If y-axis is the north pole, the latitude 30 is y=500 km
>>> print(f'{vec[:, 1].min():.2f} {vec[:, 1].max():.2f}')
500.00 500.00
"""
d = self._r * _np.sin(_np.radians(latitude_deg))
return self.cut(northpole, d, resolution=resolution)
[docs] def terminator(self, sun_direction, resolution=361):
"""
:param sun_direction: Sun direction (from the center), or in other word, subsolar point coordinate.
(Length doesn't care)
>>> import numpy as np
>>> sphere = Sphere(1000)
>>> vec = np.array(sphere.terminator([-0.1, 0, 0])) # If -x axis is the sundirection, the terminator is x=0.
>>> print(f'{vec[:, 0].min():.2f} {vec[:, 0].max():.2f}')
-0.00 0.00
"""
return self.cut(sun_direction, 0, resolution=resolution)
[docs] def sza(self, sun_direction, sza_deg, resolution=361):
"""
:param sza_deg: the solar zenith angle.
:param sun_direction: The sun direction. (Length doesn't care)
>>> import numpy as np
>>> sphere = Sphere(1000)
>>> vec = np.array(sphere.sza([0, 0, -2], 60)) # If -z axis is the sun direction, the sza 60 is z=-500 km
>>> print(f'{vec[:, 2].min():.2f} {vec[:, 2].max():.2f}')
-500.00 -500.00
"""
d = self._r * _np.cos(_np.radians(sza_deg))
return self.cut(sun_direction, d, resolution=resolution)