<commit_before>//
//  point_group.c
//  Symmetry
//
//  Created by <NAME> on 14/04/14.
//  Copyright (c) 2014 <NAME>. 
//
//  Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
//


#include <math.h>
#include <float.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "symmetry.h"
#include "linalg.h"
#include "point_group.h"
#include "permutation.h"

#define PHI 1.618033988749894848204586834
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288419716939937510582
#endif

#define CLASSIFICATION_THRESHOLD 0.01


msym_error_t symmetrizePointGroup(msym_point_group_t *ipg, msym_point_group_t **opg, msym_thresholds_t *thresholds);
msym_error_t copyPointGroup(msym_point_group_t *ipg, msym_point_group_t *opg);
msym_error_t reorientAxes(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t transformPrimary(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t transformSecondary(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t setPointGroupOrder(msym_point_group_t *pg);
msym_error_t setPointGroupName(size_t max, int n, msym_point_group_type_t type, char *name);

msym_error_t findSecondaryAxisSigma(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC2(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC3(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC4(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC2C5(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);


msym_error_t generateSymmetryOperations(msym_point_group_t *pg, msym_thresholds_t *thresholds);

msym_error_t generateSymmetryOperationsImpliedS(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t generateSymmetryOperationsImpliedCPow(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t generateSymmetryOperationsImpliedSPow(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t generateSymmetryOperationsImpliedRot(msym_point_group_t *pg, msym_thresholds_t *thresholds);

void generateSymmetryOperationsCi(msym_point_group_t *pg);
void generateSymmetryOperationsCs(msym_point_group_t *pg);
void generateSymmetryOperationsCn(msym_point_group_t *pg);
void generateSymmetryOperationsCnh(msym_point_group_t *pg);
void generateSymmetryOperationsCnv(msym_point_group_t *pg);
void generateSymmetryOperationsDn(msym_point_group_t *pg);
void generateSymmetryOperationsDnh(msym_point_group_t *pg);
void generateSymmetryOperationsDnd(msym_point_group_t *pg);
void generateSymmetryOperationsS2n(msym_point_group_t *pg);
void generateSymmetryOperationsT(msym_point_group_t *pg);
void generateSymmetryOperationsTd(msym_point_group_t *pg);
void generateSymmetryOperationsTh(msym_point_group_t *pg);
void generateSymmetryOperationsO(msym_point_group_t *pg);
void generateSymmetryOperationsOh(msym_point_group_t *pg);
void generateSymmetryOperationsI(msym_point_group_t *pg);
void generateSymmetryOperationsIh(msym_point_group_t *pg);

void addSymmetryOperationsCi(msym_point_group_t *pg);
void addSymmetryOperationsCs(msym_point_group_t *pg);
void addSymmetryOperationsCn(msym_point_group_t *pg);
void addSymmetryOperationsCnh(msym_point_group_t *pg);
void addSymmetryOperationsCnv(msym_point_group_t *pg);
void addSymmetryOperationsDn(msym_point_group_t *pg);
void addSymmetryOperationsDnh(msym_point_group_t *pg);
void addSymmetryOperationsDnd(msym_point_group_t *pg);
void addSymmetryOperationsS2n(msym_point_group_t *pg);
void addSymmetryOperationsT(msym_point_group_t *pg);
void addSymmetryOperationsTd(msym_point_group_t *pg);
void addSymmetryOperationsTh(msym_point_group_t *pg);
void addSymmetryOperationsO(msym_point_group_t *pg);
void addSymmetryOperationsOh(msym_point_group_t *pg);
void addSymmetryOperationsI(msym_point_group_t *pg);
void addSymmetryOperationsIh(msym_point_group_t *pg);

msym_error_t pointGroupFromName(char *name, msym_point_group_t *pg);

int classifySymmetryOperations(msym_point_group_t *pg);
void sortSymmetryOperations(msym_point_group_t *pg, int classes);

void print_transform(double M[3][3], double axis[3]);



msym_error_t generatePointGroup(char *name, msym_thresholds_t *thresholds, msym_point_group_t **opg){
    msym_error_t ret = MSYM_SUCCESS;
    msym_point_group_t *pg = calloc(1,sizeof(msym_point_group_t));
    if(MSYM_SUCCESS != (ret = pointGroupFromName(name,pg))) goto err;
    if(MSYM_SUCCESS != (ret = generateSymmetryOperations(pg,thresholds))) goto err;
    int classes = classifySymmetryOperations(pg);
    sortSymmetryOperations(pg,classes);
    
    if((pg->type == POINT_GROUP_Cnv && pg->n == 0) || (pg->type == POINT_GROUP_Dnh && pg->n == 0)){
        pg->perm = NULL;
    } else {
        if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations(pg->sopsl,pg->sops, thresholds, &pg->perm))) goto err;
    }
    for(msym_symmetry_operation_t *s = pg->sops;s < (pg->sops + pg->sopsl);s++){
        if(pg->primary == NULL || (s->type == PROPER_ROTATION && s->order > pg->primary->order)) pg->primary = s;
    }
    
    mleye(3,pg->transform);
    
    *opg = pg;
    return ret;
    
err:
    *opg = NULL;
    free(pg);
    return ret;
}


msym_error_t pointGroupFromName(char *name, msym_point_group_t *pg){
    msym_error_t ret = MSYM_SUCCESS;
    int n = 0, gi = 0, ri = 0;;
    char g = 0, r = 0;
    
    int map[7][6];
    
    
    struct _pg_map {
        int i;
        msym_point_group_type_t type;
    } pg_map[] = {
        {1,  POINT_GROUP_Cn},
        {2,  POINT_GROUP_Cnv},
        {3,  POINT_GROUP_Cnh},
        {4,  POINT_GROUP_Ci},
        {5,  POINT_GROUP_Cs},
        {6,  POINT_GROUP_Dn},
        {7,  POINT_GROUP_Dnh},
        {8,  POINT_GROUP_Dnd},
        {9,  POINT_GROUP_S2n},
        {10, POINT_GROUP_T},
        {11, POINT_GROUP_Th},
        {12, POINT_GROUP_Td},
        {13, POINT_GROUP_O},
        {14, POINT_GROUP_Oh},
        {15, POINT_GROUP_I},
        {16, POINT_GROUP_Ih},
        {17, POINT_GROUP_K},
        {18, POINT_GROUP_Kh}
        
    };
    
    memset(map,0,sizeof(int[7][6]));
    
    map[0][0] =  1;
    map[0][1] =  2;
    map[0][2] =  3;
    map[0][4] =  4;
    map[0][5] =  5;
    map[1][0] =  6;
    map[1][2] =  7;
    map[1][3] =  8;
    map[2][0] =  9;
    map[3][0] =  10;
    map[3][2] =  11;
    map[3][3] =  12;
    map[4][0] =  13;
    map[4][2] =  14;
    map[5][0] =  15;
    map[5][2] =  16;
    map[6][0] =  17;
    map[6][2] =  18;
    
    if(sscanf(name,"%c%d%c",&g,&n,&r) < 2 && sscanf(name,"%c%c",&g,&r) < 1){
        ret = MSYM_INVALID_POINT_GROUP;
        msymSetErrorDetails("Invalid point group name %s",name);
        goto err;
    }
    
    if(n < 0) {
        ret = MSYM_INVALID_POINT_GROUP;
        msymSetErrorDetails("Invalid point group order %d",n);
        goto err;
    }
    
    switch(g){
        case 'C' : gi = 0; break;
        case 'D' : gi = 1; break;
        case 'S' : {
            if(n < 4 || n % 2 != 0){
                ret = MSYM_INVALID_POINT_GROUP;
                msymSetErrorDetails("Improper rotation order (%d) must be even",n);
                goto err;
            }
            gi = 2;
            break;
        }
        case 'T' : gi = 3; break;
        case 'O' : gi = 4; break;
        case 'I' : gi = 5; break;
        case 'K' : gi = 6; break;
        default :
            ret = MSYM_INVALID_POINT_GROUP;
            msymSetErrorDetails("Invalid point group type %c",g);
            goto err;
    }
    
    switch(r){
        case 0   : ri = 0; break;
        case 'v' : ri = 1; break;
        case 'h' : ri = 2; break;
        case 'd' : ri = 3; break;
        case 'i' : ri = 4; break;
        case 's' : ri = 5; break;
        default :
            ret = MSYM_INVALID_POINT_GROUP;
            msymSetErrorDetails("Invalid point group subtype %c",r);
            goto err;
    }
    
    int fi, fil = sizeof(pg_map)/sizeof(pg_map[0]);
    for(fi = 0;fi < fil;fi++){
        if(pg_map[fi].i == map[gi][ri]) break;
    }
    
    if(fi == fil){
        ret = MSYM_INVALID_POINT_GROUP;
        msymSetErrorDetails("Cannot find point group %s",name);
        goto err;
    }
    
    pg->type = pg_map[fi].type;
    pg->n = n;
    if(MSYM_SUCCESS != (ret = setPointGroupOrder(pg))) goto err;
    if(MSYM_SUCCESS != (ret = setPointGroupName(sizeof(pg->name)/sizeof(char),pg->n,pg->type,pg->name))) goto err;
    
err:
    return ret;
    
}


msym_error_t setPointGroupName(size_t max, int n, msym_point_group_type_t type, char *name){
    msym_error_t ret = MSYM_SUCCESS;
    switch(type) {
        case POINT_GROUP_Ci  : snprintf(name,max,"Ci"); break;
        case POINT_GROUP_Cs  : snprintf(name,max,"Cs"); break;
        case POINT_GROUP_Cn  : snprintf(name,max,"C%d",n); break;
        case POINT_GROUP_Cnh : snprintf(name,max,"C%dh",n); break;
        case POINT_GROUP_Cnv : snprintf(name,max,"C%dv",n); break;
        case POINT_GROUP_Dn  : snprintf(name,max,"D%d",n); break;
        case POINT_GROUP_Dnh : snprintf(name,max,"D%dh",n); break;
        case POINT_GROUP_Dnd : snprintf(name,max,"D%dd",n); break;
        case POINT_GROUP_S2n : snprintf(name,max,"S%d",n); break;
        case POINT_GROUP_T   : snprintf(name,max,"T"); break;
        case POINT_GROUP_Td  : snprintf(name,max,"Td"); break;
        case POINT_GROUP_Th  : snprintf(name,max,"Th"); break;
        case POINT_GROUP_O   : snprintf(name,max,"O"); break;
        case POINT_GROUP_Oh  : snprintf(name,max,"Oh"); break;
        case POINT_GROUP_I   : snprintf(name,max,"I"); break;
        case POINT_GROUP_Ih  : snprintf(name,max,"Ih"); break;
        case POINT_GROUP_K   : snprintf(name,max,"K"); break;
        case POINT_GROUP_Kh  : snprintf(name,max,"Kh"); break;
        default :
            msymSetErrorDetails("Unknown point group when determining name");
            ret = MSYM_POINT_GROUP_ERROR;
            goto err;
    }
err:
    return ret;
}

//init point group, copy all point groups so we can free the original list later
msym_error_t createPointGroup(msym_thresholds_t *thresholds,int n, msym_point_group_type_t type, msym_symmetry_operation_t *primary, msym_symmetry_operation_t *sops, unsigned int sopsl, msym_point_group_t **rpg){
    msym_error_t ret = MSYM_SUCCESS;
    
    msym_point_group_t pg = {.n = n, .type = type, .sops = sops, .sopsl = sopsl, .primary = primary, .ct = NULL};

    if(MSYM_SUCCESS != (ret = setPointGroupOrder(&pg))) goto err;
    if(MSYM_SUCCESS != (ret = setPointGroupName(sizeof(pg.name)/sizeof(char),n,type,pg.name))) goto err;
    if(MSYM_SUCCESS != (ret = symmetrizePointGroup(&pg, rpg, thresholds))) goto err;
    if(((*rpg)->type == POINT_GROUP_Cnv && (*rpg)->n == 0) || ((*rpg)->type == POINT_GROUP_Dnh && (*rpg)->n == 0)){
        (*rpg)->perm = NULL;
    } else {
        if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations((*rpg)->sopsl,(*rpg)->sops, thresholds, &(*rpg)->perm))) goto err;
    }
    
err:
    return ret;
    
}


msym_error_t setPointGroupOrder(msym_point_group_t *pg){
    
    msym_error_t ret = MSYM_SUCCESS;
    switch (pg->type) {
        case (POINT_GROUP_Cs)  :
        case (POINT_GROUP_Ci)  : pg->order = 2; break;
        case (POINT_GROUP_Cn)  :
        case (POINT_GROUP_S2n) : pg->order = pg->n; break;
        case (POINT_GROUP_Cnh) :
        case (POINT_GROUP_Dn)  : pg->order = 2*pg->n; break;
        case (POINT_GROUP_Cnv) : pg->order = (pg->n == 0 ? 2 : 2*pg->n); break; //These will get a little tricky to symmetrize
        case (POINT_GROUP_Dnh) : pg->order = (pg->n == 0 ? 4 : 4*pg->n); break;
        case (POINT_GROUP_Dnd) : pg->order = 4*pg->n; break;
        case (POINT_GROUP_T)   : pg->order = 12; break;
        case (POINT_GROUP_Td)  :
        case (POINT_GROUP_Th)  :
        case (POINT_GROUP_O)   : pg->order = 24; break;
        case (POINT_GROUP_Oh)  : pg->order = 48; break;
        case (POINT_GROUP_I)   : pg->order = 60; break;
        case (POINT_GROUP_Ih)  : pg->order = 120; break;
        case (POINT_GROUP_K)   :
        case (POINT_GROUP_Kh)  : pg->order = 0;
        default                :
            msymSetErrorDetails("Point group has no primary axis for reorientation");
            goto err;
    }    
err:
    return ret;
}

msym_error_t findPointGroup(int sopsl, msym_symmetry_operation_t *sops, msym_thresholds_t *thresholds, msym_point_group_t **pg){
    msym_error_t ret = MSYM_SUCCESS;
    msymSetErrorDetails("");
    int inversion = 0, sigma = 0, nC[6] = {0,0,0,0,0,0}, linear = 0;
    msym_symmetry_operation_t *primary = NULL;
    msym_symmetry_operation_t *s = NULL;
    *pg = NULL;
    
    for(int i = 0;i < sopsl;i++){
        if(sops[i].type == PROPER_ROTATION && sops[i].order == 0){
            linear = 1;
            break;
        } else if (sops[i].type == PROPER_ROTATION && sops[i].order > 2){
            break;
        }
    }
    if(!linear){
        for(int i = 0;i < sopsl;i++){
            switch(sops[i].type){
                case PROPER_ROTATION :
                    if(primary == NULL || sops[i].order > primary->order) primary = &(sops[i]);
                    if(sops[i].order <= 5) nC[sops[i].order]++;
                    break;
                case INVERSION :
                    inversion = 1;
                    break;
                case REFLECTION :
                    sigma = 1;
                    break;
                case IMPROPER_ROTATION :
                    if(s == NULL || sops[i].order > s->order) s = &(sops[i]);
                    break;
                default :
                    break;

            }
        }
        if(nC[3] >= 2) {
            if(nC[5] >= 2){
                if(inversion){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Ih, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_I, primary, sops, sopsl, pg);
                }
            } else if (nC[4] >= 2) {
                if(inversion){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Oh, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_O, primary, sops, sopsl, pg);
                }
                
            } else if (sigma){
                if(inversion){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Th, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Td, primary, sops, sopsl, pg);
                }
            } else {
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_T, primary, sops, sopsl, pg);
            }
            
        } else if (primary == NULL){
            if(sigma){
                ret = createPointGroup(thresholds, 1, POINT_GROUP_Cs, primary, sops, sopsl, pg);
            } else if(inversion){
                ret = createPointGroup(thresholds, 1, POINT_GROUP_Ci, primary, sops, sopsl, pg);
            } else {
                ret = createPointGroup(thresholds, 1, POINT_GROUP_Cn, primary, NULL, 0, pg);
            }
        } else {
            int nC2 = 0;
            int sigma_h = 0;
            int nsigma_v = 0;
            
            //Special case for D2d
            if(primary->order == 2 && s != NULL && !vparallel(primary->v, s->v, thresholds->angle)){
                for(int i = 0; i < sopsl;i++){
                    if(sops[i].type == PROPER_ROTATION && sops[i].order == 2 && vparallel(sops[i].v, s->v, thresholds->angle)){
                        primary = &sops[i];
                        break;
                    }
                }
            }
            
            for(int i = 0; i < sopsl;i++){
                nC2 += sops[i].type == PROPER_ROTATION && sops[i].order == 2 && vperpendicular(primary->v,sops[i].v, thresholds->angle);
                sigma_h = sigma_h || (sops[i].type == REFLECTION && vparallel(primary->v,sops[i].v,thresholds->angle));
                nsigma_v += (sops[i].type == REFLECTION && vperpendicular(primary->v,sops[i].v,thresholds->angle));
            }            
            if(nC2){ //actually nC2 == primary->order but less is acceptable here since we can generate the rest
                if(sigma_h){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dnh, primary, sops, sopsl, pg);
                } else if (nsigma_v){ //actually nsigma_v == primary->order but less is acceptable here since we can generate the rest
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dnd, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dn, primary, sops, sopsl, pg);
                }
                
            } else if (sigma_h) {
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cnh, primary, sops, sopsl, pg);
            } else if(nsigma_v) { //actually nsigma_v == primary->order but less is acceptable here since we can generate the rest
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cnv, primary, sops, sopsl, pg);
            } else if(s != NULL){
                ret = createPointGroup(thresholds, s->order, POINT_GROUP_S2n, primary, sops, sopsl, pg);
            } else {
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cn, primary, sops, sopsl, pg);
            }
        }
    } else {
        for(int i = 0; i < sopsl;i++){
            inversion = inversion || sops[i].type == INVERSION;
            
            if(sops[i].type == PROPER_ROTATION && (sops[i].order == 0 || primary == NULL)){ //TODO The primary == NULL is needed because of a bug!
                primary = &(sops[i]);
            }
        }
        
        if(inversion){
            ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dnh, primary, sops, sopsl, pg);
        } else {
            ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cnv, primary, sops, sopsl, pg);
        }
    }
    
    return ret;
    
