
##############################################################################
#
# Copyright (c) 2003-2020 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
# Development from 2019 by School of Earth and Environmental Sciences
#
##############################################################################

from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2020 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

"""
transformations

:var __author__: name of author
:var __copyright__: copyrights
:var __license__: licence agreement
:var __url__: url entry point on documentation
:var __version__: version
:var __date__: date of the version
:var DEG: unit of degree
:var RAD: unit of radiant
"""

__author__="Lutz Gross, l.gross@uq.edu.au"

import numpy
import math

_TYPE=numpy.float64
DEG=math.pi/180.
RAD=1.
class Transformation(object):
   """
   General class to define an affine transformation *x->Ax+b*.
   """
   def __init__(self):
       """
       Creates a linear transformation.
       """
       pass

   def __call__(self,x=numpy.zeros((3,))):
       """
       Applies transformation to ``x``.
       """
       raise NotImplementeError()

class Translation(Transformation):
    """
    Defines a translation *x->x+b*.
    """
    def __init__(self,b=numpy.zeros((3,),dtype=_TYPE)):
       """
       Creates the linear transformation *x->x+b*.
       """
       super(Translation, self).__init__()
       self.__b=numpy.array(b,_TYPE)

    def __call__(self,x=numpy.zeros((3,))):
       """
       Applies translation to ``x``.
       """
       return numpy.array(x,_TYPE)+self.__b

class Rotatation(Transformation):
    """
    Defines a rotation.
    """
    def __init__(self,axis=numpy.ones((3,),dtype=_TYPE),point=numpy.zeros((3,),dtype=_TYPE),angle=0.*RAD):
       """
       Creates a rotation using an axis and a point on the axis.
       """
       self.__axis=numpy.array(axis,dtype=_TYPE)
       self.__point=numpy.array(point,dtype=_TYPE)
       lax=numpy.dot(self.__axis,self.__axis)
       if not lax>0:
          raise ValueError("points must be distinct.")
       self.__axis/=math.sqrt(lax)
       self.__angle=float(angle)

    def __call__(self,x=numpy.zeros((3,))):
       """
       Applies the rotation to ``x``.
       """
       x=numpy.array(x,_TYPE)
       z=x-self.__point
       z0=numpy.dot(z,self.__axis)
       z_per=z-z0*self.__axis
       lz_per=numpy.dot(z_per,z_per)
       if lz_per>0:
         axis1=z_per/math.sqrt(lz_per)
         axis2=_cross(axis1,self.__axis)
         lax2=numpy.dot(axis2,axis2)
         if lax2>0:
            axis2/=math.sqrt(lax2)
            return z0*self.__axis+math.sqrt(lz_per)*(math.cos(self.__angle)*axis1-math.sin(self.__angle)*axis2)+self.__point
         else:
            return x
       else:
         return x

def _cross(x, y):
    """
    Returns the cross product of ``x`` and ``y``.
    """
    return numpy.array([x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] * y[0]], _TYPE)

class Dilation(Transformation):
    """
    Defines a dilation.
    """
    def __init__(self,factor=1.,center=numpy.zeros((3,),dtype=_TYPE)):
       """
       Creates a dilation with a center and a given expansion/contraction
       factor.
       """
       if not abs(factor)>0:
          raise ValueError("factor must be non-zero.")
       self.__factor=factor
       self.__center=numpy.array(center,dtype=_TYPE)

    def __call__(self,x=numpy.zeros((3,))):
       """
       Applies dilation to ``x``.
       """
       x=numpy.array(x,_TYPE)
       return self.__factor*(x-self.__center)+self.__center

class Reflection(Transformation):
    """
    Defines a reflection on a plane.
    """
    def __init__(self,normal=numpy.ones((3,),dtype=_TYPE),offset=0.):
       """
       Defines a reflection on a plane defined in normal form.
       """
       self.__normal=numpy.array(normal,dtype=_TYPE)
       ln=math.sqrt(numpy.dot(self.__normal,self.__normal))
       if not ln>0.:
          raise ValueError("normal must have positive length.")
       self.__normal/=ln
       if isinstance(offset,float) or isinstance(offset,int):
          self.__offset=offset/ln
       else:
          self.__offset=numpy.dot(numpy.array(offset,dtype=_TYPE),self.__normal)

    def __call__(self,x=numpy.zeros((3,))):
       """
       Applies reflection to ``x``.
       """
       x=numpy.array(x,_TYPE)
       return x - 2*(numpy.dot(x,self.__normal)-self.__offset)*self.__normal


