1 //
2 // symmetry.c
3 // libmsym
4 //
5 // Created by Marcus Johansson on 12/04/14.
6 // Copyright (c) 2014 Marcus Johansson.
7 //
8 // Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
9 //
10
11 #include "symmetry.h"
12 #include "msym.h"
13 #include "linalg.h"
14 #include "symop.h"
15 #include "geometry.h"
16 #include "equivalence_set.h"
17 #include <stdlib.h>
18 #include <math.h>
19 #include <float.h>
20 #include <time.h>
21
22
23 //Special case of less than (basically with some margin)
24 #define LT(A,B,T) ((B) - (A) > (T))
25 #ifndef M_PI
26 #define M_PI 3.14159265358979323846264338327950288419716939937510582
27 #endif
28
29 int divisors(int n, int* div);
30
31 msym_error_t findEquivalenceSetSymmetryOperations(msym_equivalence_set_t *es, msym_thresholds_t *t, int *lsops, msym_symmetry_operation_t **sops);
32 msym_error_t findSymmetryLinear(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
33 msym_error_t findSymmetryPlanarRegular(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *t, int *rsopsl, msym_symmetry_operation_t **rsops);
34 msym_error_t findSymmetryPlanarIrregular(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
35 msym_error_t findSymmetryPolyhedralProlate(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
36 msym_error_t findSymmetryPolyhedralOblate(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
37 msym_error_t findSymmetrySymmetricPolyhedron(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int prim, int *rsopsl, msym_symmetry_operation_t **rsops);
38 msym_error_t findSymmetryAsymmetricPolyhedron(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
39 msym_error_t findSymmetrySpherical(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
40 msym_error_t findSymmetryCubic(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops);
41 msym_error_t findSymmetryUnknown(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *t, int *rsopsl, msym_symmetry_operation_t **rsops);
42 msym_error_t reduceSymmetry(int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, int *isopsl, msym_symmetry_operation_t **isops);
43 msym_error_t filterSymmetryOperations(int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, int *isopsl, msym_symmetry_operation_t **isops);
44
findSymmetryOperations(int esl,msym_equivalence_set_t es[esl],msym_thresholds_t * t,int * lsops,msym_symmetry_operation_t ** sops)45 msym_error_t findSymmetryOperations(int esl, msym_equivalence_set_t es[esl], msym_thresholds_t *t, int *lsops, msym_symmetry_operation_t **sops){
46 msym_error_t ret = MSYM_SUCCESS;
47 msym_symmetry_operation_t *rsops = NULL;
48 int lrsops = 0;
49 for(int i = 0; i < esl;i++){
50 int llsops = lrsops;
51 if(MSYM_SUCCESS != (ret = findEquivalenceSetSymmetryOperations(&es[i], t, &lrsops, &rsops))) goto err;
52
53 if(llsops > 0 && lrsops == 0) {
54 free(rsops);
55 rsops = NULL;
56 break;
57 }
58 }
59
60 for(int i = 0;i < lrsops;i++){
61 vnorm(rsops[i].v);
62 }
63
64 *lsops = lrsops;
65 *sops = rsops;
66 return ret;
67 err:
68 free(rsops);
69 *sops = NULL;
70 *lsops = 0;
71 return ret;
72 }
73
findEquivalenceSetSymmetryOperations(msym_equivalence_set_t * es,msym_thresholds_t * t,int * lsops,msym_symmetry_operation_t ** sops)74 msym_error_t findEquivalenceSetSymmetryOperations(msym_equivalence_set_t *es, msym_thresholds_t *t, int *lsops, msym_symmetry_operation_t **sops){
75 //function pointer syntax is a little ambiguous, this is technically less correct, but nicer to read
76
77 const struct _fmap {
78 msym_geometry_t g;
79 msym_error_t (*f)(msym_equivalence_set_t*,double[3],double[3][3],msym_thresholds_t*,int*,msym_symmetry_operation_t**);
80
81 } fmap[8] = {
82 [0] = {MSYM_GEOMETRY_SPHERICAL, findSymmetrySpherical},
83 [1] = {MSYM_GEOMETRY_LINEAR, findSymmetryLinear},
84 [2] = {MSYM_GEOMETRY_PLANAR_REGULAR, findSymmetryPlanarRegular},
85 [3] = {MSYM_GEOMETRY_PLANAR_IRREGULAR, findSymmetryPlanarIrregular},
86 [4] = {MSYM_GEOMETRY_POLYHEDRAL_PROLATE, findSymmetryPolyhedralProlate},
87 [5] = {MSYM_GEOMETRY_POLYHEDRAL_OBLATE, findSymmetryPolyhedralOblate},
88 [6] = {MSYM_GEOMETRY_ASSYMETRIC, findSymmetryAsymmetricPolyhedron},
89 [7] = {MSYM_GEOMETRY_UNKNOWN, findSymmetryUnknown}
90 };
91
92 msym_error_t ret = MSYM_SUCCESS;
93 msym_symmetry_operation_t *fsops = NULL;
94 int lfsops = 0;
95 double cm[3];
96 msym_geometry_t g;
97 double eigvec[3][3];
98 double eigval[3];
99
100
101 if(MSYM_SUCCESS != (ret = findCenterOfMass(es->length, es->elements, cm))) goto err;
102 if(MSYM_SUCCESS != (ret = findGeometry(es->length, es->elements, cm, t, &g, eigval, eigvec))) goto err;
103
104 int fi, fil = sizeof(fmap)/sizeof(fmap[0]);
105 for(fi = 0; fi < fil;fi++){
106 if(fmap[fi].g == g) {
107 if(MSYM_SUCCESS != (ret = fmap[fi].f(es,cm,eigvec,t,&lfsops,&fsops))) goto err;
108 break;
109 }
110 }
111 if(fi == fil){
112 msymSetErrorDetails("Unknown geometry of equivalence set");
113 ret = MSYM_SYMMETRY_ERROR;
114 goto err;
115 }
116
117 if(*sops == NULL){
118 *sops = fsops;
119 *lsops = lfsops;
120 } else if (lfsops != 0) {
121 if(MSYM_SUCCESS != (ret = reduceSymmetry(lfsops, fsops, t, lsops, sops))) goto err;
122 free(fsops);
123 } else if (fsops == 0 && es->length > 1) {
124 msymSetErrorDetails("No symmetry operations found in equivalence set with %d elements",es->length);
125 ret = MSYM_SYMMETRY_ERROR;
126 goto err;
127 } else {
128 free(fsops);
129 }
130
131 return ret;
132
133 err:
134 free(fsops);
135 return ret;
136
137 }
138
findSymmetryUnknown(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)139 msym_error_t findSymmetryUnknown(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
140 return MSYM_INVALID_GEOMETRY;
141 }
142
findSymmetryLinear(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)143 msym_error_t findSymmetryLinear(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
144 msym_error_t ret = MSYM_SUCCESS;
145 int prim = 0, sopsl = 0;
146 msym_symmetry_operation_t *sops = NULL;
147 if(es->length != 2){
148 msymSetErrorDetails("Expected two elements in linear equivalence set, got %d",es->length);
149 ret = MSYM_SYMMETRY_ERROR;
150 goto err;
151 }
152 if(vzero(cm,thresholds->zero)){
153 double t[3], v[3];
154 vnorm2(es->elements[0]->v,v);
155 vnorm2(es->elements[1]->v,t);
156 vadd(v,t,t);
157 vscale(0.5, t, t);
158 vsub(v,t,v);
159 vnorm(v);
160 sopsl = 3;
161 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
162 vcopy(v,sops[0].v);
163 vcopy(v,sops[1].v);
164 sops[0].type = PROPER_ROTATION;
165 sops[0].order = 0;
166 sops[0].power = 1;
167 sops[1].type = REFLECTION;
168 sops[1].order = 1;
169 sops[1].power = 1;
170 sops[2].type = INVERSION;
171 sops[2].order = 1;
172 sops[2].power = 1;
173 sops[2].v[0] = sops[2].v[1] = sops[2].v[2] = 0;
174
175 } else {
176 sopsl = 3;
177 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
178 vcopy(cm,sops[0].v);
179 vnorm(sops[0].v);
180 vcopy(ev[prim],sops[1].v);
181 vnorm(sops[1].v);
182 vcrossnorm(sops[0].v,sops[1].v,sops[2].v);
183 sops[0].type = PROPER_ROTATION;
184 sops[0].order = 2;
185 sops[0].power = 1;
186 sops[1].type = REFLECTION;
187 sops[2].type = REFLECTION;
188 }
189
190 *rsops = sops;
191 *rsopsl = sopsl;
192
193 err:
194 return ret;
195
196 }
197
findSymmetryPlanarRegular(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)198 msym_error_t findSymmetryPlanarRegular(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
199 msym_error_t ret = MSYM_SUCCESS;
200
201 int sigma_h = vzero(cm,thresholds->zero), order = es->length, prim = 2, div_len, even, inversion, n = 0, split = 0;
202 int *div, sopsl;
203 double v0[3], v0_proj[3], v_init[3], theta;
204 msym_symmetry_operation_t *sops = NULL;
205
206 vnorm2(es->elements[0]->v,v0);
207 vproj_plane(v0, ev[prim], v0_proj);
208 vnorm(v0_proj);
209
210 vcopy(v0_proj,v_init);
211
212 for(int i = 1; i < es->length;i++){
213 double vi[3], vi_proj[3];
214 vcopy(es->elements[i]->v,vi);
215 vproj_plane(vi, ev[prim], vi_proj);
216 vnorm(vi);
217 vnorm(vi_proj);
218 theta = vangle(v0_proj,vi_proj);
219
220 if(LT(theta,2*M_PI/es->length,asin(thresholds->angle)) && es->length % 2 == 0){
221 order = es->length/2;
222 split = 1;
223 vadd(v0_proj, vi_proj, v_init);
224 vnorm(v_init);
225 break;
226 }
227 }
228
229 div = malloc(order*sizeof(int));
230 div_len = divisors(order,div);
231
232 even = order % 2 == 0;
233 inversion = even && sigma_h;
234
235 sopsl = (div_len + sigma_h + order + order*sigma_h + inversion + sigma_h*(div_len - even));
236 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
237
238 for(; n < div_len; n++){
239 sops[n].type = PROPER_ROTATION;
240 sops[n].order = div[n];
241 sops[n].power = 1;
242 vcopy(ev[prim],sops[n].v);
243 }
244 if(sigma_h) {
245 sops[n].type = REFLECTION;
246 vcopy(ev[prim],sops[n].v);
247 n++;
248 for(int s = 0; s < div_len; s++){
249 if(div[s] > 2){
250 sops[n].type = IMPROPER_ROTATION;
251 sops[n].order = div[s];
252 sops[n].power = 1;
253 vcopy(ev[prim],sops[n].v);
254 n++;
255 }
256 }
257 }
258
259 if(inversion){
260 sops[n].order = 1;
261 sops[n].power = 1;
262 sops[n].type = INVERSION;
263 sops[n].v[0] = sops[n].v[1] = sops[n].v[2] = 0;
264 n++;
265 }
266
267 theta = M_PI/order;
268 for (int i = 0; i < order && n < sopsl; i++) {
269 double r[3];
270 vrotate(i*theta, v_init, ev[prim], r);
271 vnorm(r); //Not really nessecary
272 vcrossnorm(r,ev[prim],sops[n].v);
273 sops[n].type = REFLECTION;
274
275 if(!findSymmetryOperation(&sops[n], sops, n, thresholds)){
276 n++;
277 if(sigma_h){
278 vcopy(r,sops[n].v);
279 sops[n].type = PROPER_ROTATION;
280 sops[n].order = 2;
281 sops[n].power = 1;
282 n++;
283 }
284 }
285 }
286
287 free(div);
288 if(n != sopsl) {
289 msymSetErrorDetails("Unexpected number of generated symmetry operations in planar regular polygon. Got %d expected %d",n,sopsl);
290 ret = MSYM_SYMMETRY_ERROR;
291 goto err;
292 }
293
294 *rsops = sops;
295 *rsopsl = sopsl;
296
297 return ret;
298
299 err:
300 free(sops);
301 return ret;
302 }
303
findSymmetryPlanarIrregular(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)304 msym_error_t findSymmetryPlanarIrregular(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
305 msym_error_t ret = MSYM_SUCCESS;
306 int sopsl = 0;
307 msym_symmetry_operation_t *sops = NULL;
308
309 if(es->length != 4){
310 msymSetErrorDetails("Unexpected number of elements (%d) in planar irregular polygon",es->length);
311 ret = MSYM_SYMMETRY_ERROR;
312 goto err;
313 }
314
315 int iscm = vzero(cm,thresholds->zero);
316
317 if(iscm){
318 //3xC2 + 3xSigma + inversion
319 sopsl = 7;
320 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
321 } else {
322 // 2xSigma + 1xC2
323 sopsl = 3;
324 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
325 }
326
327 //The CM vector must be pointing in the same direction as the largest moment ov inertia vector, otherwise these would not be equal.
328
329 vcopy(ev[2],sops[0].v);
330 vnorm(sops[0].v);
331 sops[0].type = PROPER_ROTATION;
332 sops[0].order = 2;
333 sops[0].power = 1;
334
335 vcopy(ev[1],sops[1].v);
336 vnorm(sops[1].v);
337 sops[1].type = REFLECTION;
338
339 vcopy(ev[0],sops[2].v);
340 vnorm(sops[2].v);
341 sops[2].type = REFLECTION;
342
343 if(iscm){
344 vcopy(sops[0].v, sops[3].v);
345 sops[3].type = REFLECTION;
346
347 vcopy(sops[1].v, sops[4].v);
348 sops[4].type = PROPER_ROTATION;
349 sops[4].order = 2;
350 sops[4].power = 1;
351
352 vcopy(sops[2].v, sops[5].v);
353 sops[5].type = PROPER_ROTATION;
354 sops[5].order = 2;
355 sops[5].power = 1;
356
357 sops[6].type = INVERSION;
358 sops[6].order = 1;
359 sops[6].power = 1;
360 sops[6].v[0] = sops[6].v[1] = sops[6].v[2] = 0;
361
362 }
363
364 *rsopsl = sopsl;
365 *rsops = sops;
366
367 return ret;
368 err:
369 return ret;
370 }
371
findSymmetryPolyhedralProlate(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)372 msym_error_t findSymmetryPolyhedralProlate(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
373 return findSymmetrySymmetricPolyhedron(es,cm,ev,thresholds,0,rsopsl,rsops);
374 }
375
findSymmetryPolyhedralOblate(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)376 msym_error_t findSymmetryPolyhedralOblate(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
377 return findSymmetrySymmetricPolyhedron(es,cm,ev,thresholds,2,rsopsl,rsops);
378 }
379
findSymmetrySymmetricPolyhedron(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int prim,int * rsopsl,msym_symmetry_operation_t ** rsops)380 msym_error_t findSymmetrySymmetricPolyhedron(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int prim, int *rsopsl, msym_symmetry_operation_t **rsops){
381 msym_error_t ret = MSYM_SUCCESS;
382 int sopsl = 0;
383 msym_symmetry_operation_t *sops = NULL;
384
385 int sigma_h = 0, staggered = 0, split = 0, order = es->length/2, even, inversion, div_len;
386 int *div = NULL;
387
388 int n = 0;
389 double v0[3], v0_proj[3], v_init[3], dot0, theta, theta_C2 = 0.0, theta_sigma = 0.0;
390
391 if(!(vzero(cm,thresholds->zero))){
392 msymSetErrorDetails("Symmetric polyhedron not at center of mass. Vector length: %e > %e (zero threshold)", vabs(cm),thresholds->zero);
393 ret = MSYM_SYMMETRY_ERROR;
394 goto err;
395 }
396
397 vcopy(es->elements[0]->v,v0);
398 dot0 = vdot(v0,ev[prim]);
399
400 vproj_plane(v0, ev[prim], v0_proj);
401 vnorm(v0);
402 vnorm(v0_proj);
403 vcopy(v0_proj,v_init);
404
405
406 for(int i = 1; i < es->length;i++){
407 double vi[3], vi_proj[3], v0i[3], doti;
408 vcopy(es->elements[i]->v,vi);
409 doti = vdot(vi,ev[prim]);
410
411 vproj_plane(vi,ev[prim], vi_proj);
412 vnorm(vi);
413 vnorm(vi_proj);
414 vsub(v0,vi, v0i);
415 vnorm(v0i);
416
417 theta = vangle(v0_proj, vi_proj);
418
419 if(theta < asin(thresholds->angle)){
420 sigma_h = 1;
421 staggered = 0;
422 }
423
424 if(dot0*doti > 0.0){
425 theta_sigma = theta/2;
426 if(LT(theta, 4*M_PI/es->length, asin(thresholds->angle)) && es->length % 4 == 0){
427 vadd(v0,vi,v_init);
428 vproj_plane(v_init, ev[prim], v_init);
429 vnorm(v_init);
430 order = es->length/4;
431 even = order % 2 == 0;
432 split = 1;
433 //theta_split = theta;
434
435 }
436 } else {
437 //This can potentially find a staggered form if split by equal amounts (should be ok, TODO: check)
438 if(fabs(theta - 2*M_PI/es->length) < asin(thresholds->angle)){
439 staggered = 1;
440 }
441 // if we have not yet found that we are are staggerd/eclipsed or split this will either be over-written or this is a C2 axis for Dn
442 if(!split && !staggered && !sigma_h){
443 if(LT(theta, 2*M_PI/es->length, asin(thresholds->angle))){
444 vadd(v0_proj, vi_proj, v_init);
445 vnorm(v_init);
446 }
447 }
448 }
449 }
450
451 if(split){
452 staggered = !sigma_h;
453 }
454
455 even = order % 2 == 0;
456
457 inversion = (staggered && !even) || (sigma_h && even);
458 div = malloc(order*sizeof(int));
459 div_len = divisors(order,div);
460
461 //Symmetry operations:
462 //Primary rotation and all divisors
463 //sigma_h if eclipsed
464 //n x C2
465 //n x sigma_v if we are in staggered or eclipsed form
466 //Possible inversion (see above)
467 //Smax*2 if staggerd Sn (n > 2) for sigma_h
468
469 sopsl = (div_len + sigma_h + order + order*(sigma_h || staggered) + inversion + staggered + sigma_h*(div_len - even));
470 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
471
472 int max;
473 for(max = 0; n < div_len; n++){
474 if(div[n] > max){
475 max = div[n];
476 }
477 sops[n].type = PROPER_ROTATION;
478 sops[n].order = div[n];
479 sops[n].power = 1;
480 vcopy(ev[prim],sops[n].v);
481 }
482
483 if(sigma_h) {
484 sops[n].type = REFLECTION;
485 vcopy(ev[prim],sops[n].v);
486 n++;
487 for(int s = 0; s < div_len; s++){
488 if(div[s] > 2){
489 sops[n].type = IMPROPER_ROTATION;
490 sops[n].order = div[s];
491 sops[n].power = 1;
492 vcopy(ev[prim],sops[n].v);
493 n++;
494 }
495 }
496 }
497
498 if(inversion){
499 sops[n].type = INVERSION;
500 sops[n].order = 1;
501 sops[n].power = 1;
502 sops[n].v[0] = sops[n].v[1] = sops[n].v[2] = 0;
503 n++;
504 }
505
506 //PI/order angle between C2 & sigma_v
507 //PI/(2*order) start angle if split & staggered
508 //start angle if not split & staggered
509 //0 start angle if sigma_h
510
511 if(staggered){
512 theta_C2 = M_PI/(2*order);
513 sops[n].type = IMPROPER_ROTATION;
514 sops[n].order = 2*max;
515 sops[n].power = 1;
516 vcopy(ev[prim],sops[n].v);
517 n++;
518 }
519
520 theta = M_PI/order;
521 for (int i = 0; i < order; i++) {
522 vrotate(theta_C2 + i*theta, v_init, ev[prim], sops[n].v);
523 vnorm(sops[n].v); //Not really nessecary
524 sops[n].type = PROPER_ROTATION;
525 sops[n].order = 2;
526 sops[n].power = 1;
527 n++;
528
529 if(staggered || sigma_h){
530 vrotate(i*theta, v_init,ev[prim], sops[n].v);
531 vcrossnorm(sops[n].v,ev[prim],sops[n].v);
532 sops[n].type = REFLECTION;
533 n++;
534 }
535
536 }
537 if(n != sopsl){
538 msymSetErrorDetails("Unexpected number of generated symmetry operations in symmetric polyhedron. Got %d expected %d",n,sopsl);
539 ret = MSYM_SYMMETRY_ERROR;
540 goto err;
541 }
542
543 free(div);
544
545 *rsopsl = sopsl;
546 *rsops = sops;
547 return ret;
548
549 err:
550 free(div);
551 free(sops);
552 *rsops = NULL;
553 *rsopsl = 0;
554 return ret;
555
556 }
557
findSymmetryAsymmetricPolyhedron(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)558 msym_error_t findSymmetryAsymmetricPolyhedron(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
559 msym_error_t ret = MSYM_SUCCESS;
560 int sopsl = 0;
561 msym_symmetry_operation_t *sops = NULL;
562
563 if(es->length == 4){
564 //Only C2 axis
565 sopsl = 3;
566 } else if(es->length == 8){
567 sopsl = 7;
568 } else {
569 msymSetErrorDetails("Unexpected number of elements (%d) in asymmetric polyhedron",es->length);
570 ret = MSYM_SYMMETRY_ERROR;
571 goto err;
572 }
573
574 if(!(vzero(cm,thresholds->zero))){
575 msymSetErrorDetails("Asymmetric polyhedron not at center of mass. Vector length: %e > %e (zero threshold)", vabs(cm),thresholds->zero);
576 ret = MSYM_SYMMETRY_ERROR;
577 goto err;
578 }
579
580 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
581 vcopy(ev[0], sops[0].v);
582 vcopy(ev[1], sops[1].v);
583 vcopy(ev[2], sops[2].v);
584 vnorm(sops[0].v);
585 vnorm(sops[1].v);
586 vnorm(sops[2].v);
587 sops[0].type = PROPER_ROTATION;
588 sops[0].order = 2;
589 sops[0].power = 1;
590 sops[1].type = PROPER_ROTATION;
591 sops[1].order = 2;
592 sops[1].power = 1;
593 sops[2].type = PROPER_ROTATION;
594 sops[2].order = 2;
595 sops[2].power = 1;
596
597
598 if(es->length == 8){
599 vcopy(sops[0].v,sops[3].v);
600 vcopy(sops[1].v,sops[4].v);
601 vcopy(sops[2].v,sops[5].v);
602 sops[3].type = REFLECTION;
603 sops[3].order = 1;
604 sops[3].power = 1;
605 sops[4].type = REFLECTION;
606 sops[4].order = 1;
607 sops[4].power = 1;
608 sops[5].type = REFLECTION;
609 sops[5].order = 1;
610 sops[5].power = 1;
611 sops[6].type = INVERSION;
612 sops[6].order = 1;
613 sops[6].power = 1;
614 sops[6].v[0] = sops[6].v[1] = sops[6].v[2] = 0;
615
616 }
617
618 *rsopsl = sopsl;
619 *rsops = sops;
620 return ret;
621
622 err:
623 free(sops);
624 *rsops = NULL;
625 *rsopsl = 0;
626 return ret;
627
628 }
629
findSymmetrySpherical(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)630 msym_error_t findSymmetrySpherical(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
631 msym_error_t ret = MSYM_SUCCESS;
632 int sopsl = 0;
633 msym_symmetry_operation_t *sops = NULL;
634 if(es->length == 1){
635 if(!(vzero(cm,thresholds->zero))){
636 double t[3];
637 vcopy(es->elements[0]->v,t);
638 sopsl = 1;
639 sops = malloc(sopsl*sizeof(msym_symmetry_operation_t));
640 vcopy(t,sops[0].v);
641 vnorm(sops[0].v);
642 sops[0].type = PROPER_ROTATION;
643 sops[0].order = 0;
644 sops[0].power = 1;
645
646 }
647 *rsopsl = sopsl;
648 *rsops = sops;
649 } else {
650 ret = findSymmetryCubic(es,cm,ev,thresholds,rsopsl,rsops);
651 }
652
653 return ret;
654 //err:
655 // return ret;
656
657 }
658
659
findSymmetryCubic(msym_equivalence_set_t * es,double cm[3],double ev[3][3],msym_thresholds_t * thresholds,int * rsopsl,msym_symmetry_operation_t ** rsops)660 msym_error_t findSymmetryCubic(msym_equivalence_set_t *es, double cm[3], double ev[3][3], msym_thresholds_t *thresholds, int *rsopsl, msym_symmetry_operation_t **rsops){
661
662 msym_error_t ret = MSYM_SUCCESS;
663 /*struct _pair {
664 msym_element_t *e[2];
665 double d;
666 } *pairs = malloc(sizeof(struct _pair[150]));*/
667 double c2d = 0, c4d = 0, sigmad = 0;
668 int found = 0, nsigma = 0, esigma = 0, inversion = 0, nc[6] = {0,0,0,0,0,0}, ec[6] = {0,0,0,0,0,0}, *ncb = &nsigma, *nc3b = &nsigma;
669 msym_symmetry_operation_t **(ac[6]);
670 double thetac[6] = {0.0,M_PI,M_PI/2,M_PI/3,M_PI/4,M_PI/5};
671
672 msym_symmetry_operation_t *sops = malloc(sizeof(msym_symmetry_operation_t[120]));
673 int sopsl = 0;
674
675 double (**esv)[3] = malloc(sizeof(double (*[es->length])[3]));
676
677 msym_symmetry_operation_t **sigma = malloc(16*sizeof(msym_symmetry_operation_t*)); //only 15, but we can overflow
678
679 msym_symmetry_operation_t **c3b = NULL;
680 msym_symmetry_operation_t **cb = sigma;
681
682 msym_permutation_t perm;
683 msym_symmetry_operation_t sopc = {.type = PROPER_ROTATION, .order = 2, .power = 1};
684 msym_symmetry_operation_t sopsigma = {.type = REFLECTION, .order = 1, .power = 1};
685 msym_symmetry_operation_t sopinversion = {.type = INVERSION, .order = 1, .power = 1, .v = {0,0,0}};
686
687 //15 C2 (I), 10 C3 (I), 3 C4 (O), 6 C5 (I)
688 ac[0] = malloc((15+10+3+6)*sizeof(msym_symmetry_operation_t*));
689 ac[1] = ac[0];
690 ac[2] = ac[1];
691 ac[3] = ac[2] + 15;
692 ac[4] = ac[3] + 10;
693 ac[5] = ac[4] + 3;
694
695 for(int i = 0; i < es->length;i++){
696 esv[i] = &es->elements[i]->v;
697 }
698
699 for(int i = 0; i < es->length && !found;i++){
700 for(int j = i+1;j < es->length;j++){
701 if(vparallel(es->elements[i]->v, es->elements[j]->v, thresholds->angle)) continue;
702 double d;
703 int psopsl = sopsl;
704 vsub(es->elements[i]->v,es->elements[j]->v,sopsigma.v);
705 d = vabs(sopsigma.v);
706 vadd(es->elements[i]->v,es->elements[j]->v,sopc.v);
707 vnorm(sopc.v);
708 //vsub(es->elements[i]->v,es->elements[j]->v,sopsigma.v);
709 vnorm(sopsigma.v);
710
711 //Do this first, otherwise we might find a c2, at this distance and not look for more.
712 if(es->length % 5 != 0 && (c4d == 0 || fabs(d - c4d)/(d + c4d) < thresholds->equivalence)){
713 sopc.order = 4;
714 if(!findSymmetryOperation(&sopc, sops, sopsl, thresholds)){
715 if(MSYM_SUCCESS == findPermutation(&sopc, es->length, esv, thresholds, &perm)){
716 c4d = d;
717 copySymmetryOperation(&sops[sopsl], &sopc);
718 (ac[4])[nc[4]++] = &sops[sopsl++];
719 free(perm.p);
720 free(perm.c);
721 sopc.order = 2;
722 if(!findSymmetryOperation(&sopc, sops, sopsl, thresholds)){
723 copySymmetryOperation(&sops[sopsl], &sopc);
724 (ac[2])[nc[2]++] = &sops[sopsl++];
725
726 }
727 }
728
729 if(!findSymmetryOperation(&sopsigma, sops, sopsl, thresholds)){
730 if(MSYM_SUCCESS == findPermutation(&sopsigma, es->length, esv, thresholds, &perm)){
731 //sigmad = d; //This is a bit dangerous, but the C2 axes dhould generate the rest
732 copySymmetryOperation(&sops[sopsl], &sopsigma);
733 sigma[nsigma++] = &sops[sopsl++];
734 free(perm.p);
735 free(perm.c);
736 }
737 }
738 }
739 sopc.order = 2;
740 }
741
742 if(c2d == 0 || fabs(d - c2d)/(d + c2d) < thresholds->equivalence){
743 if(!findSymmetryOperation(&sopc, sops, sopsl, thresholds)){
744 if(MSYM_SUCCESS == findPermutation(&sopc, es->length, esv, thresholds, &perm)){
745 c2d = d;
746 copySymmetryOperation(&sops[sopsl], &sopc);
747 (ac[2])[nc[2]++] = &sops[sopsl++];
748 free(perm.p);
749 free(perm.c);
750 }
751 }
752 }
753
754 if(sigmad == 0 || fabs(d - sigmad)/(d + sigmad) < thresholds->equivalence){
755 if(!findSymmetryOperation(&sopsigma, sops, sopsl, thresholds)){
756 if(MSYM_SUCCESS == findPermutation(&sopsigma, es->length, esv, thresholds, &perm)){
757 sigmad = d;
758
759 copySymmetryOperation(&sops[sopsl], &sopsigma);
760 sigma[nsigma++] = &sops[sopsl++];
761 free(perm.p);
762 free(perm.c);
763 }
764 }
765 }
766
767 //we have generated something, we have more than 2 operations
768 if(sopsl > psopsl && sopsl >= 2 && !(nsigma == 15 && nc[2] == 0) && !(nsigma == 0 && nc[2] == 15)){
769
770 for(msym_symmetry_operation_t *sopi = sops; sopi < (sops + sopsl);sopi++){
771 for(msym_symmetry_operation_t *sopj = sops; sopj < (sops + sopsl); sopj++){
772 if(sopi == sopj) continue;
773 if(sopi->type == PROPER_ROTATION && sopj->type == PROPER_ROTATION && sopj->order == 2 && vperpendicular(sopi->v, sopj->v, thresholds->angle)){
774 copySymmetryOperation(&sops[sopsl], sopj);
775 vcrossnorm(sopi->v, sopj->v, sops[sopsl].v);
776 if(!findSymmetryOperation(&sops[sopsl], sops, sopsl,thresholds)){
777 (ac[2])[nc[2]++] = &sops[sopsl++];
778 }
779
780 } else if(!vparallel(sopi->v, sopj->v,thresholds->angle)){
781 copySymmetryOperation(&sops[sopsl], sopj);
782 applySymmetryOperation(sopi,sops[sopsl].v,sops[sopsl].v);
783 if(!findSymmetryOperation(&sops[sopsl], sops, sopsl,thresholds)){
784 if(sopj->type == REFLECTION){
785 sigma[nsigma++] = &sops[sopsl++];
786 } else {
787 (ac[sopj->order])[nc[sopj->order]++] = &sops[sopsl++];
788 }
789 }
790 }
791 if(nsigma > 15 || nc[2] > 15){
792 ret = MSYM_SYMMETRY_ERROR;
793 msymSetErrorDetails("Inconsistency when generating symmetry axes for cubic point group, %d C2 axes, %d reflection planes",nc[2],nsigma);
794 goto err;
795 }
796 }
797 }
798 }
799
800 //A little speedup, I point group may require up to 5 iterations, to check if it might be an I
801 if(((nsigma >= 3 || (nc[2] >= 3 && nsigma == 0)) && nc[4] == 0 && i > 3 && es->length % 5 != 0) ||
802 ((nsigma >= 9 || (nc[2] >= 9 && nsigma == 0)) && i > 0) ||
803 (nsigma == 15 && nc[2] == 15)){
804 found = 1;
805 break;
806 }
807 if(nsigma == 0 && nc[2] == 0 && nc[4] == 0 && i > 3 && es->length > 120){
808 ret = MSYM_SYMMETRY_ERROR;
809 msymSetErrorDetails("Found no symmetry operations in cubic group of size %d, thresholds are too high ",es->length);
810 goto err;
811 }
812 }
813 }
814
815 if(MSYM_SUCCESS == findPermutation(&sopinversion, es->length, esv, thresholds, &perm)){
816 inversion = 1;
817 copySymmetryOperation(&sops[sopsl++], &sopinversion);
818 free(perm.p);
819 free(perm.c);
820 }
821
822 esigma = (((nsigma+2) / 3) + (nsigma / 10) - (nsigma / 13))*3;
823 esigma = esigma == 1 ? 0 : esigma; //We cannot generate the remaining axes from one sigma try to generate from C2
824
825 switch(esigma) {
826 case 15 : ec[5] = 6; ec[3] = 10; ec[2] = 15; break; //Ih
827 case 9 : ec[4] = 3; ec[3] = 4; ec[2] = 9; break; //Oh
828 case 6 :
829 { //Td or Oh (a cube won't generate the pair reflection planes, but we can look for inversion)
830 if(inversion){
831 int gsigma = 0;
832 ec[4] = 3; ec[3] = 4; ec[2] = 9;
833 for(int i = 0; i < nsigma && gsigma < 3;i++){
834 for(int j = i+1; j < nsigma && gsigma < 3;j++){
835 double theta = vangle(sigma[i]->v, sigma[j]->v);
836 if(fabs(theta - M_PI/2) < asin(thresholds->angle)){
837 sops[sopsl].type = REFLECTION;
838 vadd(sigma[i]->v, sigma[j]->v, sops[sopsl].v);
839 vnorm(sops[sopsl].v);
840 if(!findSymmetryOperation(&(sops[sopsl]), sops, sopsl,thresholds)){
841 sigma[nsigma+gsigma++] = &(sops[sopsl++]);
842 //gsigma++;
843 //sopsl++;
844 }
845 if(gsigma < 3){
846 sops[sopsl].type = REFLECTION;
847 vsub(sigma[i]->v, sigma[j]->v, sops[sopsl].v);
848 vnorm(sops[sopsl].v);
849 if(!findSymmetryOperation(&(sops[sopsl]), sops, sopsl,thresholds)){
850 sigma[nsigma+gsigma++] = &(sops[sopsl++]);
851 //gsigma++;
852 //sopsl++;
853 }
854 }
855 }
856 }
857 }
858 nsigma += gsigma;
859 break;
860 } else {
861 ec[3] = 4; ec[2] = 3; break;
862 }
863 }
864 case 3 : ec[3] = 4; ec[2] = 3; c3b = sigma; nc3b = &nsigma; break; //Th
865 case 0 :
866 {
867 ec[2] = ((nc[2] > 0) + ((nc[2] > 3) << 1) + ((nc[2] > 9) << 1))*3;
868 switch(ec[2]){
869 case 3 : ec[3] = 4; ec[2] = 3; c3b = ac[2]; nc3b = &(nc[2]); break;
870 case 9 :
871 ec[4] = 3; ec[3] = 4; ec[2] = 9;
872 for(int i = 0; i < nc[2] && nc[4] < 3; i++){ //This should not happen anymore, but it may segfault on an error, so keep it for now
873 int p = 0;
874 for(int j = 0; j < nc[2] && nc[4] < 3; j++){
875 p += vperpendicular((ac[2])[i]->v, (ac[2])[j]->v,thresholds->angle);
876 }
877
878 if(p == 4){
879 sops[sopsl].type = PROPER_ROTATION;
880 sops[sopsl].order = 4;
881 sops[sopsl].power = 1;
882 vcopy((ac[2])[i]->v,sops[sopsl].v);
883 if(!findSymmetryOperation(&sops[sopsl], sops, sopsl,thresholds)){
884 (ac[4])[nc[4]] = &(sops[sopsl++]);
885 nc[4]++;
886 }
887 }
888 }
889 c3b = ac[4];
890 nc3b = &(nc[4]);
891 break;
892 //Look for C2 that have 4 other perpendicular C2 this will be a C4
893 case 15 :
894 ec[5] = 6; ec[3] = 10; ec[2] = 15;
895 ncb = &(nc[2]);
896 cb = ac[2];
897 break;
898 default :
899 msymSetErrorDetails("Unexpected number of C2 axes (%d) in cubic point group",ec[2]);
900 ret = MSYM_SYMMETRY_ERROR;
901 goto err;
902
903 }
904 break;
905 }
906
907 default :
908 msymSetErrorDetails("Unexpected number of reflection planes (%d) in cubic point group",nsigma);
909 ret = MSYM_SYMMETRY_ERROR;
910 goto err;
911 }
912
913
914 //In a Th, we might not find the C3 axes so we generated them from the 3 reflection planes.
915 //In the T, O, groups we generate them from C2 and C4 axes resp.
916 //In I they will be generated later by the c2 axis (cbase = C2 in that case)
917 if(c3b != NULL && nc[3] == 0 && *nc3b >= 3){
918 vadd(c3b[0]->v,c3b[1]->v,sops[sopsl].v);
919 vadd(sops[sopsl].v,c3b[2]->v,sops[sopsl].v);
920 vnorm(sops[sopsl].v);
921 sops[sopsl].type = PROPER_ROTATION;
922 sops[sopsl].order = 3;
923 sops[sopsl].power = 1;
924 (ac[3])[nc[3]] = &(sops[sopsl++]);
925 //sopsl++;
926 nc[3]++;
927
928 vadd(c3b[0]->v,c3b[1]->v,sops[sopsl].v);
929 vsub(sops[sopsl].v,c3b[2]->v,sops[sopsl].v);
930 vnorm(sops[sopsl].v);
931 sops[sopsl].type = PROPER_ROTATION;
932 sops[sopsl].order = 3;
933 sops[sopsl].power = 1;
934 (ac[3])[nc[3]] = &(sops[sopsl++]);
935 //sopsl++;
936 nc[3]++;
937
938 vadd(c3b[0]->v,c3b[2]->v,sops[sopsl].v);
939 vsub(sops[sopsl].v,c3b[1]->v,sops[sopsl].v);
940 vnorm(sops[sopsl].v);
941 sops[sopsl].type = PROPER_ROTATION;
942 sops[sopsl].order = 3;
943 sops[sopsl].power = 1;
944 (ac[3])[nc[3]] = &(sops[sopsl++]);
945 //sopsl++;
946 nc[3]++;
947
948 vsub(c3b[0]->v,c3b[1]->v,sops[sopsl].v);
949 vsub(sops[sopsl].v,c3b[2]->v,sops[sopsl].v);
950 vnorm(sops[sopsl].v);
951 sops[sopsl].type = PROPER_ROTATION;
952 sops[sopsl].order = 3;
953 sops[sopsl].power = 1;
954 (ac[3])[nc[3]] = &(sops[sopsl++]);
955 //sopsl++;
956 nc[3]++;
957 }
958
959 found = nc[2] == ec[2] && nc[3] == ec[3] && nc[4] == ec[4] && nc[5] == ec[5];
960
961 for(int i = 0; i < *ncb && !found && sopsl < 120;i++){
962 for(int j = i+1; j < *ncb && !found;j++){
963 double theta = fabs(vangle(cb[i]->v, cb[j]->v));
964 if(theta > M_PI/2){
965 theta = M_PI - theta;
966 }
967 for(int k = 2; k < 6;k++){
968
969 if(fabs(theta - thetac[k]) < asin(thresholds->angle)){
970 vcrossnorm(cb[i]->v, cb[j]->v, sops[sopsl].v);
971 sops[sopsl].type = PROPER_ROTATION;
972 sops[sopsl].order = k;
973 sops[sopsl].power = 1;
974 if(!findSymmetryOperation(&(sops[sopsl]), sops, sopsl,thresholds)){
975 (ac[k])[nc[k]] = &(sops[sopsl]);
976 sopsl++;
977 nc[k]++;
978 }
979 break;
980 }
981 }
982 found = nc[2] == ec[2] && nc[3] == ec[3] && nc[4] == ec[4] && nc[5] == ec[5];
983 }
984 }
985
986 int nS = 0;
987 if(nsigma > 0){
988 for(int i = 0; i < sopsl && sopsl < 120;i++){
989 if(sops[i].type == PROPER_ROTATION){
990 //Td
991 if(!inversion && sops[i].order == 2){
992 vcopy(sops[i].v, sops[sopsl].v);
993 sops[sopsl].type = IMPROPER_ROTATION;
994 sops[sopsl].order = 4;
995 sops[sopsl].power = 1;
996 sopsl++;
997 nS++;
998 }
999 else if (inversion && sops[i].order > 2) {
1000 vcopy(sops[i].v, sops[sopsl].v);
1001 sops[sopsl].type = IMPROPER_ROTATION;
1002 sops[sopsl].order = sops[i].order + (sops[i].order % 2)*sops[i].order;
1003 sops[sopsl].power = 1;
1004 sopsl++;
1005 nS++;
1006 }
1007 }
1008 }
1009 }
1010
1011 if(sopsl > 63) {
1012 msymSetErrorDetails("Generated more than 63 symmetry operations (%d) in cubic point group, not including powers",sopsl);
1013 ret = MSYM_SYMMETRY_ERROR;
1014 goto err;
1015 }
1016
1017 *rsopsl = sopsl;
1018 *rsops = sops;
1019
1020 free(esv);
1021 free(ac[0]);
1022 free(sigma);
1023
1024 return ret;
1025
1026 err:
1027 free(ac[0]);
1028 free(sigma);
1029 free(sops);
1030 free(esv);
1031 *rsopsl = 0;
1032 *rsops = NULL;
1033 return ret;
1034 }
1035
reduceSymmetry(int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,int * isopsl,msym_symmetry_operation_t ** isops)1036 msym_error_t reduceSymmetry(int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, int *isopsl, msym_symmetry_operation_t **isops){
1037 msym_error_t ret = MSYM_SUCCESS;
1038
1039 int rsopsl = *isopsl;
1040 msym_symmetry_operation_t *rsops = *isops;
1041 msym_symmetry_operation_t *cinf[2] = {NULL,NULL};
1042
1043 int inv[2] = {0,0};
1044 int inversion = 0;
1045 for(int i = 0;i < rsopsl;i++){
1046 if(!((rsops[i].type == PROPER_ROTATION && rsops[i].order == 0) || rsops[i].type == INVERSION || rsops[i].type == REFLECTION)) break;
1047 if(inv[0] && cinf[0] != NULL) break;
1048 inv[0] = inv[0] || rsops[i].type == INVERSION;
1049 if(rsops[i].type == PROPER_ROTATION && rsops[i].order == 0) cinf[0] = &rsops[i];
1050 }
1051
1052 for(int i = 0;i < sopsl;i++){
1053 if(!((sops[i].type == PROPER_ROTATION && sops[i].order == 0) || sops[i].type == INVERSION || sops[i].type == REFLECTION)) break;
1054 if(inv[1] && cinf[1] != NULL) break;
1055 inv[1] = inv[1] || sops[i].type == INVERSION;
1056 if(sops[i].type == PROPER_ROTATION && sops[i].order == 0) cinf[1] = &sops[i];
1057 }
1058
1059 inversion = inv[0] && inv[1];
1060
1061 if(cinf[0] != NULL && cinf[1] != NULL){
1062 double cross[3];
1063 int perpendicular = vperpendicular(cinf[0]->v,cinf[1]->v,thresholds->angle);
1064 int parallel = vparallel(cinf[0]->v,cinf[1]->v,thresholds->angle);
1065 vcrossnorm(cinf[0]->v, cinf[1]->v, cross);
1066 if(inversion && perpendicular){
1067 double v[3];
1068 vcopy(cinf[0]->v, v);
1069 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[7]));
1070 rsopsl = 7;
1071
1072 rsops[0].type = INVERSION;
1073 rsops[0].v[0] = rsops[0].v[1] = rsops[0].v[2] = 0;
1074 rsops[0].order = 1;
1075 rsops[0].power = 1;
1076 rsops[1].type = REFLECTION;
1077 rsops[1].order = 1;
1078 rsops[1].power = 1;
1079 rsops[2].type = REFLECTION;
1080 rsops[2].order = 1;
1081 rsops[2].power = 1;
1082 rsops[3].type = REFLECTION;
1083 rsops[3].order = 1;
1084 rsops[3].power = 1;
1085 rsops[4].type = PROPER_ROTATION;
1086 rsops[4].order = 2;
1087 rsops[4].power = 1;
1088 rsops[5].type = PROPER_ROTATION;
1089 rsops[5].order = 2;
1090 rsops[5].power = 1;
1091 rsops[6].type = PROPER_ROTATION;
1092 rsops[6].order = 2;
1093 rsops[6].power = 1;
1094
1095 vcopy(v, rsops[1].v);
1096 vcopy(cinf[1]->v, rsops[2].v);
1097 vcopy(cross, rsops[3].v);
1098 vcopy(v, rsops[4].v);
1099 vcopy(cinf[1]->v, rsops[5].v);
1100 vcopy(cross, rsops[6].v);
1101 } else if (inversion & !parallel){
1102 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[3]));
1103 rsopsl = 3;
1104 rsops[0].type = INVERSION;
1105 rsops[0].v[0] = rsops[0].v[1] = rsops[0].v[2] = 0;
1106 rsops[0].order = 1;
1107 rsops[0].power = 1;
1108 rsops[1].type = REFLECTION;
1109 rsops[1].order = 1;
1110 rsops[1].power = 1;
1111 rsops[2].type = PROPER_ROTATION;
1112 rsops[2].order = 2;
1113 rsops[2].power = 1;
1114
1115 vcopy(cross, rsops[1].v);
1116 vcopy(cross, rsops[2].v);
1117 } else if(perpendicular && !inv[0] && !inv[1]){
1118 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[1]));
1119 rsopsl = 1;
1120 rsops[0].type = REFLECTION;
1121 vcopy(cross, rsops[0].v);
1122 } else if (perpendicular){
1123 int index = inv[0] ? 0 : 1;
1124 double v[2][3];
1125 vcopy(cinf[0]->v, v[0]);
1126 vcopy(cinf[1]->v, v[1]);
1127
1128 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[3]));
1129 rsopsl = 3;
1130 rsops[0].type = REFLECTION;
1131 rsops[1].type = REFLECTION;
1132 rsops[2].type = PROPER_ROTATION;
1133 rsops[2].order = 2;
1134 rsops[2].power = 1;
1135
1136 vcopy(cross, rsops[0].v);
1137 vcopy(v[index], rsops[1].v);
1138 vcopy(v[(index+1)%2], rsops[2].v);
1139
1140 } else if (parallel){
1141 if(MSYM_SUCCESS != (ret = filterSymmetryOperations(sopsl,sops,thresholds,&rsopsl,&rsops))) goto err;
1142 } else {
1143 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[1]));
1144 rsopsl = 1;
1145 rsops[0].type = REFLECTION;
1146 vcopy(cross, rsops[0].v);
1147 }
1148 } else if (cinf[0] != NULL){
1149 double v[3];
1150 vcopy(cinf[0]->v, v);
1151 for(int i = 0;i < sopsl;i++){
1152 int add = 0;
1153 if(sops[i].type == IMPROPER_ROTATION){
1154 add = vparallel(sops[i].v,v,thresholds->angle);
1155 } else if(sops[i].type == PROPER_ROTATION){
1156 if(sops[i].order != 2){
1157 add = vparallel(sops[i].v,v,thresholds->angle);
1158 } else {
1159 add = vparallel(sops[i].v,v,thresholds->angle) || (vperpendicular(sops[i].v,cinf[0]->v,thresholds->angle) && inv[0]);
1160 }
1161 } else if(sops[i].type == REFLECTION){
1162 add = vperpendicular(sops[i].v,v,thresholds->angle);
1163 }
1164 if(add){
1165 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[rsopsl+1]));
1166 copySymmetryOperation(&rsops[rsopsl], &sops[i]);
1167 rsopsl++;
1168 }
1169 }
1170 if(MSYM_SUCCESS != (ret = filterSymmetryOperations(sopsl,sops,thresholds,&rsopsl,&rsops))) goto err;
1171 } else if (cinf[1] != NULL){
1172 for(int i = 0;i < rsopsl && rsopsl > 0;i++){
1173 msym_symmetry_operation_t *fsop = findSymmetryOperation(&rsops[i], sops, sopsl,thresholds);
1174 if(!fsop){
1175 int remove = 1;
1176 if(rsops[i].type == IMPROPER_ROTATION){
1177 remove = !vparallel(rsops[i].v,cinf[1]->v,thresholds->angle);
1178 } else if(rsops[i].type == PROPER_ROTATION){
1179 if(rsops[i].order != 2){
1180 remove = !vparallel(rsops[i].v,cinf[1]->v,thresholds->angle);
1181 } else {
1182 remove = !(vparallel(rsops[i].v,cinf[1]->v,thresholds->angle) || (vperpendicular(rsops[i].v,cinf[1]->v,thresholds->angle) && inv[1]));
1183 }
1184 } else if(rsops[i].type == REFLECTION){
1185 remove = !vperpendicular(rsops[i].v,cinf[1]->v,thresholds->angle);
1186 }
1187 if(remove){
1188 rsopsl--;
1189 copySymmetryOperation(&rsops[i], &rsops[rsopsl]);
1190 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[rsopsl]));
1191 i--;
1192 } else if(vparallel(rsops[i].v,cinf[1]->v,thresholds->angle)){
1193 if(vdot(rsops[i].v,cinf[1]->v) < 0){
1194 vsub(rsops[i].v,cinf[1]->v,rsops[i].v);
1195 } else {
1196 vadd(rsops[i].v,cinf[1]->v,rsops[i].v);
1197 }
1198 }
1199 }
1200 }
1201 } else {
1202 if(MSYM_SUCCESS != (ret = filterSymmetryOperations(sopsl,sops,thresholds, &rsopsl,&rsops))) goto err;
1203 }
1204
1205 *isopsl = rsopsl;
1206 *isops = rsops;
1207 return ret;
1208 err:
1209 return ret;
1210 }
1211
filterSymmetryOperations(int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,int * isopsl,msym_symmetry_operation_t ** isops)1212 msym_error_t filterSymmetryOperations(int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, int *isopsl, msym_symmetry_operation_t **isops){
1213 msym_error_t ret = MSYM_SUCCESS;
1214 int rsopsl = *isopsl;
1215 msym_symmetry_operation_t *rsops = *isops;
1216
1217 for(int i = 0;i < rsopsl && rsopsl > 0;i++){
1218 msym_symmetry_operation_t *fsop = findSymmetryOperation(&rsops[i], sops, sopsl,thresholds);
1219 if(!fsop){
1220 rsopsl--;
1221 copySymmetryOperation(&rsops[i], &rsops[rsopsl]);
1222 rsops = realloc(rsops, sizeof(msym_symmetry_operation_t[rsopsl]));
1223 i--;
1224 } else if (rsops[i].type == PROPER_ROTATION || rsops[i].type == IMPROPER_ROTATION || rsops[i].type == REFLECTION){
1225 if(vdot(rsops[i].v,fsop->v) < 0){
1226 vsub(rsops[i].v,fsop->v,rsops[i].v);
1227 } else {
1228 vadd(rsops[i].v,fsop->v,rsops[i].v);
1229 }
1230 }
1231 }
1232
1233 *isopsl = rsopsl;
1234 *isops = rsops;
1235 return ret;
1236 //err:
1237 // return ret;
1238
1239 }
1240
divisors(int n,int * div)1241 int divisors(int n, int* div){
1242 int max = floor(sqrt(n));
1243 div[0] = n;
1244 int l = 1;
1245 for(int i = 2; i <= max;i++) {
1246 if(n % i == 0){
1247 int fact = n/i;
1248 div[l++] = i;
1249 if(i != fact){
1250 div[l++] = n/i;
1251 }
1252 }
1253 }
1254 return l;
1255 }
1256