err:
    return ret;

}

msym_error_t findSubgroup(msym_subgroup_t *subgroup, msym_thresholds_t *thresholds){
    msym_error_t ret = MSYM_SUCCESS;
    int inversion = 0, sigma = 0, nC[6] = {0,0,0,0,0,0}, linear = 0;
    msym_symmetry_operation_t *primary = NULL;
    msym_symmetry_operation_t *s = NULL, *sop = NULL;;
    
    
    for(int i = 0;i < subgroup->sopsl;i++){
        if(subgroup->sops[i]->type == PROPER_ROTATION && subgroup->sops[i]->order == 0){
            linear = 1;
            break;
        } else if (subgroup->sops[i]->type == PROPER_ROTATION && subgroup->sops[i]->order > 2){
            break;
        }
    }
    if(!linear){
        for(int i = 0;i < subgroup->sopsl;i++){
            sop = subgroup->sops[i];
            if(sop->power > 1) continue;
            switch(subgroup->sops[i]->type){
                case PROPER_ROTATION :
                    if(primary == NULL || sop->order > primary->order) primary = sop;
                    if(sop->order <= 5) nC[sop->order]++;
                    break;
                case INVERSION :
                    inversion = 1;
                    break;
                case REFLECTION :
                    sigma = 1;
                    break;
                case IMPROPER_ROTATION :
                    if(s == NULL || sop->order > s->order) s = sop;
                    break;
                default :
                    break;
                    
            }
        }
        if(nC[3] >= 2) {
            if(nC[5] >= 2){
                if(inversion){
                    subgroup->type = POINT_GROUP_Ih;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_I;
                    subgroup->n = primary->order;
                }
            } else if (nC[4] >= 2) {
                if(inversion){
                    subgroup->type = POINT_GROUP_Oh;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_O;
                    subgroup->n = primary->order;
                }
                
            } else if (sigma){
                if(inversion){
                    subgroup->type = POINT_GROUP_Th;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_Td;
                    subgroup->n = primary->order;
                }
            } else {
                subgroup->type = POINT_GROUP_T;
                subgroup->n = primary->order;
            }
            
        } else if (primary == NULL){
            if(sigma){
                subgroup->type = POINT_GROUP_Cs;
                subgroup->n = 1;
            } else if(inversion){
                subgroup->type = POINT_GROUP_Ci;
                subgroup->n = 1;
            } else {
                subgroup->type = POINT_GROUP_Cn;
                subgroup->n = 1;
            }
        } else {
            int nC2 = 0;
            int sigma_h = 0;
            int nsigma_v = 0;
            
            if(primary->order == 2 && s != NULL && !vparallel(primary->v, s->v, thresholds->angle)){
                for(int i = 0; i < subgroup->sopsl;i++){
                    sop = subgroup->sops[i];
                    if(sop->power > 1) continue;
                    if(sop->type == PROPER_ROTATION && sop->order == 2 && vparallel(sop->v, s->v, thresholds->angle)){
                        primary = sop;
                        break;
                    }
                }
            }
            
            for(int i = 0; i < subgroup->sopsl;i++){
                sop = subgroup->sops[i];
                if(sop->power > 1) continue;
                nC2 += sop->type == PROPER_ROTATION && sop->order == 2 && vperpendicular(primary->v,sop->v, thresholds->angle);
                sigma_h = sigma_h || (sop->type == REFLECTION && vparallel(primary->v,sop->v,thresholds->angle));
                nsigma_v += (sop->type == REFLECTION && vperpendicular(primary->v,sop->v,thresholds->angle));
            }
            if(nC2 == primary->order){
                if(sigma_h){
                    subgroup->type = POINT_GROUP_Dnh;
                    subgroup->n = primary->order;
                } else if (nsigma_v == primary->order){
                    subgroup->type = POINT_GROUP_Dnd;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_Dn;
                    subgroup->n = primary->order;
                }
                
            } else if (sigma_h) {
                subgroup->type = POINT_GROUP_Cnh;
                subgroup->n = primary->order;
            } else if(nsigma_v == primary->order) {
                subgroup->type = POINT_GROUP_Cnv;
                subgroup->n = primary->order;
            } else if(s != NULL){
                subgroup->type = POINT_GROUP_S2n;
                subgroup->n = s->order;
            } else {
                subgroup->type = POINT_GROUP_Cn;
                subgroup->n = primary->order;
            }
        }
    } else {
        for(int i = 0; i < subgroup->sopsl;i++){
            inversion = inversion || subgroup->sops[i]->type == INVERSION;
            
            if(subgroup->sops[i]->type == PROPER_ROTATION && (subgroup->sops[i]->order == 0 || primary == NULL)){                 primary = subgroup->sops[i];
            }
        }
        
        if(inversion){
            subgroup->type = POINT_GROUP_Dnh;
            subgroup->n = primary->order;
        } else {
            subgroup->type = POINT_GROUP_Cnv;
            subgroup->n = primary->order;
        }
    }
    
    subgroup->primary = primary;
    if(MSYM_SUCCESS != (ret = setPointGroupName(sizeof(subgroup->name)/sizeof(char),subgroup->n,subgroup->type,subgroup->name))) goto err;
    return ret;
    
    
err:
    return ret;
    
}


