package ciss.phase_viewer.atomcoord.pmodel;

import java.util.ArrayList;

public class Pmodel {
    private static final boolean DEBUG = false;
    private static final double EPS = 1.0e-13;
    private static final double criticalDistance = 0.1; // in Angstrom

    private SpaceGroup sg;
    private Cif cif;

    private double[][] cellMatrix;
    private double[] ttl; // for calculation of distance between atoms (see
                          // 'calttl' in b_LDOS_f77.F)

    ArrayList<Atom> atomList;

    public Pmodel() {
        cif = new Cif();
        cellMatrix = new double[3][3];
        ttl = new double[6];
        atomList = new ArrayList<Atom>(100);
    }

    public int getIntTablesNum() {
        return cif.getIntTablesNum();
    }

    public int getNumSymmetryEquivPosInCif() {
        return cif.symmetryEquivPos.size();
    }

    public String getHermannMauguin() {
        return cif.getHermannMauguin();
    }

    /*
     * public void setCriticalDistance(double a) { criticalDistance = a; }
     * 
     * public double getCriticalDistance() { return criticalDistance; }
     */

    public void makeSpaceGroup() throws PmodelException {
        PointGroupOperation.PointGroup pg_ = getPG();
        sg = new SpaceGroup(pg_, cif.symmetryEquivPos);
    }

    public void makeSpaceGroup(int choice) throws PmodelException {
        sg = new SpaceGroup(getIntTablesNum(), choice);
    }

    /*
     * public static void main(String[] args) throws PmodelException,
     * IOException { double[] cellAngle; double[] cellLength; int multi; int i,
     * j;
     * 
     * if (args.length != 1) {
     * System.out.println("Usage: java Pmodel CIF_file_name"); System.exit(1); }
     * 
     * Pmodel unitCell = new Pmodel();
     * 
     * multi = Cif.getMultiplicity(args[0]); // <SELECT_MULTI> if (multi == 0) {
     * System.out.println("Warning : No 'data_' tag found.");
     * System.out.println("Read the CIF file forcibly."); } else if (multi == 1)
     * { ; // nothing to do } else { // SELECT
     * System.out.println("Multiplicity of CIF: " + multi); BufferedReader br =
     * new BufferedReader(new InputStreamReader(System.in)); String str =
     * br.readLine(); multi = Integer.parseInt(str);
     * System.out.println("Selected: " + multi); } // </SELECT_MULTI>
     * System.out.println("CIF file name: " + args[0]);
     * unitCell.readCIF(args[0], multi); // read CIF
     * 
     * int numOpe = unitCell.getNumSymmetryEquivPosInCif(); BufferedReader br =
     * new BufferedReader(new InputStreamReader(System.in)); String str;
     * 
     * cellLength = unitCell.cif.getCellLength(); cellAngle =
     * unitCell.cif.getCellAngle(); if (DEBUG) { System.out.println("<CIF>");
     * System.out.println("Cell System: " + unitCell.cif.getCrystalSystem());
     * System.out.println("symmetry");
     * //System.out.println(unitCell.cif.getHermannMauguin() + "\t" +
     * unitCell.getIntTablesNum());
     * 
     * System.out.println("length"); for (i = 0; i < 3; ++i)
     * System.out.println(cellLength[i]);
     * 
     * System.out.println("angle"); for (i = 0; i < 3; ++i)
     * System.out.println(cellAngle[i]);
     * 
     * System.out.println("atom(s): " + unitCell.cif.getNumberOfAtom()); for (i
     * = 0; i < unitCell.cif.getNumberOfAtom(); ++i)
     * System.out.println(unitCell.cif.atomList.get(i));
     * 
     * System.out.println("symmetry operation(s): " + numOpe); for (i = 0; i <
     * numOpe; ++i) System.out.println(unitCell.cif.symmetryEquivPos.get(i));
     * 
     * System.out.println("</CIF>"); str = br.readLine(); // WAIT }
     * 
     * System.out.println("<SpaceGroup>"); if (numOpe > 1) {
     * System.out.println("use equiv_pos_as_xyz"); unitCell.makeSpaceGroup(); }
     * else if (unitCell.getIntTablesNum() >= 0) {
     * System.out.println("use int_tables_num: " + unitCell.getIntTablesNum());
     * int choice = SpaceGroup.getMaxChoice(unitCell.getIntTablesNum()); //
     * <SELECT_CHOICE> if (choice > 1) { // SELECT
     * System.out.println("Maximum number of CHOICE: " + choice); str =
     * br.readLine(); // INPUT CHOICE choice = Integer.parseInt(str);
     * System.out.println("CHOICE: " + choice); } --choice; // </SELECT_CHOICE>
     * unitCell.makeSpaceGroup(choice); } else { PmodelException ex = new
     * PmodelException(); ex.message = "Cannot detect the space group.'"; throw
     * ex; } System.out.println("</SpaceGroup>");
     * 
     * // System.out.println("<CHECK>"); // unitCell.checkSymmetry(); //
     * System.out.println("</CHECK>");
     * 
     * str = br.readLine(); // WAIT
     * 
     * unitCell.makeLattice(); // EXPAND ATOM(S)
     * 
     * System.out.println("Space Group Number: " + unitCell.getIntTablesNum());
     * if (unitCell.sg.generator != null)
     * System.out.println("Number of Choice  : " +
     * unitCell.sg.generator.getNumOfChoice());
     * 
     * System.out.println("Cell Length:");
     * System.out.printf("a = %f, b = %f, c = %f\n", cellLength[0],
     * cellLength[1], cellLength[2]);
     * 
     * System.out.println("Cell Angle:");
     * System.out.printf("alpha = %f, beta = %f, gamma = %f\n", cellAngle[0],
     * cellAngle[1], cellAngle[2]);
     * 
     * System.out.println("Cell Vector:"); double[][] cellVec =
     * unitCell.getCellVec(); for (i = 0; i < 3; ++i) {
     * System.out.printf("%c_vector = ", 'a' + i); for (j = 0; j < 3; ++j)
     * System.out.printf("%17.12f", cellVec[i][j]); System.out.printf("\n"); }
     * 
     * int cpNumAtom; cpNumAtom = unitCell.getNumberOfAtoms();
     * System.out.println("Number of Atoms: " + cpNumAtom);
     * 
     * double[][] inPos; //, exPos; String[] cpSymbol; int totalAtom; inPos =
     * unitCell.getInternalPosition(); //exPos =
     * unitCell.getCartesianPosition(); totalAtom = unitCell.getNumberOfAtoms();
     * cpSymbol = unitCell.getSymbol();
     * 
     * System.out.println("Internal Position:"); for (i = 0; i < cpNumAtom; ++i)
     * System.out.printf("%7s%16.10f%16.10f%16.10f\n", cpSymbol[i], inPos[i][0],
     * inPos[i][1], inPos[i][2]);
     * 
     * if (DEBUG) { System.out.println("Cartesian Position:"); for (i = 0; i <
     * totalAtom; ++i) System.out.println(cpSymbol[i] + " " + exPos[i][0] + " "
     * + exPos[i][1] + " " + exPos[i][2]); } }
     */

