contrib/haskell/src/Hkl/Transformation.hs

   1 module Hkl.Transformation
   2        ( Transformation (..)
   3        , apply
   4        , unapply
   5        )where
   6 {-
   7     Copyright  : Copyright (C) 2014-2016 Synchrotron Soleil
   8     License    : GPL3+
   9
  10     Maintainer : picca@synchrotron-soleil.fr
  11     Stability  : Experimental
  12     Portability: GHC only?
  13 -}
  14
  15 import Prelude hiding (sqrt, sin, cos, (+), (-), (*), (**), (/))
  16 import qualified Prelude
  17
  18 import Numeric.LinearAlgebra (fromLists, Vector, Matrix,
  19                               ident, scalar, fromList,
  20                               (@>), (<>), inv)
  21
  22 import Numeric.Units.Dimensional.Prelude (_0, (-), (/~),
  23                                           Angle, sin, cos, one)
  24 import Hkl.Lattice
  25
  26 -- A Transformation which can be apply to a Vector of Double
  27 data Transformation = NoTransformation -- Doesn't transform the vector at all
  28                     | Rotation [Double] (Angle Double)
  29                     | UB Lattice
  30                     | Holder [Transformation]
  31
  32 crossprod :: Vector Double -> Matrix Double
  33 crossprod axis = fromLists [[ 0, -z,  y],
  34                             [ z,  0, -x],
  35                             [-y,  x,  0]]
  36     where
  37       x = axis @> 0
  38       y = axis @> 1
  39       z = axis @> 2
  40
  41 -- apply a transformation
  42 apply :: Transformation -> Vector Double -> Vector Double
  43 apply NoTransformation v = v
  44 apply (Rotation axis angle) v = (ident 3 Prelude.+ s Prelude.* q Prelude.+ c Prelude.* (q <> q)) <> v
  45     where
  46       ax = fromList axis
  47       c = scalar (1 Prelude.- cos angle /~ one)
  48       s = scalar (sin angle /~ one)
  49       q = crossprod ax
  50 apply (UB lattice) v = busing lattice <> v
  51 apply (Holder t) v = foldr apply v t
  52
  53 -- unapply a transformation
  54 unapply :: Vector Double -> Transformation -> Vector Double
  55 unapply v NoTransformation = v
  56 unapply v (Rotation axis angle) = apply (Rotation axis (_0 - angle)) v
  57 unapply v (UB lattice) =  inv (busing lattice) <> v
  58 unapply v (Holder t) = foldl unapply v t