msym_error_t transformAxes(msym_point_group_t *pg, msym_thresholds_t *thresholds){
    msym_error_t ret = MSYM_SUCCESS;
    switch (pg->type){
        case (POINT_GROUP_Ci)  :
        case (POINT_GROUP_K)   :
        case (POINT_GROUP_Kh)  :
            break;
        case (POINT_GROUP_Cs)  :
            pg->primary = pg->sops; //Transform sigma normal to z
        case (POINT_GROUP_Cn)  :
        case (POINT_GROUP_Cnh) :
        case (POINT_GROUP_S2n) :
            if(MSYM_SUCCESS != (ret = reorientAxes(pg,thresholds))) goto err;
            if(MSYM_SUCCESS != (ret = transformPrimary(pg,thresholds))) goto err;
            break;
        case (POINT_GROUP_Cnv) :
        case (POINT_GROUP_Dnh) :
            if(MSYM_SUCCESS != (ret = reorientAxes(pg,thresholds))) goto err;
            if(MSYM_SUCCESS != (ret = transformPrimary(pg,thresholds))) goto err;
            if(pg->n > 0){
                if(MSYM_SUCCESS != (ret = transformSecondary(pg,thresholds))) goto err;
            }
            break;
        case (POINT_GROUP_Dn)  :
        case (POINT_GROUP_Dnd) :
        case (POINT_GROUP_O)   :
        case (POINT_GROUP_Oh)  :<commit_msg>Look for Cs for subgroup selection<commit_after>//
//  point_group.c
//  Symmetry
//
//  Created by <NAME> on 14/04/14.
//  Copyright (c) 2014 <NAME>. 
//
//  Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
//


#include <math.h>
#include <float.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "symmetry.h"
#include "linalg.h"
#include "point_group.h"
#include "permutation.h"

#define PHI 1.618033988749894848204586834
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288419716939937510582
#endif

#define CLASSIFICATION_THRESHOLD 0.01


msym_error_t symmetrizePointGroup(msym_point_group_t *ipg, msym_point_group_t **opg, msym_thresholds_t *thresholds);
msym_error_t copyPointGroup(msym_point_group_t *ipg, msym_point_group_t *opg);
msym_error_t reorientAxes(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t transformPrimary(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t transformSecondary(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t setPointGroupOrder(msym_point_group_t *pg);
msym_error_t setPointGroupName(size_t max, int n, msym_point_group_type_t type, char *name);

msym_error_t findSecondaryAxisSigma(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC2(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC3(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC4(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);
msym_error_t findSecondaryAxisC2C5(msym_point_group_t *pg, double r[], msym_thresholds_t *thresholds);


msym_error_t generateSymmetryOperations(msym_point_group_t *pg, msym_thresholds_t *thresholds);

msym_error_t generateSymmetryOperationsImpliedS(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t generateSymmetryOperationsImpliedCPow(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t generateSymmetryOperationsImpliedSPow(msym_point_group_t *pg, msym_thresholds_t *thresholds);
msym_error_t generateSymmetryOperationsImpliedRot(msym_point_group_t *pg, msym_thresholds_t *thresholds);

void generateSymmetryOperationsCi(msym_point_group_t *pg);
void generateSymmetryOperationsCs(msym_point_group_t *pg);
void generateSymmetryOperationsCn(msym_point_group_t *pg);
void generateSymmetryOperationsCnh(msym_point_group_t *pg);
void generateSymmetryOperationsCnv(msym_point_group_t *pg);
void generateSymmetryOperationsDn(msym_point_group_t *pg);
void generateSymmetryOperationsDnh(msym_point_group_t *pg);
void generateSymmetryOperationsDnd(msym_point_group_t *pg);
void generateSymmetryOperationsS2n(msym_point_group_t *pg);
void generateSymmetryOperationsT(msym_point_group_t *pg);
void generateSymmetryOperationsTd(msym_point_group_t *pg);
void generateSymmetryOperationsTh(msym_point_group_t *pg);
void generateSymmetryOperationsO(msym_point_group_t *pg);
void generateSymmetryOperationsOh(msym_point_group_t *pg);
void generateSymmetryOperationsI(msym_point_group_t *pg);
void generateSymmetryOperationsIh(msym_point_group_t *pg);

void addSymmetryOperationsCi(msym_point_group_t *pg);
void addSymmetryOperationsCs(msym_point_group_t *pg);
void addSymmetryOperationsCn(msym_point_group_t *pg);
void addSymmetryOperationsCnh(msym_point_group_t *pg);
void addSymmetryOperationsCnv(msym_point_group_t *pg);
void addSymmetryOperationsDn(msym_point_group_t *pg);
void addSymmetryOperationsDnh(msym_point_group_t *pg);
void addSymmetryOperationsDnd(msym_point_group_t *pg);
void addSymmetryOperationsS2n(msym_point_group_t *pg);
void addSymmetryOperationsT(msym_point_group_t *pg);
void addSymmetryOperationsTd(msym_point_group_t *pg);
void addSymmetryOperationsTh(msym_point_group_t *pg);
void addSymmetryOperationsO(msym_point_group_t *pg);
void addSymmetryOperationsOh(msym_point_group_t *pg);
void addSymmetryOperationsI(msym_point_group_t *pg);
void addSymmetryOperationsIh(msym_point_group_t *pg);

msym_error_t pointGroupFromName(char *name, msym_point_group_t *pg);

int classifySymmetryOperations(msym_point_group_t *pg);
void sortSymmetryOperations(msym_point_group_t *pg, int classes);

void print_transform(double M[3][3], double axis[3]);



msym_error_t generatePointGroup(char *name, msym_thresholds_t *thresholds, msym_point_group_t **opg){
    msym_error_t ret = MSYM_SUCCESS;
    msym_point_group_t *pg = calloc(1,sizeof(msym_point_group_t));
    if(MSYM_SUCCESS != (ret = pointGroupFromName(name,pg))) goto err;
    if(MSYM_SUCCESS != (ret = generateSymmetryOperations(pg,thresholds))) goto err;
    int classes = classifySymmetryOperations(pg);
    sortSymmetryOperations(pg,classes);
    
    if((pg->type == POINT_GROUP_Cnv && pg->n == 0) || (pg->type == POINT_GROUP_Dnh && pg->n == 0)){
        pg->perm = NULL;
    } else {
        if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations(pg->sopsl,pg->sops, thresholds, &pg->perm))) goto err;
    }
    for(msym_symmetry_operation_t *s = pg->sops;s < (pg->sops + pg->sopsl);s++){
        if(pg->primary == NULL || (s->type == PROPER_ROTATION && s->order > pg->primary->order)) pg->primary = s;
    }
    
    mleye(3,pg->transform);
    
    *opg = pg;
    return ret;
    
err:
    *opg = NULL;
    free(pg);
    return ret;
}


msym_error_t pointGroupFromName(char *name, msym_point_group_t *pg){
    msym_error_t ret = MSYM_SUCCESS;
    int n = 0, gi = 0, ri = 0;;
    char g = 0, r = 0;
    
    int map[7][6];
    
    
    struct _pg_map {
        int i;
        msym_point_group_type_t type;
    } pg_map[] = {
        {1,  POINT_GROUP_Cn},
        {2,  POINT_GROUP_Cnv},
        {3,  POINT_GROUP_Cnh},
        {4,  POINT_GROUP_Ci},
        {5,  POINT_GROUP_Cs},
        {6,  POINT_GROUP_Dn},
        {7,  POINT_GROUP_Dnh},
        {8,  POINT_GROUP_Dnd},
        {9,  POINT_GROUP_S2n},
        {10, POINT_GROUP_T},
        {11, POINT_GROUP_Th},
        {12, POINT_GROUP_Td},
        {13, POINT_GROUP_O},
        {14, POINT_GROUP_Oh},
        {15, POINT_GROUP_I},
        {16, POINT_GROUP_Ih},
        {17, POINT_GROUP_K},
        {18, POINT_GROUP_Kh}
        
    };
    
    memset(map,0,sizeof(int[7][6]));
    
    map[0][0] =  1;
    map[0][1] =  2;
    map[0][2] =  3;
    map[0][4] =  4;
    map[0][5] =  5;
    map[1][0] =  6;
    map[1][2] =  7;
    map[1][3] =  8;
    map[2][0] =  9;
    map[3][0] =  10;
    map[3][2] =  11;
    map[3][3] =  12;
    map[4][0] =  13;
    map[4][2] =  14;
    map[5][0] =  15;
    map[5][2] =  16;
    map[6][0] =  17;
    map[6][2] =  18;
    
    if(sscanf(name,"%c%d%c",&g,&n,&r) < 2 && sscanf(name,"%c%c",&g,&r) < 1){
        ret = MSYM_INVALID_POINT_GROUP;
        msymSetErrorDetails("Invalid point group name %s",name);
        goto err;
    }
    
    if(n < 0) {
        ret = MSYM_INVALID_POINT_GROUP;
        msymSetErrorDetails("Invalid point group order %d",n);
        goto err;
    }
    
    switch(g){
        case 'C' : gi = 0; break;
        case 'D' : gi = 1; break;
        case 'S' : {
            if(n < 4 || n % 2 != 0){
                ret = MSYM_INVALID_POINT_GROUP;
                msymSetErrorDetails("Improper rotation order (%d) must be even",n);
                goto err;
            }
            gi = 2;
            break;
        }
        case 'T' : gi = 3; break;
        case 'O' : gi = 4; break;
        case 'I' : gi = 5; break;
        case 'K' : gi = 6; break;
        default :
            ret = MSYM_INVALID_POINT_GROUP;
            msymSetErrorDetails("Invalid point group type %c",g);
            goto err;
    }
    
    switch(r){
        case 0   : ri = 0; break;
        case 'v' : ri = 1; break;
        case 'h' : ri = 2; break;
        case 'd' : ri = 3; break;
        case 'i' : ri = 4; break;
        case 's' : ri = 5; break;
        default :
            ret = MSYM_INVALID_POINT_GROUP;
            msymSetErrorDetails("Invalid point group subtype %c",r);
            goto err;
    }
    
    int fi, fil = sizeof(pg_map)/sizeof(pg_map[0]);
    for(fi = 0;fi < fil;fi++){
        if(pg_map[fi].i == map[gi][ri]) break;
    }
    
    if(fi == fil){
        ret = MSYM_INVALID_POINT_GROUP;
        msymSetErrorDetails("Cannot find point group %s",name);
        goto err;
    }
    
    pg->type = pg_map[fi].type;
    pg->n = n;
    if(MSYM_SUCCESS != (ret = setPointGroupOrder(pg))) goto err;
    if(MSYM_SUCCESS != (ret = setPointGroupName(sizeof(pg->name)/sizeof(char),pg->n,pg->type,pg->name))) goto err;
    
err:
    return ret;
    
}


msym_error_t setPointGroupName(size_t max, int n, msym_point_group_type_t type, char *name){
    msym_error_t ret = MSYM_SUCCESS;
    switch(type) {
        case POINT_GROUP_Ci  : snprintf(name,max,"Ci"); break;
        case POINT_GROUP_Cs  : snprintf(name,max,"Cs"); break;
        case POINT_GROUP_Cn  : snprintf(name,max,"C%d",n); break;
        case POINT_GROUP_Cnh : snprintf(name,max,"C%dh",n); break;
        case POINT_GROUP_Cnv : snprintf(name,max,"C%dv",n); break;
        case POINT_GROUP_Dn  : snprintf(name,max,"D%d",n); break;
        case POINT_GROUP_Dnh : snprintf(name,max,"D%dh",n); break;
        case POINT_GROUP_Dnd : snprintf(name,max,"D%dd",n); break;
        case POINT_GROUP_S2n : snprintf(name,max,"S%d",n); break;
        case POINT_GROUP_T   : snprintf(name,max,"T"); break;
        case POINT_GROUP_Td  : snprintf(name,max,"Td"); break;
        case POINT_GROUP_Th  : snprintf(name,max,"Th"); break;
        case POINT_GROUP_O   : snprintf(name,max,"O"); break;
        case POINT_GROUP_Oh  : snprintf(name,max,"Oh"); break;
        case POINT_GROUP_I   : snprintf(name,max,"I"); break;
        case POINT_GROUP_Ih  : snprintf(name,max,"Ih"); break;
        case POINT_GROUP_K   : snprintf(name,max,"K"); break;
        case POINT_GROUP_Kh  : snprintf(name,max,"Kh"); break;
        default :
            msymSetErrorDetails("Unknown point group when determining name");
            ret = MSYM_POINT_GROUP_ERROR;
            goto err;
    }
err:
    return ret;
}

//init point group, copy all point groups so we can free the original list later
msym_error_t createPointGroup(msym_thresholds_t *thresholds,int n, msym_point_group_type_t type, msym_symmetry_operation_t *primary, msym_symmetry_operation_t *sops, unsigned int sopsl, msym_point_group_t **rpg){
    msym_error_t ret = MSYM_SUCCESS;
    
    msym_point_group_t pg = {.n = n, .type = type, .sops = sops, .sopsl = sopsl, .primary = primary, .ct = NULL};

    if(MSYM_SUCCESS != (ret = setPointGroupOrder(&pg))) goto err;
    if(MSYM_SUCCESS != (ret = setPointGroupName(sizeof(pg.name)/sizeof(char),n,type,pg.name))) goto err;
    if(MSYM_SUCCESS != (ret = symmetrizePointGroup(&pg, rpg, thresholds))) goto err;
    if(((*rpg)->type == POINT_GROUP_Cnv && (*rpg)->n == 0) || ((*rpg)->type == POINT_GROUP_Dnh && (*rpg)->n == 0)){
        (*rpg)->perm = NULL;
    } else {
        if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations((*rpg)->sopsl,(*rpg)->sops, thresholds, &(*rpg)->perm))) goto err;
    }
    
err:
    return ret;
    
}


msym_error_t setPointGroupOrder(msym_point_group_t *pg){
    
    msym_error_t ret = MSYM_SUCCESS;
    switch (pg->type) {
        case (POINT_GROUP_Cs)  :
        case (POINT_GROUP_Ci)  : pg->order = 2; break;
        case (POINT_GROUP_Cn)  :
        case (POINT_GROUP_S2n) : pg->order = pg->n; break;
        case (POINT_GROUP_Cnh) :
        case (POINT_GROUP_Dn)  : pg->order = 2*pg->n; break;
        case (POINT_GROUP_Cnv) : pg->order = (pg->n == 0 ? 2 : 2*pg->n); break; //These will get a little tricky to symmetrize
        case (POINT_GROUP_Dnh) : pg->order = (pg->n == 0 ? 4 : 4*pg->n); break;
        case (POINT_GROUP_Dnd) : pg->order = 4*pg->n; break;
        case (POINT_GROUP_T)   : pg->order = 12; break;
        case (POINT_GROUP_Td)  :
        case (POINT_GROUP_Th)  :
        case (POINT_GROUP_O)   : pg->order = 24; break;
        case (POINT_GROUP_Oh)  : pg->order = 48; break;
        case (POINT_GROUP_I)   : pg->order = 60; break;
        case (POINT_GROUP_Ih)  : pg->order = 120; break;
        case (POINT_GROUP_K)   :
        case (POINT_GROUP_Kh)  : pg->order = 0;
        default                :
            msymSetErrorDetails("Point group has no primary axis for reorientation");
            goto err;
    }    
err:
    return ret;
}

msym_error_t findPointGroup(int sopsl, msym_symmetry_operation_t *sops, msym_thresholds_t *thresholds, msym_point_group_t **pg){
    msym_error_t ret = MSYM_SUCCESS;
    msymSetErrorDetails("");
    int inversion = 0, sigma = 0, nC[6] = {0,0,0,0,0,0}, linear = 0;
    msym_symmetry_operation_t *primary = NULL;
    msym_symmetry_operation_t *s = NULL;
    *pg = NULL;
    
    for(int i = 0;i < sopsl;i++){
        if(sops[i].type == PROPER_ROTATION && sops[i].order == 0){
            linear = 1;
            break;
        } else if (sops[i].type == PROPER_ROTATION && sops[i].order > 2){
            break;
        }
    }
    if(!linear){
        for(int i = 0;i < sopsl;i++){
            switch(sops[i].type){
                case PROPER_ROTATION :
                    if(primary == NULL || sops[i].order > primary->order) primary = &(sops[i]);
                    if(sops[i].order <= 5) nC[sops[i].order]++;
                    break;
                case INVERSION :
                    inversion = 1;
                    break;
                case REFLECTION :
                    sigma = 1;
                    break;
                case IMPROPER_ROTATION :
                    if(s == NULL || sops[i].order > s->order) s = &(sops[i]);
                    break;
                default :
                    break;

            }
        }
        if(nC[3] >= 2) {
            if(nC[5] >= 2){
                if(inversion){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Ih, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_I, primary, sops, sopsl, pg);
                }
            } else if (nC[4] >= 2) {
                if(inversion){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Oh, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_O, primary, sops, sopsl, pg);
                }
                
            } else if (sigma){
                if(inversion){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Th, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Td, primary, sops, sopsl, pg);
                }
            } else {
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_T, primary, sops, sopsl, pg);
            }
            
        } else if (primary == NULL){
            if(sigma){
                ret = createPointGroup(thresholds, 1, POINT_GROUP_Cs, primary, sops, sopsl, pg);
            } else if(inversion){
                ret = createPointGroup(thresholds, 1, POINT_GROUP_Ci, primary, sops, sopsl, pg);
            } else {
                ret = createPointGroup(thresholds, 1, POINT_GROUP_Cn, primary, NULL, 0, pg);
            }
        } else {
            int nC2 = 0;
            int sigma_h = 0;
            int nsigma_v = 0;
            
            //Special case for D2d
            if(primary->order == 2 && s != NULL && !vparallel(primary->v, s->v, thresholds->angle)){
                for(int i = 0; i < sopsl;i++){
                    if(sops[i].type == PROPER_ROTATION && sops[i].order == 2 && vparallel(sops[i].v, s->v, thresholds->angle)){
                        primary = &sops[i];
                        break;
                    }
                }
            }
            
            for(int i = 0; i < sopsl;i++){
                nC2 += sops[i].type == PROPER_ROTATION && sops[i].order == 2 && vperpendicular(primary->v,sops[i].v, thresholds->angle);
                sigma_h = sigma_h || (sops[i].type == REFLECTION && vparallel(primary->v,sops[i].v,thresholds->angle));
                nsigma_v += (sops[i].type == REFLECTION && vperpendicular(primary->v,sops[i].v,thresholds->angle));
            }            
            if(nC2){ //actually nC2 == primary->order but less is acceptable here since we can generate the rest
                if(sigma_h){
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dnh, primary, sops, sopsl, pg);
                } else if (nsigma_v){ //actually nsigma_v == primary->order but less is acceptable here since we can generate the rest
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dnd, primary, sops, sopsl, pg);
                } else {
                    ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dn, primary, sops, sopsl, pg);
                }
                
            } else if (sigma_h) {
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cnh, primary, sops, sopsl, pg);
            } else if(nsigma_v) { //actually nsigma_v == primary->order but less is acceptable here since we can generate the rest
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cnv, primary, sops, sopsl, pg);
            } else if(s != NULL){
                ret = createPointGroup(thresholds, s->order, POINT_GROUP_S2n, primary, sops, sopsl, pg);
            } else {
                ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cn, primary, sops, sopsl, pg);
            }
        }
    } else {
        for(int i = 0; i < sopsl;i++){
            inversion = inversion || sops[i].type == INVERSION;
            
            if(sops[i].type == PROPER_ROTATION && (sops[i].order == 0 || primary == NULL)){ //TODO The primary == NULL is needed because of a bug!
                primary = &(sops[i]);
            }
        }
        
        if(inversion){
            ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Dnh, primary, sops, sopsl, pg);
        } else {
            ret = createPointGroup(thresholds, primary->order, POINT_GROUP_Cnv, primary, sops, sopsl, pg);
        }
    }
    
    return ret;
    
