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

```