Source code for irfpy.util.surfaceframe

''' 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')