package ciss.phase_viewer.atomcoord.pmodel;

final class PointGroupOperation {
    enum PointGroup {
        Oh, D6h
    };

    PointGroup pointGroup;
    private int ng0; // Oh: 48, D6h: 24
    int opMat[][][];
    int table[][];
    int inv[];

    public static void main(String args[]) {
        int i, j, n;

        PointGroupOperation pgo = new PointGroupOperation(PointGroup.Oh); // 'Oh'
                                                                          // or
                                                                          // 'D6h'

        // rotation matrix
        for (n = 0; n < pgo.ng0; ++n) {
            System.out.printf("%3d\n", n);
            for (i = 0; i < 3; ++i) {
                for (j = 0; j < 3; ++j)
                    System.out.printf("%3d", pgo.opMat[n][i][j]);
                System.out.printf("\n");
            }
            System.out.printf("\n");
        }
        // multiplication table
        for (i = 0; i < pgo.ng0; ++i) {
            for (j = 0; j < pgo.ng0; ++j)
                System.out.printf("%3d", pgo.table[i][j] + 1);
            System.out.printf("\n");
        }
        System.out.printf("\n");
        // inverse element
        for (n = 0; n < pgo.ng0; ++n)
            System.out.printf("%3d", pgo.inv[n] + 1);
        System.out.printf("\n");
    }

    PointGroupOperation(PointGroup pg) {
        pointGroup = pg;
        switch (pointGroup) {
        case Oh:
            ng0 = 48;
            opMat = new int[ng0][3][3];
            table = new int[ng0][ng0];
            inv = new int[ng0];
            this.operationMatricesOh();
            break;
        case D6h:
            ng0 = 24;
            opMat = new int[ng0][3][3];
            table = new int[ng0][ng0];
            inv = new int[ng0];
            this.operationMatricesD6h();
            break;
        default:
            // BUG
            PmodelRuntimeException ex = new PmodelRuntimeException();
            throw ex;
        }
        this.makeMultiplicationTable();
    }

    int getNg0() {
        return this.ng0;
    }

    static PointGroup getPointGroup(int n) {
        PointGroup pg_;
        if (n >= 143 && n <= 194) {
            pg_ = PointGroup.D6h;
        } else {
            pg_ = PointGroup.Oh;
        }
        return pg_;
    }

    private void operationMatricesOh() {
        int cosPhi, sinPhi, cosTheta, sinTheta, cosPsi, sinPsi;
        int i, j, n;

        int eulerOh[][] = { // [24][3]
        { 0, 0, 0 }, { 0, 2, 2 }, { 0, 2, 0 }, { 0, 0, 2 }, { 0, 1, 1 },
                { 2, 1, 3 }, { 2, 1, 1 }, { 0, 1, 3 }, { 1, 1, 2 },
                { 3, 1, 0 }, { 1, 1, 0 }, { 3, 1, 2 }, { 0, 2, 1 },
                { 0, 2, 3 }, { 0, 1, 2 }, { 1, 1, 1 }, { 2, 1, 0 },
                { 3, 1, 3 }, { 3, 1, 1 }, { 0, 1, 0 }, { 1, 0, 0 },
                { 1, 1, 3 }, { 2, 1, 2 }, { 3, 0, 0 } };

        for (n = 0; n < 24; ++n) {
            // from RIKAGAKU-JITEN
            cosPhi = cos4Oh(eulerOh[n][0]);
            sinPhi = sin4Oh(eulerOh[n][0]);
            cosTheta = cos4Oh(eulerOh[n][1]);
            sinTheta = sin4Oh(eulerOh[n][1]);
            cosPsi = cos4Oh(eulerOh[n][2]);
            sinPsi = sin4Oh(eulerOh[n][2]);

            opMat[n][0][0] = cosPsi * cosPhi * cosTheta - sinPsi * sinPhi;
            opMat[n][0][1] = -sinPsi * cosPhi * cosTheta - cosPsi * sinPhi;
            opMat[n][0][2] = cosPhi * sinTheta;
            opMat[n][1][0] = cosPsi * sinPhi * cosTheta + sinPsi * cosPhi;
            opMat[n][1][1] = -sinPsi * sinPhi * cosTheta + cosPsi * cosPhi;
            opMat[n][1][2] = sinPhi * sinTheta;
            opMat[n][2][0] = -cosPsi * sinTheta;
            opMat[n][2][1] = sinPsi * sinTheta;
            opMat[n][2][2] = cosTheta;

            for (i = 0; i < 3; ++i)
                for (j = 0; j < 3; ++j)
                    opMat[n + 24][i][j] = -opMat[n][i][j];
        }
    }

    private int cos4Oh(int angle) { // cos(PI / 2 * angle); angle = 0, 1, 2, or
                                    // 3
        switch (angle) {
        case 0:
            return 1;
        case 2:
            return -1;
        default:
            return 0;
        }
    }

