1 /* pointgroup.c */
2 /* Copyright (C) 2008 Atsushi Togo */
3 
4 #include <string.h>
5 #include <stdio.h>
6 #include "lattice.h"
7 #include "pointgroup.h"
8 #include "symmetry.h"
9 #include "mathfunc.h"
10 
11 #include "debug.h"
12 
13 #define NUM_ROT_AXES 73
14 
15 typedef struct {
16   int table[10];
17   char symbol[6];
18   Holohedry holohedry;
19   Laue laue;
20 } PointgroupType;
21 
22 static PointgroupType pointgroup_data[32] = {
23   {
24     {0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
25     "1    ",
26     TRICLI,
27     LAUE1,
28   },
29   {
30     {0, 0, 0, 0, 1, 1, 0, 0, 0, 0},
31     "-1   ",
32     TRICLI,
33     LAUE1,
34   },
35   {
36     {0, 0, 0, 0, 0, 1, 1, 0, 0, 0},
37     "2    ",
38     MONOCLI,
39     LAUE2M,
40   },
41   {
42     {0, 0, 0, 1, 0, 1, 0, 0, 0, 0},
43     "m    ",
44     MONOCLI,
45     LAUE2M,
46   },
47   {
48     {0, 0, 0, 1, 1, 1, 1, 0, 0, 0},
49     "2/m  ",
50     MONOCLI,
51     LAUE2M,
52   },
53   {
54     {0, 0, 0, 0, 0, 1, 3, 0, 0, 0},
55     "222  ",
56     ORTHO,
57     LAUEMMM,
58   },
59   {
60     {0, 0, 0, 2, 0, 1, 1, 0, 0, 0},
61     "mm2  ",
62     ORTHO,
63     LAUEMMM,
64   },
65   {
66     {0, 0, 0, 3, 1, 1, 3, 0, 0, 0},
67     "mmm  ",
68     ORTHO,
69     LAUEMMM,
70   },
71   {
72     {0, 0, 0, 0, 0, 1, 1, 0, 2, 0},
73     "4    ",
74     TETRA,
75     LAUE4M,
76   },
77   {
78     {0, 2, 0, 0, 0, 1, 1, 0, 0, 0},
79     "-4   ",
80     TETRA,
81     LAUE4M,
82   },
83   {
84     {0, 2, 0, 1, 1, 1, 1, 0, 2, 0},
85     "4/m  ",
86     TETRA,
87     LAUE4M,
88   },
89   {
90     {0, 0, 0, 0, 0, 1, 5, 0, 2, 0},
91     "422  ",
92     TETRA,
93     LAUE4MMM,
94   },
95   {
96     {0, 0, 0, 4, 0, 1, 1, 0, 2, 0},
97     "4mm  ",
98     TETRA,
99     LAUE4MMM,
100   },
101   {
102     {0, 2, 0, 2, 0, 1, 3, 0, 0, 0},
103     "-42m ",
104     TETRA,
105     LAUE4MMM,
106   },
107   {
108     {0, 2, 0, 5, 1, 1, 5, 0, 2, 0},
109     "4/mmm",
110     TETRA,
111     LAUE4MMM,
112   },
113   {
114     {0, 0, 0, 0, 0, 1, 0, 2, 0, 0},
115     "3    ",
116     TRIGO,
117     LAUE3,
118   },
119   {
120     {0, 0, 2, 0, 1, 1, 0, 2, 0, 0},
121     "-3   ",
122     TRIGO,
123     LAUE3,
124   },
125   {
126     {0, 0, 0, 0, 0, 1, 3, 2, 0, 0},
127     "32   ",
128     TRIGO,
129     LAUE3M,
130   },
131   {
132     {0, 0, 0, 3, 0, 1, 0, 2, 0, 0},
133     "3m   ",
134     TRIGO,
135     LAUE3M,
136   },
137   {
138     {0, 0, 2, 3, 1, 1, 3, 2, 0, 0},
139     "-3m  ",
140     TRIGO,
141     LAUE3M,
142   },
143   {
144     {0, 0, 0, 0, 0, 1, 1, 2, 0, 2},
145     "6    ",
146     HEXA,
147     LAUE6M,
148   },
149   {
150     {2, 0, 0, 1, 0, 1, 0, 2, 0, 0},
151     "-6   ",
152     HEXA,
153     LAUE6M,
154   },
155   {
156     {2, 0, 2, 1, 1, 1, 1, 2, 0, 2},
157     "6/m  ",
158     HEXA,
159     LAUE6M,
160   },
161   {
162     {0, 0, 0, 0, 0, 1, 7, 2, 0, 2},
163     "622  ",
164     HEXA,
165     LAUE6MMM,
166   },
167   {
168     {0, 0, 0, 6, 0, 1, 1, 2, 0, 2},
169     "6mm  ",
170     HEXA,
171     LAUE6MMM,
172   },
173   {
174     {2, 0, 0, 4, 0, 1, 3, 2, 0, 0},
175     "-62m ",
176     HEXA,
177     LAUE6MMM,
178   },
179   {
180     {2, 0, 2, 7, 1, 1, 7, 2, 0, 2},
181     "6/mmm",
182     HEXA,
183     LAUE6MMM,
184   },
185   {
186     {0, 0, 0, 0, 0, 1, 3, 8, 0, 0},
187     "23   ",
188     CUBIC,
189     LAUEM3,
190   },
191   {
192     {0, 0, 8, 3, 1, 1, 3, 8, 0, 0},
193     "m-3  ",
194     CUBIC,
195     LAUEM3,
196   },
197   {
198     {0, 0, 0, 0, 0, 1, 9, 8, 6, 0},
199     "432  ",
200     CUBIC,
201     LAUEM3M,
202   },
203   {
204     {0, 6, 0, 6, 0, 1, 3, 8, 0, 0},
205     "-43m ",
206     CUBIC,
207     LAUEM3M,
208   },
209   {
210     {0, 6, 8, 9, 1, 1, 9, 8, 6, 0},
211     "m-3m ",
212     CUBIC,
213     LAUEM3M,
214   }
215 };
216 
217 static int identity[3][3] = {
218   { 1, 0, 0},
219   { 0, 1, 0},
220   { 0, 0, 1},
221 };
222 
223 static int inversion[3][3] = {
224   {-1, 0, 0},
225   { 0,-1, 0},
226   { 0, 0,-1},
227 };
228 
229 static int rot_axes[][3] = {
230   { 1, 0, 0},
231   { 0, 1, 0},
232   { 0, 0, 1},
233   { 0, 1, 1},
234   { 1, 0, 1},
235   { 1, 1, 0},
236   { 0, 1,-1},
237   {-1, 0, 1},
238   { 1,-1, 0},
239   { 1, 1, 1}, /* 10 */
240   {-1, 1, 1},
241   { 1,-1, 1},
242   { 1, 1,-1},
243   { 0, 1, 2},
244   { 2, 0, 1},
245   { 1, 2, 0},
246   { 0, 2, 1},
247   { 1, 0, 2},
248   { 2, 1, 0},
249   { 0,-1, 2}, /* 20 */
250   { 2, 0,-1},
251   {-1, 2, 0},
252   { 0,-2, 1},
253   { 1, 0,-2},
254   {-2, 1, 0},
255   { 2, 1, 1},
256   { 1, 2, 1},
257   { 1, 1, 2},
258   { 2,-1,-1},
259   {-1, 2,-1}, /* 30 */
260   {-1,-1, 2},
261   { 2, 1,-1},
262   {-1, 2, 1},
263   { 1,-1, 2},
264   { 2,-1, 1}, /* 35 */
265   { 1, 2,-1},
266   {-1, 1, 2},
267   { 3, 1, 2},
268   { 2, 3, 1},
269   { 1, 2, 3}, /* 40 */
270   { 3, 2, 1},
271   { 1, 3, 2},
272   { 2, 1, 3},
273   { 3,-1, 2},
274   { 2, 3,-1}, /* 45 */
275   {-1, 2, 3},
276   { 3,-2, 1},
277   { 1, 3,-2},
278   {-2, 1, 3},
279   { 3,-1,-2}, /* 50 */
280   {-2, 3,-1},
281   {-1,-2, 3},
282   { 3,-2,-1},
283   {-1, 3,-2},
284   {-2,-1, 3}, /* 55 */
285   { 3, 1,-2},
286   {-2, 3, 1},
287   { 1,-2, 3},
288   { 3, 2,-1},
289   {-1, 3, 2}, /* 60 */
290   { 2,-1, 3},
291   { 1, 1, 3},
292   {-1, 1, 3},
293   { 1,-1, 3},
294   {-1,-1, 3}, /* 65 */
295   { 1, 3, 1},
296   {-1, 3, 1},
297   { 1, 3,-1},
298   {-1, 3,-1},
299   { 3, 1, 1}, /* 70 */
300   { 3, 1,-1},
301   { 3,-1, 1},
302   { 3,-1,-1},
303 };
304 
305 static int get_pointgroup_number( const Symmetry * symmetry );
306 static int get_pointgroup_class_table( int table[10],
307 				       const Symmetry * symmetry );
308 static int get_rotation_type( SPGCONST int rot[3][3] );
309 static int get_rotation_axis( SPGCONST int rot[3][3] );
310 static int get_orthogonal_axis( int ortho_axes[],
311 				SPGCONST int proper_rot[3][3],
312 				const int rot_order );
313 static int laue2m( int axes[3],
314 		   const Symmetry * symmetry );
315 
316 #ifdef DEBUG
317 static int lauemmm( int axes[3],
318 		    const Symmetry * symmetry );
319 static int laue4m( int axes[3],
320 		   const Symmetry * symmetry );
321 static int laue4mmm( int axes[3],
322 		     const Symmetry * symmetry );
323 static int laue3( int axes[3],
324 		  const Symmetry * symmetry );
325 static int laue3m( int axes[3],
326 		   const Symmetry * symmetry );
327 static int lauem3m( int axes[3],
328 		    const Symmetry * symmetry );
329 #endif
330 
331 static int laue_one_axis( int axes[3],
332 			  const Symmetry * symmetry,
333 			  const int rot_order );
334 static int lauennn( int axes[3],
335 		    const Symmetry * symmetry,
336 		    const int rot_order );
337 static int get_axes( int axes[3],
338 		     const Laue laue,
339 		     const Symmetry * symmetry );
340 static void get_proper_rotation( int prop_rot[3][3],
341 				 SPGCONST int rot[3][3] );
342 static void get_transform_matrix( int mat[3][3],
343 				  const int axes[3] );
344 static int is_exist_axis( const int axis_vec[3], const int axis_index );
345 
346 
ptg_get_pointgroup_number(const Symmetry * symmetry)347 int ptg_get_pointgroup_number( const Symmetry * symmetry )
348 {
349   return get_pointgroup_number( symmetry );
350 }
351 
ptg_get_pointgroup(const int pointgroup_number)352 Pointgroup ptg_get_pointgroup( const int pointgroup_number )
353 {
354   Pointgroup pointgroup;
355   PointgroupType pointgroup_type;
356 
357   pointgroup_type = pointgroup_data[ pointgroup_number ];
358   strcpy( pointgroup.symbol, pointgroup_type.symbol );
359   pointgroup.holohedry = pointgroup_type.holohedry;
360   pointgroup.laue = pointgroup_type.laue;
361 
362   debug_print("Point group: %s\n", pointgroup_type.symbol );
363 
364   return pointgroup;
365 }
366 
367 /* pointgroup is modified. */
ptg_get_transformation_matrix(Pointgroup * pointgroup,const Symmetry * symmetry)368 void ptg_get_transformation_matrix( Pointgroup * pointgroup,
369 				    const Symmetry * symmetry )
370 {
371   int axes[3];
372   int transform_mat[3][3];
373 
374   get_axes( axes, pointgroup->laue, symmetry );
375   get_transform_matrix( transform_mat, axes );
376   mat_copy_matrix_i3( pointgroup->transform_mat, transform_mat );
377 }
378 
get_pointgroup_number(const Symmetry * symmetry)379 static int get_pointgroup_number( const Symmetry * symmetry )
380 {
381   int i, j, pg_num, counter;
382   int table[10];
383   PointgroupType pointgroup_type;
384 
385   /* Get list of rotation part of symmetry operations */
386   if ( ! get_pointgroup_class_table( table, symmetry ) ) {
387     pg_num = -1;
388     goto end;
389   }
390 
391   pg_num = -1;
392   for ( i = 0; i < 32; i++ ) {
393     counter = 0;
394     pointgroup_type = pointgroup_data[ i ];
395     for ( j = 0; j < 10; j++ ) {
396       if ( pointgroup_type.table[j] == table[j] ) { counter++; }
397     }
398     if ( counter == 10 ) {
399       pg_num = i;
400       break;
401     }
402   }
403 
404  end:
405   return pg_num;
406 }
407 
get_pointgroup_class_table(int table[10],const Symmetry * symmetry)408 static int get_pointgroup_class_table( int table[10],
409 				       const Symmetry * symmetry )
410 {
411   /* Look-up table */
412   /* Operation   -6 -4 -3 -2 -1  1  2  3  4  6 */
413   /* Trace     -  2 -1  0  1 -3  3 -1  0  1  2 */
414   /* Determinant -1 -1 -1 -1 -1  1  1  1  1  1 */
415 
416   /* table[0] = -6 axis */
417   /* table[1] = -4 axis */
418   /* table[2] = -3 axis */
419   /* table[3] = -2 axis */
420   /* table[4] = -1 axis */
421   /* table[5] =  1 axis */
422   /* table[6] =  2 axis */
423   /* table[7] =  3 axis */
424   /* table[8] =  4 axis */
425   /* table[9] =  6 axis */
426 
427   int i, rot_type;
428 
429   for ( i = 0; i < 10; i++ ) { table[i] = 0; }
430   for ( i = 0; i < symmetry->size; i++ ) {
431     rot_type = get_rotation_type( symmetry->rot[i] );
432     if ( rot_type == -1 ) {
433       goto err;
434     } else {
435       table[rot_type]++;
436     }
437   }
438 
439   return 1;
440 
441  err:
442   warning_print("spglib: No point group symbol found ");
443   warning_print("(line %d, %s).\n", __LINE__, __FILE__);
444   return 0;
445 }
446 
get_rotation_type(SPGCONST int rot[3][3])447 static int get_rotation_type( SPGCONST int rot[3][3] )
448 {
449   int rot_type;
450 
451   if ( mat_get_determinant_i3( rot ) == -1 ) {
452     switch ( mat_get_trace_i3( rot ) ) {
453     case -2: /* -6 */
454       rot_type = 0;
455       break;
456     case -1: /* -4 */
457       rot_type = 1;
458       break;
459     case 0:  /* -3 */
460       rot_type = 2;
461       break;
462     case 1:  /* -2 */
463       rot_type = 3;
464       break;
465     case -3: /* -1 */
466       rot_type = 4;
467       break;
468     default:
469       rot_type = -1;
470       break;
471     }
472   } else {
473     switch ( mat_get_trace_i3( rot ) ) {
474     case 3:  /* 1 */
475       rot_type = 5;
476       break;
477     case -1: /* 2 */
478       rot_type = 6;
479       break;
480     case 0:  /* 3 */
481       rot_type = 7;
482       break;
483     case 1:  /* 4 */
484       rot_type = 8;
485       break;
486     case 2:  /* 6 */
487       rot_type = 9;
488       break;
489     default:
490       rot_type = -1;
491       break;
492     }
493   }
494 
495   return rot_type;
496 }
497 
get_axes(int axes[3],const Laue laue,const Symmetry * symmetry)498 static int get_axes( int axes[3],
499 		     const Laue laue,
500 		     const Symmetry * symmetry )
501 {
502   switch (laue) {
503   case LAUE1:
504     axes[0] = 0;
505     axes[1] = 1;
506     axes[2] = 2;
507     break;
508   case LAUE2M:
509     laue2m( axes, symmetry );
510     break;
511   case LAUEMMM:
512     lauennn( axes, symmetry, 2 );
513     break;
514   case LAUE4M:
515     laue_one_axis( axes, symmetry, 4 );
516     break;
517   case LAUE4MMM:
518     laue_one_axis( axes, symmetry, 4 );
519     break;
520   case LAUE3:
521     laue_one_axis( axes, symmetry, 3 );
522     break;
523   case LAUE3M:
524     laue_one_axis( axes, symmetry, 3 );
525     break;
526   case LAUE6M:
527     laue_one_axis( axes, symmetry, 3 );
528     break;
529   case LAUE6MMM:
530     laue_one_axis( axes, symmetry, 3 );
531     break;
532   case LAUEM3:
533     lauennn( axes, symmetry, 2 );
534     break;
535   case LAUEM3M:
536     lauennn( axes, symmetry, 4 );
537     break;
538   default:
539     break;
540   }
541 
542   return 1;
543 }
544 
laue2m(int axes[3],const Symmetry * symmetry)545 static int laue2m( int axes[3],
546 		   const Symmetry * symmetry )
547 {
548   int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
549   int prop_rot[3][3], t_mat[3][3];
550   int ortho_axes[NUM_ROT_AXES];
551 
552   for ( i = 0; i < symmetry->size; i++ ) {
553     get_proper_rotation( prop_rot, symmetry->rot[i] );
554 
555     /* Search two-fold rotation */
556     if ( ! ( mat_get_trace_i3( prop_rot ) == -1 ) ) {
557       continue;
558     }
559 
560     /* The first axis */
561     axes[1] = get_rotation_axis( prop_rot );
562     break;
563   }
564 
565   /* The second axis */
566   num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, 2 );
567   if ( ! num_ortho_axis ) { goto err; }
568 
569   min_norm = 8;
570   is_found = 0;
571   for ( i = 0; i < num_ortho_axis; i++ ) {
572     norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
573     if ( norm < min_norm ) {
574       min_norm = norm;
575       axes[0] = ortho_axes[i];
576       is_found = 1;
577     }
578   }
579   if ( ! is_found ) { goto err; }
580 
581   /* The third axis */
582   min_norm = 8;
583   is_found = 0;
584   for ( i = 0; i < num_ortho_axis; i++ ) {
585     norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
586     if ( norm < min_norm && ( ! ( ortho_axes[i] == axes[0] ) ) ) {
587       min_norm = norm;
588       axes[2] = ortho_axes[i];
589       is_found = 1;
590     }
591   }
592   if ( ! is_found ) { goto err; }
593 
594   get_transform_matrix( t_mat, axes );
595   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
596     tmpval = axes[0];
597     axes[0] = axes[2];
598     axes[2] = tmpval;
599   }
600 
601   return 1;
602 
603  err:
604   return 0;
605 }
606 
607 #ifdef DEBUG
lauemmm(int axes[3],const Symmetry * symmetry)608 static int lauemmm( int axes[3],
609 		    const Symmetry * symmetry )
610 {
611   int i, count, axis, tmpval;
612   int prop_rot[3][3], t_mat[3][3];
613 
614   for ( i = 0; i < 3; i++ ) { axes[i] = -1; }
615   count = 0;
616   for ( i = 0; i < symmetry->size; i++ ) {
617     get_proper_rotation( prop_rot, symmetry->rot[i] );
618 
619     /* Search two-fold rotation */
620     if ( ! ( mat_get_trace_i3( prop_rot ) == -1 ) ) {
621       continue;
622     }
623 
624     axis = get_rotation_axis( prop_rot );
625     if ( ! ( ( axis == axes[0] ) ||
626 	     ( axis == axes[1] ) ||
627 	     ( axis == axes[2] ) ) ) {
628       axes[count] = axis;
629       count++;
630     }
631   }
632 
633   get_transform_matrix( t_mat, axes );
634   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
635     tmpval = axes[0];
636     axes[0] = axes[1];
637     axes[1] = tmpval;
638   }
639 
640   return 1;
641 }
642 
laue4m(int axes[3],const Symmetry * symmetry)643 static int laue4m( int axes[3],
644 		   const Symmetry * symmetry )
645 {
646   int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
647   int axis_vec[3];
648   int prop_rot[3][3], t_mat[3][3];
649   int ortho_axes[NUM_ROT_AXES];
650 
651   for ( i = 0; i < symmetry->size; i++ ) {
652     get_proper_rotation( prop_rot, symmetry->rot[i] );
653 
654     /* Search foud-fold rotation */
655     if (  mat_get_trace_i3( prop_rot ) == 1 ) {
656       /* The first axis */
657       axes[2] = get_rotation_axis( prop_rot );
658       break;
659     }
660   }
661 
662   /* The second axis */
663   num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, 4 );
664   if ( ! num_ortho_axis ) { goto err; }
665 
666   min_norm = 8;
667   is_found = 0;
668   for ( i = 0; i < num_ortho_axis; i++ ) {
669     norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
670     if ( norm < min_norm ) {
671       min_norm = norm;
672       axes[0] = ortho_axes[i];
673       is_found = 1;
674     }
675   }
676   if ( ! is_found ) { goto err; }
677 
678   /* The third axis */
679   mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
680   is_found = 0;
681   for ( i = 0; i < NUM_ROT_AXES; i++ ) {
682     if ( is_exist_axis( axis_vec, i ) ) {
683       is_found = 1;
684       axes[1] = i;
685       break;
686     }
687   }
688   if ( ! is_found ) { goto err; }
689 
690   get_transform_matrix( t_mat, axes );
691   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
692     tmpval = axes[0];
693     axes[0] = axes[1];
694     axes[1] = tmpval;
695   }
696 
697   return 1;
698 
699  err:
700   return 0;
701 }
702 
laue4mmm(int axes[3],const Symmetry * symmetry)703 static int laue4mmm( int axes[3],
704 		     const Symmetry * symmetry )
705 {
706   int i, is_found, tmpval, axis;
707   int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
708   int axis_vec[3];
709 
710   for ( i = 0; i < symmetry->size; i++ ) {
711     get_proper_rotation( prop_rot, symmetry->rot[i] );
712 
713     /* Search foud-fold rotation */
714     if ( mat_get_trace_i3( prop_rot ) == 1 ) {
715       /* The first axis */
716       axes[2] = get_rotation_axis( prop_rot );
717       break;
718     }
719   }
720 
721   is_found = 0;
722   for ( i = 0; i < symmetry->size; i++ ) {
723     get_proper_rotation( prop_rot2, symmetry->rot[i] );
724 
725     /* Search two-fold rotation */
726     if ( ! ( mat_get_trace_i3( prop_rot2 ) == -1 ) ) {
727       continue;
728     }
729 
730     /* The second axis */
731     axis = get_rotation_axis( prop_rot2 );
732     if ( ! ( axis == axes[2] ) ) {
733       axes[0] = axis;
734       is_found = 1;
735       break;
736     }
737   }
738   if ( ! is_found ) { goto err; }
739 
740 
741   /* The third axis */
742   mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
743   is_found = 0;
744   for ( i = 0; i < NUM_ROT_AXES; i++ ) {
745     if ( is_exist_axis( axis_vec, i ) ) {
746       is_found = 1;
747       axes[1] = i;
748       break;
749     }
750   }
751   if ( ! is_found ) { goto err; }
752 
753   get_transform_matrix( t_mat, axes );
754   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
755     tmpval = axes[0];
756     axes[0] = axes[1];
757     axes[1] = tmpval;
758   }
759 
760   return 1;
761 
762  err:
763   return 0;
764 }
765 
laue3(int axes[3],const Symmetry * symmetry)766 static int laue3( int axes[3],
767 		  const Symmetry * symmetry )
768 {
769   int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
770   int prop_rot[3][3], t_mat[3][3];
771   int axis_vec[3];
772   int ortho_axes[NUM_ROT_AXES];
773 
774   for ( i = 0; i < symmetry->size; i++ ) {
775     get_proper_rotation( prop_rot, symmetry->rot[i] );
776 
777     /* Search thee-fold rotation */
778     if ( mat_get_trace_i3( prop_rot ) == 0 ) {
779       /* The first axis */
780       axes[2] = get_rotation_axis( prop_rot );
781       break;
782     }
783   }
784 
785   /* The second axis */
786   num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, 3 );
787   if ( ! num_ortho_axis ) { goto err; }
788   min_norm = 8;
789   is_found = 0;
790   for ( i = 0; i < num_ortho_axis; i++ ) {
791     norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
792     if ( norm < min_norm ) {
793       min_norm = norm;
794       axes[0] = ortho_axes[i];
795       is_found = 1;
796     }
797   }
798   if ( ! is_found ) { goto err; }
799 
800   /* The third axis */
801   mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
802   is_found = 0;
803   for ( i = 0; i < NUM_ROT_AXES; i++ ) {
804     is_found = is_exist_axis( axis_vec, i );
805     if ( is_found == 1 ) {
806       axes[1] = i;
807       break;
808     }
809     if ( is_found == -1 ) {
810       axes[1] = i + NUM_ROT_AXES;
811       break;
812     }
813   }
814   if ( ! is_found ) { goto err; }
815 
816   get_transform_matrix( t_mat, axes );
817   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
818     tmpval = axes[0];
819     axes[0] = axes[1];
820     axes[1] = tmpval;
821   }
822 
823   return 1;
824 
825  err:
826   return 0;
827 }
828 
laue3m(int axes[3],const Symmetry * symmetry)829 static int laue3m( int axes[3],
830 		   const Symmetry * symmetry )
831 {
832   int i, is_found, tmpval, axis;
833   int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
834   int axis_vec[3];
835 
836   for ( i = 0; i < symmetry->size; i++ ) {
837     get_proper_rotation( prop_rot, symmetry->rot[i] );
838 
839     /* Search three-fold rotation */
840     if ( mat_get_trace_i3( prop_rot ) == 0 ) {
841       /* The first axis */
842       axes[2] = get_rotation_axis( prop_rot );
843       debug_print("laue3m prop_rot\n");
844       debug_print_matrix_i3( prop_rot );
845       break;
846     }
847   }
848 
849   is_found = 0;
850   for ( i = 0; i < symmetry->size; i++ ) {
851     get_proper_rotation( prop_rot2, symmetry->rot[i] );
852 
853     /* Search two-fold rotation */
854     if ( ! ( mat_get_trace_i3( prop_rot2 ) == -1 ) ) {
855       continue;
856     }
857 
858     /* The second axis */
859     axis = get_rotation_axis( prop_rot2 );
860     if ( ! ( axis == axes[2] ) ) {
861       axes[0] = axis;
862       is_found = 1;
863       break;
864     }
865   }
866   if ( ! is_found ) { goto err; }
867 
868   /* The third axis */
869   mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
870   is_found = 0;
871   for ( i = 0; i < NUM_ROT_AXES; i++ ) {
872     is_found = is_exist_axis( axis_vec, i );
873     if ( is_found == 1 ) {
874       axes[1] = i;
875       break;
876     }
877     if ( is_found == -1 ) {
878       axes[1] = i + NUM_ROT_AXES;
879       break;
880     }
881   }
882   if ( ! is_found ) { goto err; }
883 
884   get_transform_matrix( t_mat, axes );
885   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
886     tmpval = axes[0];
887     axes[0] = axes[1];
888     axes[1] = tmpval;
889   }
890 
891   return 1;
892 
893  err:
894   return 0;
895 }
896 
lauem3m(int axes[3],const Symmetry * symmetry)897 static int lauem3m( int axes[3],
898 		    const Symmetry * symmetry )
899 {
900   int i, count, axis, tmpval;
901   int prop_rot[3][3], t_mat[3][3];
902 
903   for ( i = 0; i < 3; i++ ) { axes[i] = -1; }
904   count = 0;
905   for ( i = 0; i < symmetry->size; i++ ) {
906     get_proper_rotation( prop_rot, symmetry->rot[i] );
907 
908     /* Search four-fold rotation */
909     if ( ! ( mat_get_trace_i3( prop_rot ) == 1 ) ) {
910       continue;
911     }
912 
913     axis = get_rotation_axis( prop_rot );
914     if ( ! ( ( axis == axes[0] ) ||
915 	     ( axis == axes[1] ) ||
916 	     ( axis == axes[2] ) ) ) {
917       axes[count] = axis;
918       count++;
919     }
920   }
921 
922   get_transform_matrix( t_mat, axes );
923   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
924     tmpval = axes[0];
925     axes[0] = axes[1];
926     axes[1] = tmpval;
927   }
928 
929   return 1;
930 }
931 #endif
932 
laue_one_axis(int axes[3],const Symmetry * symmetry,const int rot_order)933 static int laue_one_axis( int axes[3],
934 			  const Symmetry * symmetry,
935 			  const int rot_order )
936 {
937   int i, j, num_ortho_axis, det, min_det, is_found, tmpval;
938   int axis_vec[3], tmp_axes[3];
939   int prop_rot[3][3], t_mat[3][3];
940   int ortho_axes[NUM_ROT_AXES];
941 
942   for ( i = 0; i < symmetry->size; i++ ) {
943     get_proper_rotation( prop_rot, symmetry->rot[i] );
944 
945     /* Search foud-fold rotation */
946     if ( rot_order == 4 ) {
947       if (  mat_get_trace_i3( prop_rot ) == 1 ) {
948 	/* The first axis */
949 	axes[2] = get_rotation_axis( prop_rot );
950 	break;
951       }
952     }
953 
954     /* Search three-fold rotation */
955     if ( rot_order == 3 ) {
956       if (  mat_get_trace_i3( prop_rot ) == 0 ) {
957 	/* The first axis */
958 	axes[2] = get_rotation_axis( prop_rot );
959 	break;
960       }
961     }
962   }
963 
964   /* Candidates of the second axis */
965   num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, rot_order );
966   if ( ! num_ortho_axis ) { goto err; }
967 
968   tmp_axes[2] = axes[2];
969   min_det = 4;
970   is_found = 0;
971   for ( i = 0; i < num_ortho_axis; i++ ) {
972     tmp_axes[0] = ortho_axes[i];
973     mat_multiply_matrix_vector_i3( axis_vec, prop_rot,
974 				   rot_axes[tmp_axes[0]] );
975     for ( j = 0; j < num_ortho_axis; j++ ) {
976       is_found = is_exist_axis( axis_vec, ortho_axes[j] );
977       if ( is_found == 1 ) {
978 	tmp_axes[1] = ortho_axes[j];
979 	break;
980       }
981       if ( is_found == -1 ) {
982 	tmp_axes[1] = ortho_axes[j] + NUM_ROT_AXES;
983 	break;
984       }
985     }
986 
987     get_transform_matrix( t_mat, tmp_axes );
988     det = mat_get_determinant_i3( t_mat );
989     if ( det < 0 ) { det = -det; }
990     if ( det < min_det ) {
991       min_det = det;
992       axes[0] = tmp_axes[0];
993       axes[1] = tmp_axes[1];
994     }
995   }
996   if ( ! is_found ) { goto err; }
997 
998   get_transform_matrix( t_mat, axes );
999   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
1000     tmpval = axes[0];
1001     axes[0] = axes[1];
1002     axes[1] = tmpval;
1003   }
1004 
1005   return 1;
1006 
1007  err:
1008   return 0;
1009 }
1010 
lauennn(int axes[3],const Symmetry * symmetry,const int rot_order)1011 static int lauennn( int axes[3],
1012 		    const Symmetry * symmetry,
1013 		    const int rot_order )
1014 {
1015   int i, count, axis, tmpval;
1016   int prop_rot[3][3], t_mat[3][3];
1017 
1018   for ( i = 0; i < 3; i++ ) { axes[i] = -1; }
1019   count = 0;
1020   for ( i = 0; i < symmetry->size; i++ ) {
1021     get_proper_rotation( prop_rot, symmetry->rot[i] );
1022 
1023     /* Search two- or four-fold rotation */
1024     if ( ( mat_get_trace_i3( prop_rot ) == -1 && rot_order == 2) ||
1025 	 ( mat_get_trace_i3( prop_rot ) == 1 && rot_order == 4 ) ) {
1026       axis = get_rotation_axis( prop_rot );
1027       if ( ! ( ( axis == axes[0] ) ||
1028 	       ( axis == axes[1] ) ||
1029 	       ( axis == axes[2] ) ) ) {
1030 	axes[count] = axis;
1031 	count++;
1032       }
1033     }
1034   }
1035 
1036   get_transform_matrix( t_mat, axes );
1037   if ( mat_get_determinant_i3( t_mat ) < 0 ) {
1038     tmpval = axes[0];
1039     axes[0] = axes[1];
1040     axes[1] = tmpval;
1041   }
1042 
1043   return 1;
1044 }
1045 
get_rotation_axis(SPGCONST int proper_rot[3][3])1046 static int get_rotation_axis( SPGCONST int proper_rot[3][3] )
1047 {
1048   int i, axis = -1;
1049   int vec[3];
1050 
1051   /* No specific axis for I and -I */
1052   if ( mat_check_identity_matrix_i3( proper_rot, identity ) ) {
1053     goto end;
1054   }
1055 
1056   /* Look for eigenvector = rotation axis */
1057   for ( i = 0; i < NUM_ROT_AXES; i++ ) {
1058     mat_multiply_matrix_vector_i3( vec, proper_rot, rot_axes[i] );
1059     if ( vec[0] == rot_axes[i][0] &&
1060 	 vec[1] == rot_axes[i][1] &&
1061 	 vec[2] == rot_axes[i][2] ) {
1062       axis = i;
1063       break;
1064     }
1065   }
1066 
1067  end:
1068 #ifdef DEBUG
1069   if ( axis == -1 ) {
1070     printf("rotation axis cound not found.\n");
1071   }
1072 #endif
1073 
1074   return axis;
1075 }
1076 
get_orthogonal_axis(int ortho_axes[],SPGCONST int proper_rot[3][3],const int rot_order)1077 static int get_orthogonal_axis( int ortho_axes[],
1078 				SPGCONST int proper_rot[3][3],
1079 				const int rot_order )
1080 {
1081   int i, num_ortho_axis;
1082   int vec[3];
1083   int sum_rot[3][3], rot[3][3];
1084 
1085   num_ortho_axis = 0;
1086 
1087   mat_copy_matrix_i3( sum_rot, identity );
1088   mat_copy_matrix_i3( rot, identity );
1089   for ( i = 0; i < rot_order-1; i++ ) {
1090     mat_multiply_matrix_i3( rot, proper_rot, rot );
1091     mat_add_matrix_i3( sum_rot, rot, sum_rot );
1092   }
1093 
1094   for ( i = 0; i < NUM_ROT_AXES; i++ ) {
1095     mat_multiply_matrix_vector_i3( vec, sum_rot, rot_axes[i] );
1096     if ( vec[0] == 0 &&
1097 	 vec[1] == 0 &&
1098 	 vec[2] == 0 ) {
1099       ortho_axes[ num_ortho_axis ] = i;
1100       num_ortho_axis++;
1101     }
1102   }
1103 
1104   return num_ortho_axis;
1105 }
1106 
get_proper_rotation(int prop_rot[3][3],SPGCONST int rot[3][3])1107 static void get_proper_rotation( int prop_rot[3][3],
1108 				 SPGCONST int rot[3][3] )
1109 {
1110   if ( mat_get_determinant_i3( rot ) == -1 ) {
1111     mat_multiply_matrix_i3( prop_rot, inversion, rot );
1112   } else {
1113     mat_copy_matrix_i3( prop_rot, rot );
1114   }
1115 }
1116 
get_transform_matrix(int mat[3][3],const int axes[3])1117 static void get_transform_matrix( int mat[3][3],
1118 				  const int axes[3] )
1119 {
1120   int i, j, s[3];
1121 
1122   for ( i = 0; i < 3; i++ ) {
1123     if ( axes[i] < NUM_ROT_AXES ) {
1124       s[i] = 1;
1125     } else {
1126       s[i] = -1;
1127     }
1128   }
1129   for ( i = 0; i < 3; i++ ) {
1130     for ( j = 0; j < 3; j++ ) {
1131       mat[i][j] = s[j] * rot_axes[axes[j]%NUM_ROT_AXES][i];
1132     }
1133   }
1134 }
1135 
is_exist_axis(const int axis_vec[3],const int axis_index)1136 static int is_exist_axis( const int axis_vec[3], const int axis_index )
1137 {
1138   if ( ( axis_vec[0] == rot_axes[axis_index][0] ) &&
1139        ( axis_vec[1] == rot_axes[axis_index][1] ) &&
1140        ( axis_vec[2] == rot_axes[axis_index][2] ) ) { return 1; }
1141   if ( ( axis_vec[0] == -rot_axes[axis_index][0] ) &&
1142        ( axis_vec[1] == -rot_axes[axis_index][1] ) &&
1143        ( axis_vec[2] == -rot_axes[axis_index][2] ) ) { return -1; }
1144   return 0;
1145 }
1146 
1147