''' An implementation of conversion to surface-interaction frame.
.. codeauthor:: Yoshifumi Futaana
Surface interaction frame (SIF) is defined as follows.
Inpinging angle, ti, is defined by the angle between
the negate of the inpinging velocity vector (-vi) and the surface normal, ns.
::
ti = angle(-vi, ns)
From the outgoing velocity vector (vo), two angles are defined.
The outgoing elevation angle, to, is defined as ::
to = angle(vo, ns)
The azimuthal angle, po, is the angle of the outgoing vector projected
on the surface from that of the incoming angle
This can be calculated by ::
|po| = angle( cross(vi, ns), cross(vo, ns) )
The +/- of po should also be considered, i.e. ::
Sign( ns dot cross(cross(vi, ns), cross(vo, ns) )
'''
import math
import logging
from irfpy.util.vector3d import Vector3d
[docs]def tosif(invec, outvec, normvec):
''' Calculate the SIF values from the given vectors.
Three input parameters (3-element array-like or Vector3d instance)
are needed. The instances will not be changed (const).
@param invec A impinging vector in arbitrary coordinate system.
@param outvec A outgoing vector in arbitrary coordinate system.
@param normvec A normal vector in arbitrary coordinate system.
@retval ti, to, po A tuple of three angles in degrees.
If the invec is anti-parallel to normvec, ti=0, po=0 is returned.
The condition ti > 90 or to > 90 is not acceptable
from a physical point of view, but the calculated values are
returned (with warning) since mathematically they can be calculated.
>>> t0, t1, p1 = tosif([-1, 0, -1],
... [-1/math.sqrt(2), 1/math.sqrt(2), 1], [0, 0, 1])
>>> print('%.1f' % t0)
45.0
>>> print('%.1f' % t1)
45.0
>>> print('%.1f' % p1)
-45.0
If the outgoing direction is opposite, p1 should be positive.
>>> t0, t1, p1 = tosif([-1, 0, -1],
... [-1/math.sqrt(2), -1/math.sqrt(2), 1], [0, 0, 1])
>>> print('%.1f' % t0)
45.0
>>> print('%.1f' % t1)
45.0
>>> print('%.1f' % p1)
45.0
In case of the flux coming from inside surface, warning is displayed.
>>> t0, t1, p1 = tosif([-1, 0, 1],
... [-1/math.sqrt(2), -1/math.sqrt(2), 1], [0, 0, 1])
>>> print('%.1f' % t0)
135.0
In case the flux is from antiparallel to the norm
>>> t0, t1, p1 = tosif([0, -1, 0], [1, 1, 1], [0, 1, 0])
>>> print('%.1f' % t0)
0.0
>>> print('%.1f' % p1)
0.0
'''
logger = logging.getLogger(__name__ + '.tosif')
vin = Vector3d(invec[0], invec[1], invec[2])
vin.normalize()
logger.debug(vin)
vout = Vector3d(outvec[0], outvec[1], outvec[2])
vout.normalize()
logger.debug(vout)
norm = Vector3d(normvec[0], normvec[1], normvec[2])
norm.normalize()
logger.debug(norm)
## Parallel flow?
if vin.angle(norm) == 0:
t0 = 180
else:
vin.negate() # vin is now -vin.
if vin.angle(norm) == 0:
t0 = 0
else:
t0 = norm.angle(vin) * 180. / math.pi
vin.negate() # vin is now vin.
if t0 > 90:
logger.warning('Input vector %s comes below surface (%d)? Normal %s.'
% (str(vin), t0, str(norm)))
t1 = norm.angle(vout) * 180. / math.pi
ixn = Vector3d()
ixn.cross(vin, norm)
oxn = Vector3d()
oxn.cross(vout, norm)
if ixn.lengthSquared() == 0 or oxn.lengthSquared() == 0:
p1 = 0
else:
p1 = ixn.angle(oxn) * 180. / math.pi
sp1 = Vector3d()
sp1.cross(ixn, oxn)
sign = sp1.dot(norm)
if sign < 0:
p1 = -p1
return t0, t1, p1
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')