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