''' Calculate intersection between LoS and sphere.
.. codeauthor:: Yoshifumi Futaana
*Line of sight and sphere*
Given a observation point, *P*, and the looking vector, *V*, whether it crosses
a sphere with radius of R at the point C.
Use :func:`los_sphere` function.
To get limb position, use :func:`los_limb_sphere`.
'''
import math
import logging
_logger = logging.getLogger(__name__)
import numpy as _np
import numpy as np
from irfpy.util.vector3d import Vector3d
from irfpy.util.cone import Cone as _Cone
[docs]def los_sphere(pos, view, cent, radius, insphere='exception'):
r''' Returns the intersection of the LOS at the sphere surface.
::
* pos
view \
\ /---\
intersection * \
/ \
| * |
\ cent /
\ /
\---/ sphere with radius
:param pos: Position of observer
:type pos: ``numpy.array``
:param view: Viewing direction
:type view: ``numpy.array``
:param cent: Center of the target sphere
:type cent: ``numpy.array``
:param radius: Radius of the sphere
:type radius: ``float``
:keyword insphere: If the observer position is in sphere, how to handle?
"exception" will raise RuntimeError, "none" will return None.
:return: Position of the intersection as ``vector3d.Vector3d`` instance.
``None`` will be returned if no intersection is found.
:rtype: ``vector3d.Vector3d`` instance or ``None``.
>>> print(los_sphere([10,0,0], [-1,0,0], [-5,0,0], 3))
Vector3d( -2, 0, 0 )
>>> print(los_sphere([10,0,0], [0,1,0], [-5,0,0], 3))
None
>>> print(los_sphere([10, 0, 0], [1, 0, 0], [-5, 0, 0], 3))
None
If the given pos is in the sphere, RuntimeError will be raised as default.
>>> try:
... los_sphere([10,0,0], [-1,0,0], [-5,0,0], 30)
... print('RuntimeError should have been raised.')
... except RuntimeError as e:
... print('OK')
... except:
... print('Not correct exception.')
OK
'''
_logger = logging.getLogger(__name__ + '.los_sphere')
vpos = Vector3d(pos[0], pos[1], pos[2])
vview = Vector3d(view[0], view[1], view[2])
vcent = Vector3d(cent[0], cent[1], cent[2])
if vview.lengthSquared() == 0:
raise RuntimeError('View vector length is 0')
vview.normalize()
# Conversion of the coordinates. Centerize vpos.
vcent.sub(vpos)
# logger.debug('Vpos = %s (centerized)' % str(vpos))
_logger.debug('Vview = %s' % str(vview))
_logger.debug('Vcent = %s (Vpos centered)' % str(vcent))
# Distance**2 between P and C.
L_P_C2 = vcent.lengthSquared()
R2 = radius * radius
if (L_P_C2 <= R2):
if insphere.lower() in ('none',):
return None
else: # elif insphere.lower() in ('exception',):
raise RuntimeError('Given position P is in the sphere.')
# Limb cone angle is asin (R / L_P_C)
L_P_C = math.sqrt(L_P_C2)
limbang = radius / L_P_C
if limbang > 1:
# They may not happen but sometimes occurred by numerical error.
_logger.warning('limbangle out of domain (%f). Set to 1.' % limbang)
limbang = 1
elif limbang < -1:
_logger.warning('limbangle out of domain (%f). Set to -1.' % limbang)
limbang = -1
limbang = math.asin(limbang)
_logger.debug('Limb angle = %f' % (limbang * 180. / math.pi))
# Angle between view and P-C
viewangle = vview.angle(vcent)
_logger.debug('View angle = %f' % (viewangle * 180. / math.pi))
if viewangle > limbang:
return None
# The length to the sphere is, L cos t - sqrt( L^2 cos^2 t - (L^2 - R^2) )
cosviewangle = math.cos(viewangle)
det = L_P_C2 * cosviewangle ** 2 - (L_P_C2 - R2)
if det < 0:
_logger.warning('Determinant negative (%g). Set to 0.' % det)
det = 0
_logger.debug('det = %f' % det)
k = L_P_C * cosviewangle - math.sqrt(det) # Nearer.
vview.scale(k)
vview.add(vpos)
return vview
[docs]def intersections_sphere(pos, view, cent, radius):
r""" Returnt the intersection coordinates.
::
* pos
view \
\ /----\
intersection1* \
/ \ cent \
| \ * |
\ \ /
\ \ /
\--*-/ sphere with radius
intersections2
:param pos: (3,) or (N, 3) shaped array of positions
:param view: (3,) or (N, 3) shaped array of viweing directions (or velocity vector)
:param cent: The center of the reference sphere
:param radius: The radius of the reference spehre
:return: A tuple of (is1, is2, k). is1 and is2 is array with shape of (3,) or (N, 3). k is array with shape of (2,) or (N,2)
Samples. Let's consider the sphere cetnered at ``(1, 1, 0)`` with a radius of 1.
>>> cent = [1, 1, 0]
>>> r = 1
If one places at ``(-1, 1, 0)``, looking toward ``(1, 0, 0)``, then one sees the interesections
at ``(0, 1, 0)`` and ``(2, 1, 0)``.
>>> pos = [-1, 1, 0]
>>> view = [1, 0, 0]
>>> is1, is2, k = intersections_sphere(pos, view, cent, r)
>>> print(is1)
[0. 1. 0.]
>>> print(is2)
[2. 1. 0.]
Let's see the next sample. For the same shere as above, but let's see from the different position, ``(-1, 2, 0)``.
The view vector is assumed to be ``(2, 0, 0)``. In this geometry, the intersection is 1 point, (1, 2, 0).
The k-values are 1 (remmeber again the k-values are normalized by the view vector, in this case ``(2, 0, 0)``.).
>>> pos = [-1, 2, 0]
>>> view = [2, 0, 0]
>>> is3, is4, k = intersections_sphere(pos, view, cent, r)
>>> print(is3)
[1. 2. 0.]
>>> print(is4)
[1. 2. 0.]
What if one cannot see the sphere. Set the position ``(3, -1, 0)`` and looking toward ``(0, 1, 0)``.
``nan`` will be returned, while some warning messages are shown.
>>> pos = [3, -1, 0]
>>> view = [0, 1, 0]
>>> p = intersections_sphere(pos, view, cent, r)
>>> print(p)
(array([nan, nan, nan]), array([nan, nan, nan]), array([nan, nan]))
The last example is if the position is inside the sphere, e.g. at ``(1.5, 1, 0)``, looking toward ``(0, 1, 0)`` with a length of square root of 3.
In this case, k- becomes negative. But still you can calculate the intersection points.
>>> pos = [1.5, 1, 0]
>>> view = [0, np.sqrt(3), 0]
>>> is5, is6, k = intersections_sphere(pos, view, cent, r)
>>> print(is5)
[1.5 0.1339746 0. ]
>>> print(is6)
[1.5 1.8660254 0. ]
*Vectorized sample*
The above example can be vectorized, which has much higher performance.
>>> pos = [[-1, 1, 0], [-1, 2, 0,], [3, -1, 0], [1.5, 1, 0]] # (N=4, 3) shape
>>> view = [[1, 0, 0], [2, 0, 0], [0, 1, 0], [0, np.sqrt(3), 0]] # (N=4, 3) shape
>>> ism, isp, k = intersections_sphere(pos, view, cent, r)
>>> print(ism.shape) # (N=4, 3) shape
(4, 3)
>>> print(ism)
[[0. 1. 0. ]
[1. 2. 0. ]
[ nan nan nan]
[1.5 0.1339746 0. ]]
>>> print(isp) # (4, 3) shape
[[2. 1. 0. ]
[1. 2. 0. ]
[ nan nan nan]
[1.5 1.8660254 0. ]]
"""
k = lengths_to_intersections(pos, view, cent, radius) # Here ``k`` is (N, 2) shape
km, kp = np.atleast_2d(k).T # Now km, kp is scalar or (N,) array
p = np.atleast_2d(pos) # (N, 3)
v = np.atleast_2d(view) # (N, 3)
rm = np.squeeze(p + (v.T * km).T)
rp = np.squeeze(p + (v.T * kp).T)
return (rm, rp, k)
[docs]def lengths_to_intersections(pos, view, cent, radius):
r"""Return the lengths to the intersections.
::
* pos
view \
\ /----\
intersection1* \
/ \ cent \
| \ * |
\ \ /
\ \ /
\--*-/ sphere with radius
intersections2
:param pos: (3,) or (N, 3) shaped array of positions
:param view: (3,) or (N, 3) shaped array of viweing directions (or velocity vector)
:param cent: The center of the reference sphere
:param radius: The radius of the reference spehre
:return: An array of (2,) or (N, 2), where 2 refers to the length of intersection1 and intersection2
normalized by the ``view``. Always small for the first, large (or equal) for the second.
If no intersection, (np.nan, np.nan) is returned.
The returned is two floating point. To calculate the interesection1 and 2, you may use
:func:`intersections_sphere` function.
Samples. Let's consider the sphere cetnered at ``(1, 1, 0)`` with a radius of 1.
>>> cent = [1, 1, 0]
>>> r = 1
If one places at ``(-1, 1, 0)``, looking toward ``(1, 0, 0)``, then one sees the interesections
at ``(0, 1, 0)`` and ``(2, 1, 0)``; then ``k- = 1`` and ``k+ =3``. Note that the ``k``-values are
normalized by the length of view vector.
>>> pos = [-1, 1, 0]
>>> view = [1, 0, 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
>>> print(km, kp)
1.0 3.0
Let's see the next sample. For the same shere as above, but let's see from the different position, ``(-1, 2, 0)``.
The view vector is assumed to be ``(2, 0, 0)``. In this geometry, the intersection is 1 point, (1, 2, 0).
The k-values are 1 (remmeber again the k-values are normalized by the view vector, in this case ``(2, 0, 0)``.).
>>> pos = [-1, 2, 0]
>>> view = [2, 0, 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
>>> print(km, kp)
1.0 1.0
What if one cannot see the sphere. Set the position ``(3, -1, 0)`` and looking toward ``(0, 1, 0)``.
``nan`` will be returned, while some warning messages are shown.
>>> pos = [3, -1, 0]
>>> view = [0, 1, 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
The last example is if the position is inside the sphere, e.g. at ``(1.5, 1, 0)``, looking toward ``(0, 1, 0)`` with a length of square root of 3.
In this case, k- becomes negative. Still you can calculate the intersection points as above.
>>> pos = [1.5, 1, 0]
>>> view = [0, np.sqrt(3), 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
>>> print(km, kp)
-0.5 0.5
*Vectorized sample*
The above example can be vectorized, which has much higher performance.
>>> pos = [[-1, 1, 0], [-1, 2, 0,], [3, -1, 0], [1.5, 1, 0]]
>>> view = [[1, 0, 0], [2, 0, 0], [0, 1, 0], [0, np.sqrt(3), 0]]
>>> k = lengths_to_intersections(pos, view, cent, r)
>>> print(k[:, 0])
[ 1. 1. nan -0.5]
>>> print(k[:, 1])
[3. 1. nan 0.5]
"""
rvec = np.atleast_2d(pos) # (N, 3)
vvec = np.atleast_2d(view) # (N, 3)
ovec = np.atleast_2d(cent) # (N, 3)
rovec = rvec - ovec
v2 = (vvec ** 2).sum(axis=1) # (N,)
ro2 = (rovec ** 2).sum(axis=1) # (N,)
r2 = radius ** 2 # scalar
vdotro = (vvec * rovec).sum(axis=1) # (N,)
det = vdotro ** 2 - v2 * (ro2 - r2) # (N,)
rdet = np.sqrt(det) # For negative determinant, rdet becomes np.nan. It happens when no-intersections exist, and ok for most of the cases.
kminus = (-vdotro - rdet) / v2 # (N,)
kplus = (-vdotro + rdet) / v2 # (N,)
k = np.array([kminus, kplus]).transpose()
k = k.squeeze() # (N,) or (1,)
return k
[docs]def los_limb_sphere(pos, view, cent, radius):
r''' Return the intersection of LOS or limb.
This function returns the
intersection between ``view`` and the sphere.
If the configuration results in the intersection,
the result is the same as :func:`los_sphere`.
But if the intersection is not found
(view direction is into space not sphere),
this function returns the limb position of the sphere
in the same plane defined by the vectors of ``pos`` and ``view``.
::
* pos
| view
/---\ |
/ \ |
/ \ |
| * *----| No intersection, returning the neighboring surface
\ cent / |
\ /
\---/ sphere with radius
:param pos: Position of observer
:type pos: ``numpy.array``
:param view: Viewing direction
:type view: ``numpy.array``
:param cent: Center of the target sphere
:type cent: ``numpy.array``
:param radius: Radius of the sphere
:type radius: ``float``
:return: Position of the intersection or limb
as ``vector3d.Vector3d`` instance.
:rtype: ``vector3d.Vector3d`` instance or ``None``.
>>> print(los_limb_sphere([10,0,0], [-1,0,0], [-5,0,0], 3))
Vector3d( -2, 0, 0 )
>>> print(los_limb_sphere([10,0,0], [0,1,0], [-5,0,0], 3))
Vector3d( -4.4, 2.93939, 0 )
>>> print(los_limb_sphere([10,0,0], [-1,1,0], [-5,0,0], 3))
Vector3d( -4.4, 2.93939, 0 )
>>> print(los_limb_sphere([15,0,0], [-1,1,0], [0,0,0], 3))
Vector3d( 0.6, 2.93939, 0 )
*Note for developer*
The ``los_sphere`` is developed in vector3d.Vector3d class,
but ``los_limb_sphere`` is written basically using np.array class.
'''
los = los_sphere(pos, view, cent, radius)
if los is not None:
return los
### If no intersection was found.
# Convert to sphere-centered frame. All the calculation will be done.
rvec = np.array(pos, dtype=float) - np.array(cent, dtype=float)
kvec = np.array(view, dtype=float)
r2 = (rvec * rvec).sum() # r^2
k2 = (kvec * kvec).sum() # k^2
p2 = (rvec * kvec).sum() # r.k, inner product
R2 = radius * radius # R^2
# The position is expressed as tvec = a * rvec + b * kvec.
# b = sqrt ( (R^4 - R^2) / (p^4 - k^2 r^2) )
b = R2 * (R2 - r2) / (p2 * p2 - k2 * r2)
b = np.sqrt(b)
# a = (R^2 - p^2 * b) / r^2
a = (R2 - p2 * b) / r2
tvec = a * rvec + b * kvec
# Revert to the original frame
tvec = tvec + np.array(cent)
return Vector3d(tvec[0], tvec[1], tvec[2])
[docs]def isvisible(observer_xyz, target_xyz, sphere_center_xyz, radius):
""" Return if the target is visible from the observer if there is a sphere.
:param observer_xyz: Numpy array of observer location. Must be (3,) shape.
:param target_xyz: Numpy array of target location. Must be (3,) shape.
:param sphere_center_xyz: Numpy array of sphere center. Must be (3,) shape.
:param radius: Radius.
:returns: *True* if the target is visible.
*Use case*
It is commonly used if a spacecraft (observer) can see specific location (target) if a planetary body exist.
*Theory*
The observer can only see the specific target either when
- No intersection of its view angle and the sphere
- Within a virtual sphere centered at the observer with radius of the distance between the observer and the limb of the sphere
- But not inside the sphere
*Sample*
Suppose a sphere at (-5, 0, 0) with radius of 3. Observer is assumed to be at (5, 0, 0) in this example.
Point (-1.5, 0, 0) should be visible from the observer at (5, 0, 0).
>>> print(isvisible([5, 0, 0], [-1.5, 0, 0,], [-5., 0, 0], 3))
True
Point (-2.5, 0, 0) should not be visible from any observer because the point is inside the sphere
>>> print(isvisible([5, 0, 0], [-2.5, 0, 0], [-5, 0, 0], 3))
False
Point (10, 0, 0) should be visible, since the sphere is opposite direction for the observer.
>>> print(isvisible([5, 0, 0], [10, 0, 0], [-5, 0, 0], 3))
True
Point (-5, 3.001, 0) should not be visible, because it is below the horizon
>>> print(isvisible([5, 0, 0], [-5, 3.001, 0], [-5, 0, 0], 3))
False
Point (-10, 0, 0) should not be visible, because it is behind the sphere.
>>> print(isvisible([5, 0, 0], [-10, 0, 0], [-5, 0, 0], 3))
False
"""
observer_arr = np.array(observer_xyz)
target_arr = np.array(target_xyz)
sphere_center_arr = np.array(sphere_center_xyz)
if np.linalg.norm(target_arr - sphere_center_arr) < radius: # Target is in sphere
return False
view_vector = target_arr - observer_arr
intersection = los_sphere(observer_arr, view_vector, sphere_center_arr, radius, insphere='none')
if not intersection: # No intersection, meaning it is visible.
return True
distance_to_target = np.linalg.norm(target_arr - observer_arr)
visible_distance = np.sqrt(np.linalg.norm(sphere_center_arr - observer_arr) ** 2 - radius ** 2)
return distance_to_target <= visible_distance
[docs]def arevisible(observer_xyz, targets_xyz, sphere_center_xyz, radius):
""" Return if the target is visible from the observer if there is a sphere.
:param observer_xyz: Numpy array of observer location. Must be (3,) shape.
:param targets_xyz: Numpy array of target location. Must be (N, 3) shape.
:param sphere_center_xyz: Numpy array of sphere center. Must be (3,) shape.
:param radius: Radius.
:returns: *True* if the target is visible.
This is for evaluating multiple particles of visibility, i.e., parallelized version of :func:`isvisible`.
*Sample*
Suppose a sphere at (-5, 0, 0) with radius of 3. Observer is assumed to be at (5, 0, 0) in this example.
- Point (-1.5, 0, 0) should be visible from the observer at (5, 0, 0).
- Point (-2.5, 0, 0) should not be visible from any observer because the point is inside the sphere
- Point (10, 0, 0) should be visible, since the sphere is opposite direction for the observer.
- Point (-5, 3.001, 0) should not be visible, because it is below the horizon
- Point (-10, 0, 0) should not be visible, because it is behind the sphere.
>>> print(isvisible([5, 0, 0], [-1.5, 0, 0,], [-5., 0, 0], 3))
True
>>> print(isvisible([5, 0, 0], [-2.5, 0, 0], [-5, 0, 0], 3))
False
>>> print(isvisible([5, 0, 0], [10, 0, 0], [-5, 0, 0], 3))
True
>>> print(isvisible([5, 0, 0], [-5, 3.001, 0], [-5, 0, 0], 3))
False
>>> print(isvisible([5, 0, 0], [-10, 0, 0], [-5, 0, 0], 3))
False
>>> xyz = np.array([[-1.5, 0, 0], [-2.5, 0, 0], [10., 0, 0, ], [-5, 3.001, 0], [-10, 0, 0]])
>>> print(arevisible([5, 0, 0], xyz, [-5, 0, 0], 3))
[ True False True False False]
"""
observer = np.array(observer_xyz)
observer_arr = np.array(observer_xyz)[np.newaxis, :] # (1, 3) for broadcast
targets_arr = np.array(targets_xyz) # (N, 3)
sphere_center = np.array(sphere_center_xyz)
sphere_center_arr = np.array(sphere_center_xyz)[np.newaxis, :] # (1, 3) for broadcast
# In sphere validation
isoutside = (((targets_arr - sphere_center_arr) ** 2).sum(axis=1) >= radius ** 2)
# In cone or not
coneangle = np.arcsin(radius / np.linalg.norm(sphere_center - observer))
cone = _Cone(observer, sphere_center - observer, np.rad2deg(coneangle))
outofcone = ~cone.within(targets_arr) # To be visible, the axis should be outside the cone
# Inside sphere
distance_to_target = ((targets_arr - observer_arr) ** 2).sum(axis=1)
visible_distance = np.linalg.norm(sphere_center_arr - observer_arr) ** 2 - radius ** 2
insphere = (distance_to_target <= visible_distance)
return outofcone | (isoutside & insphere)
def _vec3d_to_array(v):
return np.array([v.x, v.y, v.z])
def _array_to_vec3d(a):
return Vector3d(a[0], a[1], a[2])
[docs]def sphere_limb_plane_projection(sphere_location, radius, resolution=360):
""" Return a (theta, phi) list of a sphere projected on a plane.
:param sphere_location: (x, y, z) of the sphere.
:param radius: The radius of the sphere
:returns: A pair (theta, phi). theta (and phi) is a numpy array with a size of (360,), or specified by the parameter resolution. These values are with the unit of radians.
theta ranges [0, pi], and phi ranges [-pi, pi].
:raises ValueError: if the sphere radius is larger than
the distance to the sphere.
(i.e., the observer at origin is inside the sphere)
It is used for example, to make a all sky map for ions
or ENAs, but also to plot the planet.
>>> theta, phi = sphere_limb_plane_projection([300, 0, 1000], 300)
>>> print(phi.shape)
(360,)
This gives the angles (theta and phi) of the planet at (300, 0, 1000) [km] with radius of 300 [km].
You get 360 points by default to be connected (or scatter plotted).
The returned angles are in radians. The 100'th sampling point is shown as
>>> print('{:.4f}'.format(phi[100]))
0.8962
>>> print('{:.4f}'.format(theta[100]))
0.3705
You can of course plot it
>>> plt.plot(phi, theta, '.') # doctest: +SKIP
*Theory*
It is a combination of frame rotations.
Sphere loated at (X, Y, Z) with radius R.
Observer at the origin.
Then, the projection to all sky map (theta, phi) is represented by
theta = acos(z)
phi = atan2(y, x)
(x, y, z) = ((cos PHI, -sin PHI, 0), (sin PHI, cos PHI, 0), (0, 0, 1)) . ((cos THETA, 0, sin THETA), (0, 1, 0), (-sin THETA, 0, cos THETA)) . (sin theta0 cos alpha, sin theta0 sin alpha, cos theta0)
PHI = atan2(Y, X)
THETA = acos(Z / L)
L = sqrt(X^2 + Y^2 + Z^2)
theta0 = asin(R / L)
alpha is the parameter to range from -pi to pi.
"""
X, Y, Z = sphere_location
R = radius
L = _np.linalg.norm([X, Y, Z])
if L <= R:
raise ValueError("Radius is bigger than the distance.")
theta0 = _np.arcsin(R / L) # Apparent angle of the sphere
THETA = _np.arccos(Z / L) # Sphere center in polar frame
PHI = _np.arctan2(Y, X)
# Parameter
alpha = _np.linspace(0, 2 * np.pi, resolution)
# The limb positions in the (x'', y'', z'') frame,
# where z'' is from the origin to the center of sphere.
xdd = _np.sin(theta0) * _np.cos(alpha)
ydd = _np.sin(theta0) * _np.sin(alpha)
zdd = _np.cos(theta0) * _np.ones([resolution])
rdd = _np.array([xdd, ydd, zdd]) # (3, 360)
# Conversion matrix from (x'', y'', z'') to (x, y, z), Euler angle
cos = _np.cos
sin = _np.sin
m1 = _np.array([[cos(THETA), 0, sin(THETA)],
[0, 1, 0],
[-sin(THETA), 0, cos(THETA)]])
_logger.info(m1)
m2 = _np.array([[cos(PHI), -sin(PHI), 0],
[sin(PHI), cos(PHI), 0],
[0, 0, 1]])
_logger.info(m2)
m = m2.dot(m1) # (3, 3) matrix
r = m.dot(rdd) # (3, 360) martix
_logger.info(rdd.shape)
_logger.info(r.shape)
theta = np.arccos(r[2, :])
phi = np.arctan2(r[1, :], r[0, :])
return (theta, phi)