err:
    return ret;

}

msym_error_t findSubgroup(msym_subgroup_t *subgroup, msym_thresholds_t *thresholds){
    msym_error_t ret = MSYM_SUCCESS;
    int inversion = 0, sigma = 0, nC[6] = {0,0,0,0,0,0}, linear = 0;
    msym_symmetry_operation_t *primary = NULL;
    msym_symmetry_operation_t *s = NULL, *sop = NULL;;
    
    
    for(int i = 0;i < subgroup->sopsl;i++){
        if(subgroup->sops[i]->type == PROPER_ROTATION && subgroup->sops[i]->order == 0){
            linear = 1;
            break;
        } else if (subgroup->sops[i]->type == PROPER_ROTATION && subgroup->sops[i]->order > 2){
            break;
        }
    }
    if(!linear){
        for(int i = 0;i < subgroup->sopsl;i++){
            sop = subgroup->sops[i];
            if(sop->power > 1) continue;
            switch(subgroup->sops[i]->type){
                case PROPER_ROTATION :
                    if(primary == NULL || sop->order > primary->order) primary = sop;
                    if(sop->order <= 5) nC[sop->order]++;
                    break;
                case INVERSION :
                    inversion = 1;
                    break;
                case REFLECTION :
                    sigma = 1;
                    break;
                case IMPROPER_ROTATION :
                    if(s == NULL || sop->order > s->order) s = sop;
                    break;
                default :
                    break;
                    
            }
        }
        if(nC[3] >= 2) {
            if(nC[5] >= 2){
                if(inversion){
                    subgroup->type = POINT_GROUP_Ih;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_I;
                    subgroup->n = primary->order;
                }
            } else if (nC[4] >= 2) {
                if(inversion){
                    subgroup->type = POINT_GROUP_Oh;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_O;
                    subgroup->n = primary->order;
                }
                
            } else if (sigma){
                if(inversion){
                    subgroup->type = POINT_GROUP_Th;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_Td;
                    subgroup->n = primary->order;
                }
            } else {
                subgroup->type = POINT_GROUP_T;
                subgroup->n = primary->order;
            }
            
        } else if (primary == NULL){
            if(sigma){
                subgroup->type = POINT_GROUP_Cs;
                subgroup->n = 1;
            } else if(inversion){
                subgroup->type = POINT_GROUP_Ci;
                subgroup->n = 1;
            } else {
                subgroup->type = POINT_GROUP_Cn;
                subgroup->n = 1;
            }
        } else {
            int nC2 = 0;
            int sigma_h = 0;
            int nsigma_v = 0;
            
            if(primary->order == 2 && s != NULL && !vparallel(primary->v, s->v, thresholds->angle)){
                for(int i = 0; i < subgroup->sopsl;i++){
                    sop = subgroup->sops[i];
                    if(sop->power > 1) continue;
                    if(sop->type == PROPER_ROTATION && sop->order == 2 && vparallel(sop->v, s->v, thresholds->angle)){
                        primary = sop;
                        break;
                    }
                }
            }
            
            for(int i = 0; i < subgroup->sopsl;i++){
                sop = subgroup->sops[i];
                if(sop->power > 1) continue;
                nC2 += sop->type == PROPER_ROTATION && sop->order == 2 && vperpendicular(primary->v,sop->v, thresholds->angle);
                sigma_h = sigma_h || (sop->type == REFLECTION && vparallel(primary->v,sop->v,thresholds->angle));
                nsigma_v += (sop->type == REFLECTION && vperpendicular(primary->v,sop->v,thresholds->angle));
            }
            if(nC2 == primary->order){
                if(sigma_h){
                    subgroup->type = POINT_GROUP_Dnh;
                    subgroup->n = primary->order;
                } else if (nsigma_v == primary->order){
                    subgroup->type = POINT_GROUP_Dnd;
                    subgroup->n = primary->order;
                } else {
                    subgroup->type = POINT_GROUP_Dn;
                    subgroup->n = primary->order;
                }
                
            } else if (sigma_h) {
                subgroup->type = POINT_GROUP_Cnh;
                subgroup->n = primary->order;
            } else if(nsigma_v == primary->order) {
                subgroup->type = POINT_GROUP_Cnv;
                subgroup->n = primary->order;
            } else if(s != NULL){
                subgroup->type = POINT_GROUP_S2n;
                subgroup->n = s->order;
            } else {
                subgroup->type = POINT_GROUP_Cn;
                subgroup->n = primary->order;
            }
        }
    } else {
        for(int i = 0; i < subgroup->sopsl;i++){
            inversion = inversion || subgroup->sops[i]->type == INVERSION;
            
            if(subgroup->sops[i]->type == PROPER_ROTATION && (subgroup->sops[i]->order == 0 || primary == NULL)){                 primary = subgroup->sops[i];
            }
        }
        
        if(inversion){
            subgroup->type = POINT_GROUP_Dnh;
            subgroup->n = primary->order;
        } else {
            subgroup->type = POINT_GROUP_Cnv;
            subgroup->n = primary->order;
        }
    }
    
    subgroup->primary = primary;
    if(MSYM_SUCCESS != (ret = setPointGroupName(sizeof(subgroup->name)/sizeof(char),subgroup->n,subgroup->type,subgroup->name))) goto err;
    return ret;
    
    
err:
    return ret;
    
}