    public void readCIF(String fileName, int ndb) throws PmodelException {
        cif.read(fileName, ndb);
    }

    public PointGroupOperation.PointGroup getPG() {
        final double eps = 0.00005;

        double[] cellAngle = cif.getCellAngle();
        double[] cellLength = cif.getCellLength();

        // by '_symmetry_cell_setting' in CIF
        switch (cif.getCrystalSystem()) {
        case TRICLINIC:
        case MONOCLINIC:
        case ORTHORHOMBIC:
        case TETRAGONAL:
        case CUBIC:
            return PointGroupOperation.PointGroup.Oh;
        case TRIGONAL:
        case RHOMBOHEDRAL:
        case HEXAGONAL:
            return PointGroupOperation.PointGroup.D6h;
        default:
            ; // undefined crystal system
        }
        // by '_symmetry_Int_Tables_number' in CIF
        int sgn = getIntTablesNum();
        if (sgn > 0)
            return PointGroupOperation.getPointGroup(sgn);

        // by cell-shape
        PointGroupOperation.PointGroup pg_ = PointGroupOperation.PointGroup.D6h;
        if (cellLength[0] != cellLength[1])
            pg_ = PointGroupOperation.PointGroup.Oh;
        if (Math.abs(cellAngle[0] - 90.0) > eps)
            pg_ = PointGroupOperation.PointGroup.Oh;
        if (Math.abs(cellAngle[1] - 90.0) > eps)
            pg_ = PointGroupOperation.PointGroup.Oh;
        if (Math.abs(cellAngle[2] - 120.0) > eps)
            pg_ = PointGroupOperation.PointGroup.Oh;
        return pg_;
    }

    /*
     * public void checkSymmetry() throws PmodelException { }
     */

    public void makeLattice() throws PmodelException {
        double[] cellAngle = cif.getCellAngle();
        double[] cellLength = cif.getCellLength();
        int[] nc = new int[1];

        int i, j, k;
        double temp;

        // make cellMatrix
        double alpha = cellAngle[0] * Math.PI / 180.0;
        double beta = cellAngle[1] * Math.PI / 180.0;
        double gamma = cellAngle[2] * Math.PI / 180.0;
        // rotation
        for (i = 0; i < 3; ++i)
            for (j = 0; j < 3; ++j)
                cellMatrix[i][j] = 0.0; // clear
        cellMatrix[0][0] = 1.;
        cellMatrix[0][1] = Math.cos(gamma);
        cellMatrix[0][2] = Math.cos(beta);
        cellMatrix[1][1] = Math.sin(gamma);
        temp = (Math.cos(alpha) - Math.cos(beta) * Math.cos(gamma))
                / Math.sin(gamma);
        cellMatrix[1][2] = temp;
        cellMatrix[2][2] = Math.sqrt(Math.sin(beta) * Math.sin(beta) - temp
                * temp);
        // magnification
        for (i = 0; i < 3; ++i)
            for (j = 0; j < 3; ++j)
                cellMatrix[i][j] *= cellLength[j];

        // make ttl
        for (i = 0; i < 6; ++i)
            ttl[i] = 0.0;
        for (i = 0; i < 3; ++i) {
            k = (i + 1) % 3;
            for (j = 0; j < 3; ++j) {
                ttl[i] += cellMatrix[j][i] * cellMatrix[j][i];
                ttl[i + 3] += cellMatrix[j][i] * cellMatrix[j][k];
            }
        }

        // Expand Atoms by Symmetry Operation
        expandAtoms();
        eliminateDuplicateAtoms();
    }

