binoculars-ng/src/Hkl/Geometry.hs

   1 {-# LANGUAGE DeriveAnyClass    #-}
   2 {-# LANGUAGE DeriveGeneric     #-}
   3 {-# LANGUAGE FlexibleInstances #-}
   4 {-# LANGUAGE GADTs             #-}
   5 {-# LANGUAGE OverloadedStrings #-}
   6
   7 module Hkl.Geometry
   8        ( Factory(..)
   9        , Geometry(..)
  10        , GeometryState(..)
  11        , fixed
  12        , sixsMedHGisaxs
  13        , sixsMedVGisaxs
  14        , sixsUhvGisaxs
  15        , zaxis
  16        , withGeometry
  17        ) where
  18
  19 import           Data.Aeson            (FromJSON (..), ToJSON (..))
  20 import           Data.Tree             (Tree (..), foldTree)
  21 import qualified Data.Vector.Storable  as V
  22 import           Foreign               (ForeignPtr, Ptr, newForeignPtr, nullPtr,
  23                                         withForeignPtr)
  24 import           Foreign.C             (CDouble (..), newCString, withCString)
  25 import           GHC.Generics          (Generic)
  26 import           Numeric.LinearAlgebra (Vector)
  27
  28 import           Prelude               hiding (max, min)
  29
  30 import           Hkl.C.Hkl
  31 import           Hkl.Lattice
  32 import           Hkl.Orphan            ()
  33
  34 -------------
  35 -- Factory --
  36 -------------
  37
  38 data Factory = K6c | Fixe | Uhv | Mars | MedH | MedV | SoleilSiriusKappa
  39   deriving (Generic, FromJSON, ToJSON)
  40
  41 instance Show Factory where
  42   show K6c               = "K6C"
  43   show Fixe              = undefined
  44   show Uhv               = "ZAXIS"
  45   show Mars              = "SOLEIL MARS"
  46   show MedH              = "SOLEIL SIXS MED1+2"
  47   show MedV              = "SOLEIL SIXS MED2+3"
  48   show SoleilSiriusKappa = "SOLEIL SIRIUS KAPPA"
  49
  50 -- factoryFromString :: String -> Factory
  51 -- factoryFromString s
  52 --   | s == "K6C"  = K6c
  53 --   | s == undefined = Fixe
  54 --   | s == "ZAXIS" = Uhv
  55 --   | s == "SOLEIL MARS" = Mars
  56 --   | s == undefined = MedH
  57 --   | s == "SOLEIL SIXS MED2+3" = MedV
  58 --   | s == "SOLEIL SIRIUS KAPPA" = SoleilSiriusKappa
  59 --   | otherwise   = error $ "unknown diffractometer type:" ++ s
  60
  61 newFactory :: Factory -> IO (Ptr C'HklFactory)
  62 newFactory f = Foreign.C.withCString (show f) $ \cname -> c'hkl_factory_get_by_name cname nullPtr
  63
  64 --------------
  65 -- Geometry --
  66 --------------
  67
  68 data Transformation
  69   = NoTransformation
  70   | Rotation Double Double Double
  71   | Translation Double Double Double
  72   deriving (Generic, FromJSON, Show, ToJSON)
  73
  74 data Axis
  75   = Axis String Transformation Unit
  76   deriving (Generic, FromJSON, Show, ToJSON)
  77
  78 data GeometryState
  79     = GeometryState Double (Vector CDouble)
  80       deriving (Generic, FromJSON, Show, ToJSON)
  81
  82 data Geometry
  83   = Geometry'Custom (Tree Axis) (Maybe GeometryState)
  84   | Geometry'Factory Factory (Maybe GeometryState)
  85     deriving (Generic, FromJSON, Show, ToJSON)
  86
  87 fixed :: Geometry
  88 fixed
  89   = Geometry'Custom
  90     (Node (Axis "mu" (Rotation 0 0 1) Unit'Angle'Degree) [ Node (Axis "omega" (Rotation 0 (-1) 0) Unit'Angle'Degree) [] ])
  91     Nothing
  92
  93
  94 sixsMedHGisaxs :: Geometry
  95 sixsMedHGisaxs
  96   = Geometry'Custom
  97     ( Node (Axis "" NoTransformation Unit'NoUnit)
  98       [ Node (Axis "beta" (Rotation 0 (-1) 0) Unit'Angle'Degree) [Node (Axis "mu" (Rotation 0 0 1) Unit'Angle'Degree) [] ]
  99       , Node (Axis "eix" (Translation 0 0 (-1)) Unit'Length'MilliMeter) [Node (Axis "eiz" (Translation 0 1 0) Unit'Length'MilliMeter) [] ]
 100       ]
 101     )
 102     Nothing
 103
 104 sixsMedVGisaxs :: Geometry
 105 sixsMedVGisaxs
 106   = Geometry'Custom
 107     ( Node (Axis "" NoTransformation Unit'NoUnit)
 108       [ Node (Axis "beta" (Rotation 0 (-1) 0) Unit'Angle'Degree) [Node (Axis "mu" (Rotation 0  0 1) Unit'Angle'Degree) [Node (Axis "omega" (Rotation 0 (-1) 0) Unit'Angle'Degree) [] ] ]
 109       , Node (Axis "eix" (Translation 0 0 (-1)) Unit'Length'MilliMeter) [Node (Axis "eiz" (Translation 0 1 0) Unit'Length'MilliMeter) [] ]
 110       ]
 111     )
 112     Nothing
 113
 114 sixsUhvGisaxs :: Geometry
 115 sixsUhvGisaxs
 116   = Geometry'Custom
 117     ( Node (Axis "" NoTransformation Unit'NoUnit)
 118       [ Node (Axis "mu" (Rotation 0 0 1) Unit'Angle'Degree) [ Node (Axis "omega" (Rotation 0 (-1) 0) Unit'Angle'Degree) [] ]
 119       , Node (Axis "eix" (Translation 0 (-1) 0) Unit'Length'MilliMeter) [Node (Axis "eiz" (Translation 0 0 1) Unit'Length'MilliMeter) []]
 120       ]
 121     )
 122     Nothing
 123
 124 zaxis :: Geometry
 125 zaxis
 126   = Geometry'Custom
 127     ( Node (Axis "mu" (Rotation 0 0 1) Unit'Angle'Degree)
 128       [ Node (Axis "omega" (Rotation 0 (-1) 0) Unit'Angle'Degree) []
 129       , Node (Axis "delta" (Rotation 0 (-1) 0) Unit'Angle'Degree) [ Node (Axis "gamma" (Rotation 0 0 1) Unit'Angle'Degree) [] ]
 130       ]
 131     )
 132     Nothing
 133
 134 pokeGeometry :: ForeignPtr C'HklGeometry -> GeometryState -> IO ()
 135 pokeGeometry fptr (GeometryState lw vs) =
 136   withForeignPtr fptr $ \ptr -> do
 137   -- set the source
 138   let wavelength = CDouble lw
 139   c'hkl_geometry_wavelength_set ptr wavelength c'HKL_UNIT_USER nullPtr
 140
 141   -- set the axes
 142   let n = toEnum . V.length $ vs
 143   V.unsafeWith vs $ \values ->
 144     c'hkl_geometry_axis_values_set ptr values n c'HKL_UNIT_USER nullPtr
 145
 146 newGeometry :: Geometry -> IO (ForeignPtr C'HklGeometry)
 147 newGeometry (Geometry'Custom g ms)
 148   = do
 149   gPtr <- c'hkl_geometry_new nullPtr nullPtr
 150   mapM_ (addHolder gPtr) (axes g)
 151   fptr <- newForeignPtr p'hkl_geometry_free gPtr
 152   case ms of
 153     Nothing  -> return ()
 154     (Just s) -> pokeGeometry fptr s
 155   return fptr
 156     where
 157       addHolder :: Ptr C'HklGeometry -> [Axis] -> IO ()
 158       addHolder gPtr axs = do
 159         hPtr <- c'hkl_geometry_add_holder gPtr
 160         mapM_ (addAxis hPtr) axs
 161
 162       addAxis :: Ptr C'HklHolder -> Axis -> IO ()
 163       addAxis h (Axis n t u) = case t of
 164         NoTransformation -> return ()
 165         Rotation x y z -> do
 166           c'n <- newCString n
 167           c'hkl_holder_add_rotation h c'n (CDouble x) (CDouble y) (CDouble z) (unitToHklUnit u)
 168         Translation x y z -> do
 169           c'n <- newCString n
 170           c'hkl_holder_add_translation h c'n (CDouble x) (CDouble y) (CDouble z) (unitToHklUnit u)
 171
 172       axes :: Tree Axis -> [[Axis]]
 173       axes g' = foldTree (\x xsss -> if null xsss then [[x]] else concatMap (map (x:)) xsss) g'
 174 newGeometry (Geometry'Factory f ms)
 175     = do fptr <- newForeignPtr p'hkl_geometry_free =<< c'hkl_factory_create_new_geometry =<< newFactory f
 176          case ms of
 177            Nothing  -> return ()
 178            (Just s) -> pokeGeometry fptr s
 179          return fptr
 180
 181 withGeometry :: Geometry -> (Ptr C'HklGeometry -> IO r) -> IO r
 182 withGeometry g f
 183     = do fptr <- newGeometry g
 184          withForeignPtr fptr f