contrib/haskell/src/Hkl/Diffractometer.hs

   1 module Hkl.Diffractometer
   2        ( e4c
   3        , Mode (..)
   4        , computeAngles
   5        , computeHkl
   6        , fromMode
   7        ) where
   8 {-
   9     Copyright  : Copyright (C) 2014-2016 Synchrotron Soleil
  10     License    : GPL3+
  11
  12     Maintainer : picca@synchrotron-soleil.fr
  13     Stability  : Experimental
  14     Portability: GHC only?
  15 -}
  16
  17 import Data.Tree
  18 import Prelude hiding (sqrt, sin, cos, (+), (-), (*), (**), (/))
  19 import qualified Prelude
  20
  21 import Numeric.LinearAlgebra (Vector, Matrix,
  22                               fromList,
  23                               toList)
  24
  25 import Numeric.GSL.Root (root, RootMethod (Hybrids))
  26
  27 -- import Text.Printf (printf)
  28
  29 import Numeric.Units.Dimensional.Prelude (degree,
  30                                           Angle,
  31                                           (*~~), (/~~))
  32
  33 import Hkl.Lattice
  34 import Hkl.Source
  35 import Hkl.Transformation
  36
  37 -- diffractometer
  38 -- data Diffractometer = Diffractometer [Angle Double -> Transformation] [Angle Double -> Transformation]
  39 newtype Diffractometer = Tree (Angle Double -> Transformation)
  40 data Mode = ModeHklE4CConstantPhi
  41
  42 e4c :: Tree (Angle Double -> Transformation)
  43 e4c = Node base
  44       [ Node omega
  45         [ Node chi
  46           [ Node phi
  47             [ Node rx
  48               [ Node ry
  49                 [ Node rz [] ]]]]]
  50       , Node tth []
  51       ]
  52   where
  53     base = NoTransformation
  54     rx = Rotation [1, 0, 0]
  55     ry = Rotation [0, 1, 0]
  56     rz = Rotation [0, 0, 1]
  57     omega = Rotation [0, -1, 0]
  58     chi = Rotation [1, 0, 0]
  59     phi = Rotation [0, -1, 0]
  60     tth = Rotation [0, -1, 0]
  61
  62 getSample :: Diffractometer -> [Angle Double -> Transformation]
  63 getSample (Tree _ [x:xs]) = flatten x
  64
  65 getDetector :: Diffractometer -> [Angle Double -> Transformation]
  66 getDetector (Tree _ [x:xs]) = flatten xs
  67
  68 computeHkl :: Diffractometer -> [Angle Double] -> Lattice -> Vector Double
  69 computeHkl diffractometer values lattice =
  70     unapply q (Holder (zipWith ($) sample s ++ [UB lattice]))
  71         where
  72           sample = getSample diffractometer
  73           detector = getDetector diffractometer
  74           (s, d) = splitAt (length sample) values
  75           kf = apply (Holder (zipWith ($) detector d)) ki
  76           q = kf Prelude.- ki
  77
  78 fromMode :: Mode -> [Double] -> [Angle Double] -> [Angle Double]
  79 fromMode ModeHklE4CConstantPhi fitted angles =
  80     newAngles
  81         where
  82           (_vs, _d) = splitAt 2 fitted
  83           (_cs, _) = splitAt 6 angles
  84           (_, _ccs) = splitAt 2 _cs
  85           newAngles = _vs *~~ degree ++ _ccs ++ _d *~~ degree
  86
  87 toMode :: Mode -> [Angle Double] -> [Double]
  88 toMode ModeHklE4CConstantPhi angles =
  89     v
  90         where
  91           (_s, _d) = splitAt 6 angles
  92           (_ss, _) = splitAt 2 _s
  93           v = (_ss ++ _d) /~~ degree
  94
  95 computeAngles' :: Diffractometer -> [Angle Double] -> Lattice -> Mode -> [Double] -> [Double] -> [Double]
  96 computeAngles' diffractometer angles lattice mode hkl fitted =
  97     toList (computeHkl diffractometer newAngles lattice Prelude.- fromList hkl)
  98         where
  99           newAngles = fromMode mode fitted angles
 100
 101 computeAngles :: Diffractometer -> [Angle Double] -> Lattice -> Mode -> [Double] -> ([Double], Matrix Double)
 102 computeAngles diffractometer angles lattice mode hkl =
 103     root Hybrids 1E-7 30 f guess
 104          where
 105            f = computeAngles' diffractometer angles lattice mode hkl
 106            guess = toMode mode angles
 107