    private int sin4Oh(int angle) { // sin(PI / 2 * angle)
        switch (angle) {
        case 1:
            return 1;
        case 3:
            return -1;
        default:
            return 0;
        }
    }

    private void operationMatricesD6h() {
        int opMat1[][] = new int[3][3];
        int opMat2[][] = new int[3][3];
        int i, j, k, n;

        int eulerD6h[][] = { // [12][3]
        { 0, 0, 0 }, { 1, 0, 0 }, { 2, 0, 0 }, { 3, 0, 0 }, { 4, 0, 0 },
                { 5, 0, 0 }, { 0, 3, 0 }, { 4, 3, 0 }, { 2, 3, 0 },
                { 3, 3, 0 }, { 1, 3, 0 }, { 5, 3, 0 } };

        for (n = 0; n < 12; ++n) {
            // clear matrices
            for (i = 0; i < 3; ++i)
                for (j = 0; j < 3; ++j) {
                    opMat1[i][j] = 0;
                    opMat2[i][j] = 0;
                    opMat[n][i][j] = 0;
                }

            // the First rotation (around z-axis)
            opMat1[2][2] = 1;
            switch (eulerD6h[n][0] % 3) {
            case 0:
                opMat1[0][0] = 1;
                // opMat1[1][0] = 0;
                opMat1[1][1] = 1;
                break;
            case 1:
                opMat1[0][0] = 1;
                opMat1[1][0] = 1;
                // opMat1[1][1] = 0;
                break;
            case 2:
                // opMat1[0][0] = 0;
                opMat1[1][0] = 1;
                opMat1[1][1] = -1;
                break;
            }
            if (eulerD6h[n][0] >= 3) {
                opMat1[0][0] = -opMat1[0][0];
                opMat1[1][0] = -opMat1[1][0];
                opMat1[1][1] = -opMat1[1][1];
            }
            opMat1[0][1] = -opMat1[1][0];

            // the Second rotation (around y-axis)
            opMat2[1][1] = 1;
            switch (eulerD6h[n][1]) {
            case 0:
                opMat2[0][0] = 1;
                // opMat2[0][1] = 0;
                opMat2[2][2] = 1;
                break;
            case 3:
                opMat2[0][0] = -1;
                opMat2[0][1] = 1;
                opMat2[2][2] = -1;
                break;
            default:
                // ERROR
                PmodelRuntimeException ex = new PmodelRuntimeException();
                throw ex;
            }

            // the Third rotation (around z-axis)
            ; // nothing to do.

            // Matrix Multiplication
            for (i = 0; i < 3; ++i)
                for (j = 0; j < 3; ++j)
                    for (k = 0; k < 3; ++k)
                        opMat[n][i][j] += opMat2[i][k] * opMat1[k][j];
        }

        for (n = 0; n < 12; ++n)
            for (i = 0; i < 3; ++i)
                for (j = 0; j < 3; ++j)
                    opMat[n + 12][i][j] = -opMat[n][i][j];
    }

    private void makeMultiplicationTable() {
        int tempMat[][] = new int[3][3];
        int i, j, k, l, m;

        // Clear tables
        for (i = 0; i < ng0; ++i) {
            // inv[i] = -1;
            for (j = 0; j < ng0; ++j)
                table[i][j] = -1;
        }

        // Make multiplication table
        for (i = 0; i < ng0; ++i)
            for (j = 0; j < ng0; ++j) {
                // Matrix multiplication
                for (k = 0; k < 3; ++k)
                    for (l = 0; l < 3; ++l) {
                        tempMat[k][l] = 0;
                        for (m = 0; m < 3; ++m)
                            tempMat[k][l] += opMat[j][k][m] * opMat[i][m][l];
                    }
                // Search
                for (k = 0; k < ng0; ++k)
                    if (this.isIdentical(k, tempMat)) {
                        table[i][j] = k;
                        break;
                    }
            }

        // Make inverse element table
        for (i = 0; i < ng0; ++i) {
            for (j = 0; j < ng0; ++j)
                if (table[i][j] == 0)
                    break;
            inv[i] = j;
        }
    }

    private boolean isIdentical(int n, int mat[][]) {
        int i, j;

        if (n < 0 || n >= ng0) {
            PmodelRuntimeException ex = new PmodelRuntimeException();
            ex.message = "Illigal value of operation.";
            throw ex;
        }
        // Check
        for (i = 0; i < 3; ++i)
            for (j = 0; j < 3; ++j)
                if (mat[i][j] != opMat[n][i][j])
                    return false;
        return true;
    }

    Atom rotate(Atom atom, int ope) { // -> calss ATOM
        Atom a1 = atom.copy();
        for (int i = 0; i < 3; ++i) {
            a1.position[i] = 0.0;
            for (int j = 0; j < 3; ++j)
                a1.position[i] += opMat[ope][i][j] * atom.position[j];
        }

        return a1;
    }
}