[((3484, 3590), 'numpy.array', 'numpy.array', (['[x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] *\n    y[0]]', '_TYPE'], {}), '([x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y\n    [1] - x[1] * y[0]], _TYPE)\n', (3495, 3590), False, 'import numpy\n'), ((1623, 1640), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {}), '((3,))\n', (1634, 1640), False, 'import numpy\n'), ((1851, 1881), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {'dtype': '_TYPE'}), '((3,), dtype=_TYPE)\n', (1862, 1881), False, 'import numpy\n'), ((2015, 2036), 'numpy.array', 'numpy.array', (['b', '_TYPE'], {}), '(b, _TYPE)\n', (2026, 2036), False, 'import numpy\n'), ((2061, 2078), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {}), '((3,))\n', (2072, 2078), False, 'import numpy\n'), ((2286, 2315), 'numpy.ones', 'numpy.ones', (['(3,)'], {'dtype': '_TYPE'}), '((3,), dtype=_TYPE)\n', (2296, 2315), False, 'import numpy\n'), ((2321, 2351), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {'dtype': '_TYPE'}), '((3,), dtype=_TYPE)\n', (2332, 2351), False, 'import numpy\n'), ((2472, 2502), 'numpy.array', 'numpy.array', (['axis'], {'dtype': '_TYPE'}), '(axis, dtype=_TYPE)\n', (2483, 2502), False, 'import numpy\n'), ((2522, 2553), 'numpy.array', 'numpy.array', (['point'], {'dtype': '_TYPE'}), '(point, dtype=_TYPE)\n', (2533, 2553), False, 'import numpy\n'), ((2564, 2599), 'numpy.dot', 'numpy.dot', (['self.__axis', 'self.__axis'], {}), '(self.__axis, self.__axis)\n', (2573, 2599), False, 'import numpy\n'), ((2695, 2709), 'math.sqrt', 'math.sqrt', (['lax'], {}), '(lax)\n', (2704, 2709), False, 'import math\n'), ((2768, 2785), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {}), '((3,))\n', (2779, 2785), False, 'import numpy\n'), ((2857, 2878), 'numpy.array', 'numpy.array', (['x', '_TYPE'], {}), '(x, _TYPE)\n', (2868, 2878), False, 'import numpy\n'), ((2912, 2937), 'numpy.dot', 'numpy.dot', (['z', 'self.__axis'], {}), '(z, self.__axis)\n', (2921, 2937), False, 'import numpy\n'), ((2981, 3004), 'numpy.dot', 'numpy.dot', (['z_per', 'z_per'], {}), '(z_per, z_per)\n', (2990, 3004), False, 'import numpy\n'), ((3698, 3728), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {'dtype': '_TYPE'}), '((3,), dtype=_TYPE)\n', (3709, 3728), False, 'import numpy\n'), ((3974, 4006), 'numpy.array', 'numpy.array', (['center'], {'dtype': '_TYPE'}), '(center, dtype=_TYPE)\n', (3985, 4006), False, 'import numpy\n'), ((4031, 4048), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {}), '((3,))\n', (4042, 4048), False, 'import numpy\n'), ((4116, 4137), 'numpy.array', 'numpy.array', (['x', '_TYPE'], {}), '(x, _TYPE)\n', (4127, 4137), False, 'import numpy\n'), ((4314, 4343), 'numpy.ones', 'numpy.ones', (['(3,)'], {'dtype': '_TYPE'}), '((3,), dtype=_TYPE)\n', (4324, 4343), False, 'import numpy\n'), ((4461, 4493), 'numpy.array', 'numpy.array', (['normal'], {'dtype': '_TYPE'}), '(normal, dtype=_TYPE)\n', (4472, 4493), False, 'import numpy\n'), ((4878, 4895), 'numpy.zeros', 'numpy.zeros', (['(3,)'], {}), '((3,))\n', (4889, 4895), False, 'import numpy\n'), ((4965, 4986), 'numpy.array', 'numpy.array', (['x', '_TYPE'], {}), '(x, _TYPE)\n', (4976, 4986), False, 'import numpy\n'), ((2154, 2175), 'numpy.array', 'numpy.array', (['x', '_TYPE'], {}), '(x, _TYPE)\n', (2165, 2175), False, 'import numpy\n'), ((3118, 3141), 'numpy.dot', 'numpy.dot', (['axis2', 'axis2'], {}), '(axis2, axis2)\n', (3127, 3141), False, 'import numpy\n'), ((4513, 4552), 'numpy.dot', 'numpy.dot', (['self.__normal', 'self.__normal'], {}), '(self.__normal, self.__normal)\n', (4522, 4552), False, 'import numpy\n'), ((3045, 3062), 'math.sqrt', 'math.sqrt', (['lz_per'], {}), '(lz_per)\n', (3054, 3062), False, 'import math\n'), ((3180, 3195), 'math.sqrt', 'math.sqrt', (['lax2'], {}), '(lax2)\n', (3189, 3195), False, 'import math\n'), ((4806, 4838), 'numpy.array', 'numpy.array', (['offset'], {'dtype': '_TYPE'}), '(offset, dtype=_TYPE)\n', (4817, 4838), False, 'import numpy\n'), ((5007, 5034), 'numpy.dot', 'numpy.dot', (['x', 'self.__normal'], {}), '(x, self.__normal)\n', (5016, 5034), False, 'import numpy\n'), ((3230, 3247), 'math.sqrt', 'math.sqrt', (['lz_per'], {}), '(lz_per)\n', (3239, 3247), False, 'import math\n'), ((3249, 3271), 'math.cos', 'math.cos', (['self.__angle'], {}), '(self.__angle)\n', (3257, 3271), False, 'import math\n'), ((3278, 3300), 'math.sin', 'math.sin', (['self.__angle'], {}), '(self.__angle)\n', (3286, 3300), False, 'import math\n')]