    private void expandAtoms() {
        if (DEBUG)
            System.out.println(cif.atomList.size() + "  "
                    + sg.pointGroupOperation.getNg0());
        for (int i = 0; i < cif.atomList.size(); ++i) { // atom
            for (int j = 0; j < sg.pointGroupOperation.getNg0(); ++j) { // rotation
                if (!sg.isActive(j))
                    continue;
                Atom atm = sg.pointGroupOperation
                        .rotate(cif.atomList.get(i), j);
                for (int k = 0; k < 4; ++k) { // primitive translation
                    if (sg.transVec[j][k] == null)
                        continue;
                    atomList.add(sg.transVec[j][k].shift(atm));
                }
            }
        }
        if (DEBUG)
            System.out.println("Number of Atoms (including duplications): "
                    + atomList.size());
    }

    private void eliminateDuplicateAtoms() {
        int i = 0, j, count = 0;

        while (i < atomList.size() - 1) {
            j = i + 1;
            while (j < atomList.size()) {
                if (isDuplicate(i, j)) {
                    atomList.remove(j);
                    ++count;
                } else
                    ++j;
            }
            ++i;
        }
        if (DEBUG)
            System.out.println("Number of Eliminated atom(s): " + count);

        return;
    }

    private boolean isDuplicate(int m, int n) {
        // check atom symbol
        ;
        // check distance
        return (distanceSq(m, n) < criticalDistance * criticalDistance);
    }

    private double distanceSq(int m, int n) {
        double[] diff = new double[3];
        double a, b, temp, dis;
        int i, j;

        for (i = 0; i < 3; ++i) {
            a = atomList.get(m).position[i]; // internal position
            b = atomList.get(n).position[i]; // internal position
            temp = a - b;
            while (temp < -0.5)
                temp += 1.0;
            while (temp > 0.5)
                temp -= 1.0;
            diff[i] = temp;
        }

        dis = 0.0;
        for (i = 0; i < 3; ++i) {
            j = (i + 1) % 3;
            dis += ttl[i] * diff[i] * diff[i] + 2.0 * ttl[i + 3] * diff[i]
                    * diff[j];
        }

        return dis;
    }

    public double[] getCellAngle() {
        return cif.getCellAngle();
    }

    public double[] getCellLength() {
        return cif.getCellLength();
    }

    public int getNumberOfAtoms() {
        return atomList.size();
    }

    public String[] getSymbol() {

        int numAtom = getNumberOfAtoms();
        String[] symbol = new String[numAtom];

        for (int i = 0; i < numAtom; ++i)
            symbol[i] = atomList.get(i).symbol;

        return symbol;
    }

    public String[] getAlias() {

        int numAtom = getNumberOfAtoms();
        String[] alias = new String[numAtom];

        for (int i = 0; i < numAtom; ++i)
            alias[i] = atomList.get(i).alias;

        return alias;
    }

    public double[][] getInternalPosition() {

        int numAtom = getNumberOfAtoms();
        double[][] pos = new double[numAtom][3];

        for (int i = 0; i < numAtom; ++i)
            for (int j = 0; j < 3; ++j)
                pos[i][j] = atomList.get(i).position[j];

        return pos;
    }

    /*
     * private double[][] getCartesianPosition() { int i, j, n;
     * 
     * int numAtom = getNumberOfAtoms(); double[][] pos = new
     * double[numAtom][3];
     * 
     * // operate the matrix for (n = 0; n < numAtom; ++n) for (i = 0; i < 3;
     * ++i) { pos[n][i] = 0.; for (j = 0; j < 3; ++j) pos[n][i] +=
     * cellMatrix[i][j] * atomList.get(n).position[j]; }
     * 
     * for (n = 0; n < numAtom; ++n) for (i = 0; i < 3; ++i) if
     * (Math.abs(pos[n][i]) < EPS) pos[n][i] = 0.0;
     * 
     * return pos; }
     */

    public double[][] getCellVec() {

        double[][] cellVec = new double[3][3];
        for (int i = 0; i < 3; ++i)
            for (int j = 0; j < 3; ++j) {
                cellVec[i][j] = cellMatrix[j][i];
                if (Math.abs(cellVec[i][j]) < EPS)
                    cellVec[i][j] = 0.;
            }

        return cellVec;
    }
}
