Merge pull request #33 from Baharis/property_at_focus

Property at focus

docs/conf.py

@@ -22,7 +22,7 @@
 author = u'Daniel Tchoń'
 
 # The full version, including alpha/beta/rc tags
-release = '0.2.1'
+release = '0.2.2'
 
 
 # -- General configuration ---------------------------------------------------

hikari/dataframes/hkl.py

@@ -571,7 +571,7 @@ def fill(self, radius=2.0):
         def _make_hkl_ball(i=max_index):
             hkl_grid = np.mgrid[-i:i:2j*i+1j, -i:i:2j*i+1j, -i:i:2j*i+1j]
             hkls = np.stack(hkl_grid, -1).reshape(-1, 3)
-            xyz = np.matrix((self.a_w, self.b_w, self.c_w))
+            xyz = np.array((self.a_w, self.b_w, self.c_w))
             return hkls[lin.norm(hkls @ xyz, axis=1) <= radius]
         hkl = _make_hkl_ball()
 

hikari/scripts/angular_explorer.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 from numpy import linalg as lin
+import pandas as pd
 
 from hikari.dataframes import HklFrame, LstFrame
 from hikari.symmetry import SG, Group
@@ -184,7 +185,7 @@ def _draw_map(self, artist, extension):
         a.x_axis = self.hkl_frame.a_w / lin.norm(self.hkl_frame.a_w)
         a.y_axis = self.hkl_frame.b_w / lin.norm(self.hkl_frame.b_w)
         a.z_axis = self.hkl_frame.c_w / lin.norm(self.hkl_frame.c_w)
-        a.focus = self.focus
+        a.focus = self.focus_cartesian
         a.heat_limits = self.prop_limits
         a.heat_palette = self.axis
         a.histogram = self.histogram
@@ -243,7 +244,7 @@ def lg(self):
         return self._lg
 
     @property
-    def focus(self):
+    def focus_cartesian(self):
         _focus = []
         if self.orientation is not None:
             if len(self.orientation.shape) == 1:
@@ -259,6 +260,28 @@ def focus(self):
                     _focus.append(v / lin.norm(v))
         return _focus
 
+    @property
+    def focus_spherical_closest(self):
+        _focus = []
+        if self.orientation is not None:
+            for focus_vector in self.focus_cartesian:
+                _, th, ph = np.rad2deg(cart2sph(*focus_vector))
+                th_closest = min(self.th_range, key=lambda x: abs(x - th))
+                ph_closest = min(self.ph_range, key=lambda x: abs(x - ph))
+                _focus.append([1.0, th_closest, ph_closest])
+        return _focus
+
+    @property
+    def prop_at_focus_closest(self):
+        df = pd.DataFrame({'ph': self.data_dict['ph'],
+                           'th': self.data_dict['th'],
+                           'p': self.data_dict[self.property_name]})
+        properties_ = []
+        for _, th, ph in self.focus_spherical_closest:
+            row = df.loc[(df['th'] == th) & (df['ph'] == ph)].iloc[0]
+            properties_.append(row['p'])
+        return properties_
+
     @property
     def th_range(self):
         return self.th_limits.arange(step=self.angle_res)
@@ -323,6 +346,15 @@ def descriptive_statistics_string(self):
             f'# avg ={avg_p:8.5f}\n'
         return s
 
+    @property
+    def prop_at_focus_string(self):
+        s = f'# value of {self.property_name} at focus points:\n'
+        for focus_, prop_ in zip(self.focus_spherical_closest,
+                                 self.prop_at_focus_closest):
+            _, th, ph = focus_
+            s += f'# val ={prop_:8.5f} at th ={th :6.1f} ph ={ph :6.1f}\n'
+        return s + '\n'
+
 
 class AngularPotencyExplorer(AngularPropertyExplorer):
     hkl_is_read_not_generated = False
@@ -358,6 +390,8 @@ def explore(self):
             lst.write('\n')
 
         np.savetxt(dat_path, potency_mesh)
+        if self.orientation is not None:
+            lst.write(self.prop_at_focus_string)
         lst.write(self.descriptive_statistics_string)
         lst.close()
 
@@ -408,6 +442,8 @@ def explore(self):
             lst.write('\n')
 
         np.savetxt(dat_path, r1_mesh)
+        if self.orientation is not None:
+            lst.write(self.prop_at_focus_string)
         lst.write(self.descriptive_statistics_string)
         lst.close()
 

setup.py

@@ -25,7 +25,7 @@
 
 setup(
     name='hikari-toolkit',
-    version='0.2.1',
+    version='0.2.2',
     author='Daniel Tchoń',
     author_email='[email protected]',
     packages=find_packages(exclude=('legacy', )),

test/test_scripts.py

@@ -2,7 +2,9 @@
 import tempfile
 import unittest
 
-from hikari.scripts import calculate_similarity_indices
+import numpy as np
+
+from hikari.scripts import calculate_similarity_indices, potency_map
 
 
 nacl_cif_path = str(pathlib.Path(__file__).parent.joinpath('NaCl.cif'))
@@ -65,5 +67,22 @@ def test_calculate_similarity_index_for_negative_adps(self):
             self.assertAlmostEqual(float(last_line_contents[-3]), 50.0)
 
 
+class TestHklPotencyScripts(unittest.TestCase):
+    temp_dir = tempfile.TemporaryDirectory()
+    hkl_path = str(pathlib.Path(temp_dir.name) / 'potency.hkl')
+
+    def test_potency_map_simple(self):
+        potency_map(a=10, b=10, c=10, al=90, be=90, ga=90, space_group='P1',
+                    path=self.hkl_path, output_quality=2, wavelength='MoKa')
+
+    def test_potency_map_with_focus(self):
+        ori = np.array([[-0.08263, +0.00536, -0.03667],
+                        [-0.01671, -0.03679, -0.03238],
+                        [+0.06030, +0.02438, -0.04088]])
+        potency_map(a=10, b=10, c=10, al=90, be=90, ga=90, space_group='Pm-3m',
+                    path=self.hkl_path, output_quality=2, wavelength='MoKa',
+                    orientation=ori)
+
+
 if __name__ == '__main__':
     unittest.main()