""" Defines surface interaction frame
The definition:
Sufrace interaction frame is the Cartesian, defined by the two vector.
- The primary vector, which is usually the normal of surface, is referred by z_s.
- The secondary vector, which is usually the velocity vector of the projectile, defines (+x_s, z_s).
"""
from collections import namedtuple as _namedtuple
import numpy as _np
from irfpy.util import nptools as _nt
[docs]def conversion_matrix_surf2glob(primary_vector, secondary_vector):
""" Get the conversion matrix from the surface frame to global frame
:primary_vector: The primary vector, usually the surface normal vector,
expressed in the global coordinates.
:secondary_vector: The secondary vector, usually the incoming projectile velocity vector,
expressed in the global coordinates.
An example follows.
Assume the location of the north pole. The projectile as solar wind.
>>> norpol = [0, 0, 2400]
>>> swvel = [-400, 0, 0]
>>> matrix = conversion_matrix_surf2glob(norpol, swvel)
In this case, the conversion matrix is [[-1, 0, 0], [0, -1, 0], [0, 0, 1]].
>>> print(matrix) # doctest: +NORMALIZE_WHITESPACE
[[-1. 0. 0.]
[ 0. -1. 0.]
[ 0. 0. 1.]]
Let the location at the -45 deg in the noon meridian, with the projectile as solar wind.
>>> pos = [1, 0, -1]
>>> matrix = conversion_matrix_surf2glob(pos, swvel)
The surface normal vector ([0, 0, 1] in the surface coordinates) will be parallel to the position:
>>> print(matrix.dot([0, 0, 1]))
[ 0.70710678 0. -0.70710678]
The x axis in the surface coordiantes shall be z_glob=(0,0,0), while x_glob=(-.707, 0, -.707), and y_glob=0.
>>> print(matrix.dot([1, 0, 0]))
[-0.70710678 0. -0.70710678]
"""
zs = _np.array(primary_vector)
zs = zs / _np.linalg.norm(zs)
ys = _np.cross(zs, secondary_vector)
ys = ys / _np.linalg.norm(ys)
xs = _np.cross(ys, zs)
xs = xs / _np.linalg.norm(xs)
return _np.array([xs, ys, zs]).T
[docs]def conversion_matrix_s2g_normal_sun_reference(longitude, latitude):
""" Return the conversion matrix normal-sun reference.
The normal-sun-reference frame is, regardless of the name,
is widely used frame defined as follows.
- z axis as the local zenith
- x axis as opposite from the "Sun" (more precisely, [-1, 0, 0] in the global frame)
- y axis completes.
:param latitude: Latitude in degrees.
:param longitude: Logitude in degrees.
Let the observer at the north pole.
>>> cMat = conversion_matrix_s2g_normal_sun_reference(0, 90)
>>> print(cMat)
[[-1.000000e+00 0.000000e+00 6.123234e-17]
[ 0.000000e+00 -1.000000e+00 0.000000e+00]
[ 6.123234e-17 0.000000e+00 1.000000e+00]]
Let's put the observer at S45 Lon=0.
>>> cMat = conversion_matrix_s2g_normal_sun_reference(0, -45)
>>> print(cMat.dot([0, 0, 1]))
[ 0.70710678 0. -0.70710678]
>>> print(cMat.dot([1, 0, 0]))
[-0.70710678 0. -0.70710678]
Don't put the observer at the exact subsolar !
"""
lat = _np.radians(latitude)
lon = _np.radians(longitude)
prim = _np.array([_np.cos(lon) * _np.cos(lat),
_np.sin(lon) * _np.cos(lat),
_np.sin(lat)])
seco = _np.array([-1, 0, 0])
return conversion_matrix_surf2glob(prim, seco)
lonlattuple = _namedtuple('lonlattuple', ['lon', 'lat'])
meshtuple = _namedtuple('meshtuple', ['c', 'b', 'd', 'dOmega'])
""" Named tuple for mesh.
"""
[docs]def meshes(nLon=360, nLat=180, radian=False, flatten=False, centerLon=0):
""" Return 2d arrays, representing the center, bound, and delta for the longitude/latitude mesh in degrees.
:keyword nLon: Number of longitude bin. Longitude covers from -180 to 180 deg.
:keyword nLat: Number of latitude bin. Latitude covers from -90 to +90 deg.
:keyword radian: Return the value in radian if *True*.
:keyword flatten: Flatten the results into 1D.
:keyword centerLon: The center longitude. Usually it is 0.
:returns: A namedtuple for mesh grid (c, d, b, dOmega).
- c is the center of the grid, lonlattuple object.
- b is the boundary of the grid, lonlattuple object.
- d is the delta, lonlattuple object.
- dOmega is the steradian.
You can get the mesh as follows.
>>> cLon = meshes().c.lon # Obtain the center longtidue mesh.
>>> print(cLon.shape)
(360, 180)
>>> bLat = meshes().b.lat # Obtains the boundary latitude mesh.
>>> print(bLat.shape)
(361, 181)
You can also obtain the data at once.
>>> (cLon, cLat), (bLon, bLat), (dLon, dLat), dOmega = meshes()
>>> print(cLon.shape)
(360, 180)
>>> print(cLat.shape)
(360, 180)
>>> print(bLon.shape)
(361, 181)
>>> print(bLat.shape)
(361, 181)
>>> print(f'{dOmega.sum():.2f}') # Should be 4pi
12.57
Using the radian instead of degrees, and the array is flattened to 1D array.
>>> (cLon, cLat), (bLon, bLat), (dLon, dLat), dOmega = meshes(radian=True, flatten=True)
>>> print(cLon.shape)
(64800,)
>>> print(bLon.shape)
(65341,)
If you want to make the grid with longitude range (0, 360), give centerLon parameter.
>>> mesh = meshes(centerLon=180)
>>> mesh[0][0].min()
0.5
>>> mesh[0][0].max()
359.5
>>> import numpy as np
>>> mesh = meshes(centerLon=np.pi, radian=True)
>>> mesh[0][0].min()
0.008726646259971648
>>> mesh[0][0].max()
6.274458660919615
"""
if radian:
centerLon = _np.rad2deg(centerLon) # CenterLon shall be first converted to deg.
cLon, bLon, dLon = _nt.lingrids(-180, 180, nLon)
cLat, bLat, dLat = _nt.lingrids(-90, 90, nLat)
cLon = cLon + centerLon
bLon = bLon + centerLon
cLon, cLat = _np.meshgrid(cLon, cLat, indexing='ij')
bLon, bLat = _np.meshgrid(bLon, bLat, indexing='ij')
dLon, dLat = _np.meshgrid(dLon, dLat, indexing='ij')
dOmega = (_np.sin(_np.radians(bLat[1:, 1:])) - _np.sin(_np.radians(bLat[-1:, :-1]))) * _np.radians(dLon) # bLat forst argument does not matter.
if radian:
cLon = _np.radians(cLon)
cLat = _np.radians(cLat)
bLon = _np.radians(bLon)
bLat = _np.radians(bLat)
dLon = _np.radians(dLon)
dLat = _np.radians(dLat)
if flatten:
cLon = cLon.flatten()
cLat = cLat.flatten()
bLon = bLon.flatten()
bLat = bLat.flatten()
dLon = dLon.flatten()
dLat = dLat.flatten()
dOmega = dOmega.flatten()
c = lonlattuple(cLon, cLat)
b = lonlattuple(bLon, bLat)
d = lonlattuple(dLon, dLat)
mesh = meshtuple(c, b, d, dOmega)
return mesh