msym_error_t transformAxes(msym_point_group_t *pg, msym_thresholds_t *thresholds){
    msym_error_t ret = MSYM_SUCCESS;
    switch (pg->type){
        case (POINT_GROUP_Ci)  :
        case (POINT_GROUP_K)   :
        case (POINT_GROUP_Kh)  :
            break;
        case (POINT_GROUP_Cs)  :
            for(pg->primary = pg->sops; pg->primary < (pg->primary + pg->sopsl) && pg->primary->type != REFLECTION; pg->primary++){};
        case (POINT_GROUP_Cn)  :
        case (POINT_GROUP_Cnh) :
        case (POINT_GROUP_S2n) :
            if(MSYM_SUCCESS != (ret = reorientAxes(pg,thresholds))) goto err;
            if(MSYM_SUCCESS != (ret = transformPrimary(pg,thresholds))) goto err;
            break;
        case (POINT_GROUP_Cnv) :
        case (POINT_GROUP_Dnh) :
            if(MSYM_SUCCESS != (ret = reorientAxes(pg,thresholds))) goto err;
            if(MSYM_SUCCESS != (ret = transformPrimary(pg,thresholds))) goto err;
            if(pg->n > 0){
                if(MSYM_SUCCESS != (ret = transformSecondary(pg,thresholds))) goto err;
            }
            break;
        case (POINT_GROUP_Dn)  :
        case (POINT_GROUP_Dnd) :
        case (POINT_GROUP_O)   :
        case (POINT_GROUP_Oh)  :
{'fraction_non_alphanumeric': 0.12583678517442431, 'fraction_numerical': 0.0136241386756709, 'mean_word_length': 2.145803742357809, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 471, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6493422', 'n_tokens_mistral': 17939, 'n_tokens_neox': 17035, 'n_words': 3372}