contrib/hkl.hs

   1 {-# LANGUAGE ScopedTypeVariables #-}
   2 {-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
   3 {-
   4     Copyright  : Copyright (C) 2014-2015 Synchrotron Soleil
   5     License    : GPL3+
   6
   7     Maintainer : picca@synchrotron-soleil.fr
   8     Stability  : Experimental
   9     Portability: GHC only?
  10 -}
  11
  12 import Prelude hiding (sqrt, sin, cos, (+), (-), (*), (**), (/))
  13 import qualified Prelude
  14
  15 import Numeric.LinearAlgebra (fromLists, Vector, Matrix,
  16                               ident, scalar, fromList,
  17                               (@>), (<>), vecdisp, disps, inv,
  18                               toList, dispf)
  19
  20 import Numeric.GSL.Root (root, RootMethod (Hybrids))
  21
  22 -- import Text.Printf (printf)
  23
  24 import Numeric.Units.Dimensional.Prelude (_0, _1, _2, nano, meter, degree,
  25                                           (*~), (/~), (+), (-), (*), (**), (/),
  26                                           Length, Angle, sin, cos, one, sqrt,
  27                                           Dimensionless, (*~~), (/~~))
  28
  29 import Options.Applicative hiding ((<>))
  30
  31 data Lattice = Cubic (Length Double) -- a = b = c, alpha = beta = gamma = 90
  32              | Tetragonal (Length Double) (Length Double) -- a = b != c, alpha = beta = gamma = 90
  33              | Orthorhombic (Length Double) (Length Double) (Length Double) -- a != b != c,  alpha = beta = gamma = 90
  34              | Rhombohedral (Length Double) (Angle Double) -- a = b = c, alpha = beta = gamma != 90
  35              | Hexagonal (Length Double) (Length Double) -- a = b != c, alpha = beta = 90, gamma = 120
  36              | Monoclinic (Length Double) (Length Double) (Length Double) (Angle Double) -- a != b != c, alpha = gamma = 90, beta != 90
  37              | Triclinic (Length Double) (Length Double) (Length Double) (Angle Double) (Angle Double) (Angle Double) -- a != b != c, alpha != beta != gamma != 90
  38                deriving (Show)
  39
  40 busing' :: Length Double -> Length Double -> Length Double -> Angle Double -> Angle Double-> Angle Double -> Matrix Double
  41 busing' a b c alpha beta gamma = fromLists [[b00 /~ (one / meter), b01/~ (one / meter), b02/~ (one / meter)],
  42                                             [0  , b11 /~ (one / meter), b12 /~ (one / meter)],
  43                                             [0  , 0  , b22 /~ (one / meter)]]
  44     where
  45       b00 = tau * sin alpha / (a * d)
  46       b01 = b11 / d * (cos alpha * cos beta - cos gamma)
  47       b02 = tmp / d * (cos gamma * cos alpha - cos beta)
  48       b11 = tau / (b * sin alpha)
  49       b12 = tmp / (sin beta * sin gamma) * (cos beta * cos gamma - cos alpha)
  50       b22 = tau / c
  51       d = sqrt(_1 - cos alpha ** _2 - cos beta ** _2 - cos gamma ** _2 + _2 * cos alpha * cos beta * cos gamma)
  52       tmp = b22 / sin alpha;
  53
  54 busing :: Lattice -> Matrix Double
  55 busing (Cubic a) = busing' a a a (90 *~ degree) (90 *~ degree) (90 *~ degree)
  56 busing (Tetragonal a c) = busing' a a c (90 *~ degree) (90 *~ degree) (90 *~ degree)
  57 busing (Orthorhombic a b c) = busing' a b c (90 *~ degree) (90 *~ degree) (90 *~ degree)
  58 busing (Rhombohedral a alpha)= busing' a a a alpha alpha alpha
  59 busing (Hexagonal a c) = busing' a a c (90 *~ degree) (90 *~ degree) (120 *~ degree)
  60 busing (Monoclinic a b c beta) = busing' a b c (90 *~ degree) beta (90 *~ degree)
  61 busing (Triclinic a b c alpha beta gamma) = busing' a b c alpha beta gamma
  62
  63 -- A Transformation which can be apply to a Vector of Double
  64 data Transformation = NoTransformation -- Doesn't transform the vector at all
  65                     | Rotation [Double] (Angle Double)
  66                     | UB Lattice
  67                     | Holder [Transformation]
  68
  69 tau :: Dimensionless Double
  70 tau = _1 -- 1 or 2*pi
  71
  72 crossprod :: Vector Double -> Matrix Double
  73 crossprod axis = fromLists [[ 0, -z,  y],
  74                             [ z,  0, -x],
  75                             [-y,  x,  0]]
  76     where
  77       x = axis @> 0
  78       y = axis @> 1
  79       z = axis @> 2
  80
  81 -- apply a transformation
  82 apply :: Transformation -> Vector Double -> Vector Double
  83 apply NoTransformation v = v
  84 apply (Rotation axis angle) v = (ident 3 Prelude.+ s Prelude.* q Prelude.+ c Prelude.* (q <> q)) <> v
  85     where
  86       ax = fromList axis
  87       c = scalar (1 Prelude.- cos angle /~ one)
  88       s = scalar (sin angle /~ one)
  89       q = crossprod ax
  90 apply (UB lattice) v = busing lattice <> v
  91 apply (Holder t) v = foldr apply v t
  92
  93 -- unapply a transformation
  94 unapply :: Vector Double -> Transformation -> Vector Double
  95 unapply v NoTransformation = v
  96 unapply v (Rotation axis angle) = apply (Rotation axis (_0 - angle)) v
  97 unapply v (UB lattice) =  inv (busing lattice) <> v
  98 unapply v (Holder t) = foldl unapply v t
  99
 100 -- source
 101 lambda :: Length Double
 102 lambda = 1.54 *~ nano meter
 103
 104 ki :: Vector Double
 105 ki = fromList [(tau / lambda) /~ (one / meter), 0, 0]
 106
 107 -- diffractometer
 108 data Diffractometer = Diffractometer [Angle Double -> Transformation] [Angle Double -> Transformation]
 109 data Mode = ModeHklE4CConstantPhi
 110
 111 e4c :: Diffractometer
 112 e4c = Diffractometer [omega, chi, phi, rx, ry, rz] [tth]
 113       where
 114         rx = Rotation [1, 0, 0]
 115         ry = Rotation [0, 1, 0]
 116         rz = Rotation [0, 0, 1]
 117         omega = Rotation [0, -1, 0]
 118         chi = Rotation [1, 0, 0]
 119         phi = Rotation [0, -1, 0]
 120         tth = Rotation [0, -1, 0]
 121
 122 computeHkl :: Diffractometer -> [Angle Double] -> Lattice -> Vector Double
 123 computeHkl (Diffractometer sample detector) values lattice =
 124     unapply q (Holder (zipWith ($) sample s ++ [UB lattice]))
 125         where
 126           (s, d) = splitAt (length sample) values
 127           kf = apply (Holder (zipWith ($) detector d)) ki
 128           q = kf Prelude.- ki
 129
 130 fromMode :: Mode -> [Double] -> [Angle Double] -> [Angle Double]
 131 fromMode ModeHklE4CConstantPhi fitted angles =
 132     newAngles
 133         where
 134           (_vs, _d) = splitAt 2 fitted
 135           (_cs, _) = splitAt 6 angles
 136           (_, _ccs) = splitAt 2 _cs
 137           newAngles = _vs *~~ degree ++ _ccs ++ _d *~~ degree
 138
 139 toMode :: Mode -> [Angle Double] -> [Double]
 140 toMode ModeHklE4CConstantPhi angles =
 141     v
 142         where
 143           (_s, _d) = splitAt 6 angles
 144           (_ss, _) = splitAt 2 _s
 145           v = (_ss ++ _d) /~~ degree
 146
 147 computeAngles' :: Diffractometer -> [Angle Double] -> Lattice -> Mode -> [Double] -> [Double] -> [Double]
 148 computeAngles' diffractometer angles lattice mode hkl fitted =
 149     toList (computeHkl diffractometer newAngles lattice Prelude.- fromList hkl)
 150         where
 151           newAngles = fromMode mode fitted angles
 152
 153 computeAngles :: Diffractometer -> [Angle Double] -> Lattice -> Mode -> [Double] -> ([Double], Matrix Double)
 154 computeAngles diffractometer angles lattice mode hkl =
 155     root Hybrids 1E-7 30 f guess
 156          where
 157            f = computeAngles' diffractometer angles lattice mode hkl
 158            guess = toMode mode angles
 159
 160 dispv :: Vector Double -> IO ()
 161 dispv = putStr . vecdisp (disps 2)
 162
 163 disp :: Matrix Double -> IO ()
 164 disp = putStr . dispf 3
 165
 166 -- command parsing
 167 data Command
 168     = Ca Double Double Double -- ca command
 169
 170 data Options
 171     = Options Command
 172
 173 withInfo :: Parser a -> String -> ParserInfo a
 174 withInfo opts desc = info (helper <*> opts) $ progDesc desc
 175
 176 parseCa :: Parser Command
 177 parseCa = Ca
 178     <$> argument auto (metavar "H")
 179     <*> argument auto (metavar "K")
 180     <*> argument auto (metavar "L")
 181
 182 parseCommand :: Parser Command
 183 parseCommand = subparser $
 184     command "ca"   (parseCa `withInfo` "compute angles for the given hkl")
 185
 186 parseOptions :: Parser Options
 187 parseOptions = Options <$> parseCommand
 188
 189 -- Actual program logic
 190 run :: Options -> IO ()
 191 run (Options cmd) =
 192     case cmd of
 193         Ca h k l-> do
 194            print (solution /~~ degree)
 195            dispv (computeHkl e4c solution lattice)
 196            disp path
 197                where
 198                  (sol, path) = computeAngles e4c angles lattice mode [h, k, l]
 199                  s = [30.0, 0.0, 0.0, 0.0, 10.0, 0.0]
 200                  d = [60.0]
 201                  angles = (s ++ d) *~~ degree
 202                  solution = fromMode mode sol angles
 203                  lattice = Cubic (1.54 *~ nano meter)
 204                  mode = ModeHklE4CConstantPhi
 205
 206 main :: IO ()
 207 main = run =<< execParser
 208     (parseOptions `withInfo` "Interact with hkl API")