_rotation_groups.py 4.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  1. import numpy as np
  2. from scipy.constants import golden as phi
  3. def icosahedral(cls):
  4. g1 = tetrahedral(cls).as_quat()
  5. a = 0.5
  6. b = 0.5 / phi
  7. c = phi / 2
  8. g2 = np.array([[+a, +b, +c, 0],
  9. [+a, +b, -c, 0],
  10. [+a, +c, 0, +b],
  11. [+a, +c, 0, -b],
  12. [+a, -b, +c, 0],
  13. [+a, -b, -c, 0],
  14. [+a, -c, 0, +b],
  15. [+a, -c, 0, -b],
  16. [+a, 0, +b, +c],
  17. [+a, 0, +b, -c],
  18. [+a, 0, -b, +c],
  19. [+a, 0, -b, -c],
  20. [+b, +a, 0, +c],
  21. [+b, +a, 0, -c],
  22. [+b, +c, +a, 0],
  23. [+b, +c, -a, 0],
  24. [+b, -a, 0, +c],
  25. [+b, -a, 0, -c],
  26. [+b, -c, +a, 0],
  27. [+b, -c, -a, 0],
  28. [+b, 0, +c, +a],
  29. [+b, 0, +c, -a],
  30. [+b, 0, -c, +a],
  31. [+b, 0, -c, -a],
  32. [+c, +a, +b, 0],
  33. [+c, +a, -b, 0],
  34. [+c, +b, 0, +a],
  35. [+c, +b, 0, -a],
  36. [+c, -a, +b, 0],
  37. [+c, -a, -b, 0],
  38. [+c, -b, 0, +a],
  39. [+c, -b, 0, -a],
  40. [+c, 0, +a, +b],
  41. [+c, 0, +a, -b],
  42. [+c, 0, -a, +b],
  43. [+c, 0, -a, -b],
  44. [0, +a, +c, +b],
  45. [0, +a, +c, -b],
  46. [0, +a, -c, +b],
  47. [0, +a, -c, -b],
  48. [0, +b, +a, +c],
  49. [0, +b, +a, -c],
  50. [0, +b, -a, +c],
  51. [0, +b, -a, -c],
  52. [0, +c, +b, +a],
  53. [0, +c, +b, -a],
  54. [0, +c, -b, +a],
  55. [0, +c, -b, -a]])
  56. return cls.from_quat(np.concatenate((g1, g2)))
  57. def octahedral(cls):
  58. g1 = tetrahedral(cls).as_quat()
  59. c = np.sqrt(2) / 2
  60. g2 = np.array([[+c, 0, 0, +c],
  61. [0, +c, 0, +c],
  62. [0, 0, +c, +c],
  63. [0, 0, -c, +c],
  64. [0, -c, 0, +c],
  65. [-c, 0, 0, +c],
  66. [0, +c, +c, 0],
  67. [0, -c, +c, 0],
  68. [+c, 0, +c, 0],
  69. [-c, 0, +c, 0],
  70. [+c, +c, 0, 0],
  71. [-c, +c, 0, 0]])
  72. return cls.from_quat(np.concatenate((g1, g2)))
  73. def tetrahedral(cls):
  74. g1 = np.eye(4)
  75. c = 0.5
  76. g2 = np.array([[c, -c, -c, +c],
  77. [c, -c, +c, +c],
  78. [c, +c, -c, +c],
  79. [c, +c, +c, +c],
  80. [c, -c, -c, -c],
  81. [c, -c, +c, -c],
  82. [c, +c, -c, -c],
  83. [c, +c, +c, -c]])
  84. return cls.from_quat(np.concatenate((g1, g2)))
  85. def dicyclic(cls, n, axis=2):
  86. g1 = cyclic(cls, n, axis).as_rotvec()
  87. thetas = np.linspace(0, np.pi, n, endpoint=False)
  88. rv = np.pi * np.vstack([np.zeros(n), np.cos(thetas), np.sin(thetas)]).T
  89. g2 = np.roll(rv, axis, axis=1)
  90. return cls.from_rotvec(np.concatenate((g1, g2)))
  91. def cyclic(cls, n, axis=2):
  92. thetas = np.linspace(0, 2 * np.pi, n, endpoint=False)
  93. rv = np.vstack([thetas, np.zeros(n), np.zeros(n)]).T
  94. return cls.from_rotvec(np.roll(rv, axis, axis=1))
  95. def create_group(cls, group, axis='Z'):
  96. if not isinstance(group, str):
  97. raise ValueError("`group` argument must be a string")
  98. permitted_axes = ['x', 'y', 'z', 'X', 'Y', 'Z']
  99. if axis not in permitted_axes:
  100. raise ValueError("`axis` must be one of " + ", ".join(permitted_axes))
  101. if group in ['I', 'O', 'T']:
  102. symbol = group
  103. order = 1
  104. elif group[:1] in ['C', 'D'] and group[1:].isdigit():
  105. symbol = group[:1]
  106. order = int(group[1:])
  107. else:
  108. raise ValueError("`group` must be one of 'I', 'O', 'T', 'Dn', 'Cn'")
  109. if order < 1:
  110. raise ValueError("Group order must be positive")
  111. axis = 'xyz'.index(axis.lower())
  112. if symbol == 'I':
  113. return icosahedral(cls)
  114. elif symbol == 'O':
  115. return octahedral(cls)
  116. elif symbol == 'T':
  117. return tetrahedral(cls)
  118. elif symbol == 'D':
  119. return dicyclic(cls, order, axis=axis)
  120. elif symbol == 'C':
  121. return cyclic(cls, order, axis=axis)
  122. else:
  123. assert False