/*
!=======================================================================
!
!  PROGRAM  PHASE-Viewer  (PHASE-Viewer 2014.01 ver.3.3.0)
!
!  Created on ----
!  AUTHOR(S): KOGA, Junichiro
!  File : BondCalculator.java
!  
!  Contact address :  Phase System Consortium
!                     E-mail: phase_system@nims.go.jp URL https://azuma.nims.go.jp
!
!
!   Since 2002, this program set had been intensively developed as a part of the following 
!  national projects supported by the Ministry of Education, Culture, Sports, Science and
!  Technology (MEXT) of Japan; "Frontier Simulation Software for Industrial Science
!  (FSIS)" from 2002 to 2005, "Revolutionary Simulation Software (RSS21)" from 2006 to
!  2008. "Research and Development of Innovative Simulation Software (RISS)" from 2008
!  to 2013. These projects is lead by the Center for Research on Innovative Simulation 
!  Software (CISS), the Institute of Industrial Science (IIS), the University of Tokyo.
!   Since 2013, this program set has been further developed centering on PHASE System
!  Consortium. 
!   The activity of development of this program set has been supervised by Takahisa Ohno.
!
!=======================================================================
 */
package ciss.phase_viewer.acviewer.scenegraphelements.bond;

import java.util.Arrays;
import java.util.Vector;

import org.apache.log4j.Logger;

import ciss.phase_viewer.acviewer.ConfigDataManager;
import ciss.phase_viewer.acviewer.J3DDataManager;
import ciss.phase_viewer.common.ElementDef;
import ciss.phase_viewer.common.ElementInfo;

public class BondCalculator {
    private static Logger logger = Logger.getLogger(BondCalculator.class
            .getName());

    private int NumAt;

    private int NumBondMax;

    private int NumBond;

    private double rbond;

    private J3DDataManager mACVD;

    private ConfigDataManager mCD;

    private double[][] xyz;

    public BondCalculator() {
    }

    /**
     * @param mACVD
     *            原子配置表示方法の設定が保持されているオブジェクト.
     * @param mCD
     *            原子配置データが保持されているオブジェクト.
     */
    public BondCalculator(J3DDataManager mACVD, ConfigDataManager mCD) {
        this.mACVD = mACVD;
        this.mCD = mCD;
        this.NumAt = mCD.getNumAt();
        this.NumBondMax = NumAt * (NumAt - 1) / 2;
    }

    //
    // public void calculateBond(double[][] xyz) {
    // this.xyz = xyz;
    // calcNearestNeighbor();
    // }

    /**
     * 総ボンド数を取得するメソッド. getNearestNeighborメソッドで ボンドの計算を行う最中にカウントされるので,
     * 必ずgetNearestNeighborメソッド を呼んだ後に実行する.
     * 
     * @return 総ボンド数
     */
    public int getNumBond() {
        return this.NumBond;
    }

    public int getNumBondMax() {
        return this.NumBondMax;
    }

    private double[][] Hatoms;

    private double[][][] HbondCandidates;

    private double[] HbondRmax;

    private double[][] hxyz;

    private int[] map;

    private boolean prepareForHbond() {
        map = new int[mCD.getNumAt()];
        ElementInfo einfo = ElementInfo.getElementInfo();
        ElementDef[] candidates = einfo.getHBondDonorAcceptorCandidates();
        if (candidates == null || candidates.length == 0) {
            // logger.error("no H-bond candidates.");
            return false;
        }

        String[] elements = mCD.getElements();
        Vector vec = new Vector();
        int count = 0;
        for (int i = 0; i < elements.length; i++) {
            if (elements[i].equals("H")) {
                vec.addElement(xyz[i]);
                map[i] = count;
                count++;
            }
        }
        if (vec.size() == 0) {
            logger.debug("found no H atoms");
            return false;
        }
        Hatoms = new double[vec.size()][];
        logger.debug("num. H atoms: " + Hatoms.length);
        vec.copyInto(Hatoms);

        vec.clear();
        HbondCandidates = new double[candidates.length][][];
        HbondRmax = new double[candidates.length];
        boolean found = false;
        rbondMax = -10000000;
        double lenmax = mCD.getLenMax();
        for (int i = 0; i < candidates.length; i++) {
            for (int j = 0; j < elements.length; j++) {
                if (candidates[i].getSymbol().equals(elements[j])) {
                    vec.addElement(xyz[j]);
                    map[j] = count;
                    count++;
                }
            }
            if (vec.size() != 0) {
                HbondCandidates[i] = new double[vec.size()][];
                HbondRmax[i] = candidates[i].getMaxHbondDistance() / lenmax;
                if (HbondRmax[i] > rbondMax)
                    rbondMax = HbondRmax[i];
                vec.copyInto(HbondCandidates[i]);
                found = true;
            }
            vec.clear();
        }

        return found;
    }

    private void preProcessForHbond() {
        // 座標を囲むセルを決めておく
        double[] min = new double[3];
        double[] max = new double[3];
        for (int i = 0; i < 3; i++) {
            min[i] = 10000000;
            max[i] = -10000000;
        }
        int numat2consider = 0;
        for (int i = 0; i < HbondCandidates.length; i++) {
            if (HbondCandidates[i] == null)
                continue;
            for (int j = 0; j < HbondCandidates[i].length; j++) {
                for (int k = 0; k < 3; k++) {
                    if (HbondCandidates[i][j][k] > max[k]) {
                        max[k] = HbondCandidates[i][j][k];
                    }
                    if (HbondCandidates[i][j][k] < min[k]) {
                        min[k] = HbondCandidates[i][j][k];
                    }
                }
                numat2consider++;
            }
        }

        for (int i = 0; i < Hatoms.length; i++) {
            for (int j = 0; j < Hatoms[i].length; j++) {
                if (Hatoms[i][j] > max[j]) {
                    max[j] = Hatoms[i][j];
                }
                if (Hatoms[i][j] < min[j]) {
                    min[j] = Hatoms[i][j];
                }
            }
            numat2consider++;
        }

        bigCell = new double[3];
        for (int i = 0; i < 3; i++) {
            bigCell[i] = max[i] - min[i] + 0.01; // ほんの少しだけ大きくした方が便利
        }

        NumAt = numat2consider;
        xyz_buff = new double[numat2consider][3];
        hxyz = new double[numat2consider][3];
        covalentRadius = new double[numat2consider];
        int count = 0;
        for (int i = 0; i < Hatoms.length; i++) {
            for (int j = 0; j < 3; j++) {
                xyz_buff[count][j] = Hatoms[i][j] - min[j];
                hxyz[count][j] = Hatoms[i][j];
                covalentRadius[count] = 0;
            }
            count++;
        }
        for (int i = 0; i < HbondCandidates.length; i++) {
            if (HbondCandidates[i] == null)
                continue;
            for (int j = 0; j < HbondCandidates[i].length; j++) {
                for (int k = 0; k < 3; k++) {
                    xyz_buff[count][k] = HbondCandidates[i][j][k] - min[k];
                    hxyz[count][k] = HbondCandidates[i][j][k];
                    covalentRadius[count] = HbondRmax[i];
                }
                count++;
            }
        }

        // 小さなセルの大きさ, division数を決めておく.
        ndiv = new int[3];
        smallCell = new double[3];
        logger.debug("rbondMax: " + rbondMax);
        for (int i = 0; i < 3; i++) {
            ndiv[i] = (int) (bigCell[i] / rbondMax) + 1;
            smallCell[i] = bigCell[i] / ndiv[i];
        }

        logger.debug("ndiv: " + ndiv[0] + ", " + ndiv[1] + ", " + ndiv[2]);
        logger.debug("smallCell: " + smallCell[0] + ", " + smallCell[1] + ", "
                + smallCell[2]);

        // セルと原子の対応を決めておく
        cell_atom_map = new Vector[ndiv[0]][ndiv[1]][ndiv[2]];
        for (int i1 = 0; i1 < ndiv[0]; i1++) {
            for (int i2 = 0; i2 < ndiv[1]; i2++) {
                for (int i3 = 0; i3 < ndiv[2]; i3++) {
                    cell_atom_map[i1][i2][i3] = null;
                }
            }
        }

        for (int i = 0; i < xyz_buff.length; i++) {
            int i1 = (int) (xyz_buff[i][0] / smallCell[0]);
            int i2 = (int) (xyz_buff[i][1] / smallCell[1]);
            int i3 = (int) (xyz_buff[i][2] / smallCell[2]);
            if (cell_atom_map[i1][i2][i3] == null) {
                cell_atom_map[i1][i2][i3] = new Vector();
            }
            logger.debug("i1, i2, i3: " + i1 + ", " + i2 + ", " + i3);
            cell_atom_map[i1][i2][i3].addElement(new Integer(i));
        }

    }

    private double[] bigCell;

    private double[] smallCell;

    private int[] ndiv;

    private double[][] xyz_buff;

    private double[] covalentRadius;

    private double rbondMax;

    private Vector[][][] cell_atom_map;

    private String[] elements;

    private void preProcess() {
        int i;
        int j;
        this.NumAt = mCD.getNumAt();

        // 座標を囲むセルを決めておく
        double[] min = new double[3];
        double[] max = new double[3];
        for (i = 0; i < 3; i++) {
            min[i] = 10000000;
            max[i] = -10000000;
        }
        for (i = 0; i < NumAt; i++) {
            for (j = 0; j < 3; j++) {
                if (xyz[i][j] > max[j]) {
                    max[j] = xyz[i][j];
                }
                if (xyz[i][j] < min[j]) {
                    min[j] = xyz[i][j];
                }
            }
        }
        bigCell = new double[3];
        for (i = 0; i < 3; i++) {
            bigCell[i] = max[i] - min[i] + 0.01; // ほんの少しだけ大きくした方が便利
        }
        xyz_buff = new double[NumAt][3];
        for (i = 0; i < NumAt; i++) {
            for (j = 0; j < 3; j++) {
                xyz_buff[i][j] = xyz[i][j] - min[j];
            }
        }

        // "最大ボンド長"を決めておく ついでに共有半径を決めておく
        covalentRadius = new double[NumAt];
        int NumEl = mACVD.getNumEl();
        elements = mCD.getElements();

        String[] elementString = new String[NumEl];
        double[] covRad = new double[NumEl];
        double lenMax = mCD.getLenMax();
        double bondFactor = mACVD.getBondFactor();

        for (int jele = 0; jele < NumEl; jele++) {
            elementString[jele] = mACVD.getElement()[jele];
            covRad[jele] = mACVD.getCovRad()[jele];
        }

        rbondMax = -100000;
        Vector covrads = new Vector();
        for (i = 0; i < NumAt; i++) {
            covalentRadius[i] = 0;
            for (int jele = 0; jele < NumEl; jele++) {
                if (elements[i] == null || elements[i].trim().length() == 0)
                    continue;
                if (elementString[jele].equals(elements[i])) {
                    covalentRadius[i] = bondFactor * covRad[jele] / lenMax;
                }
            }
        }
        for (i = 0; i < covalentRadius.length; i++) {
            if (covalentRadius[i] > rbondMax) {
                rbondMax = covalentRadius[i];
            }
        }

        if (rbondMax <= 0) { // 滅多にないはずだけど.
            return;
        }

        rbondMax *= 2;

        // 小さなセルの大きさ, division数を決めておく.
        ndiv = new int[3];
        smallCell = new double[3];
        logger.debug("rbondMax: " + rbondMax);
        for (i = 0; i < 3; i++) {
            ndiv[i] = (int) (bigCell[i] / rbondMax) + 1;
            smallCell[i] = bigCell[i] / ndiv[i];
        }

        logger.debug("ndiv: " + ndiv[0] + ", " + ndiv[1] + ", " + ndiv[2]);
        logger.debug("smallCell: " + smallCell[0] + ", " + smallCell[1] + ", "
                + smallCell[2]);

        // セルと原子の対応を決めておく
        cell_atom_map = new Vector[ndiv[0]][ndiv[1]][ndiv[2]];
        for (int i1 = 0; i1 < ndiv[0]; i1++) {
            for (int i2 = 0; i2 < ndiv[1]; i2++) {
                for (int i3 = 0; i3 < ndiv[2]; i3++) {
                    cell_atom_map[i1][i2][i3] = null;
                }
            }
        }

        for (int ii = 0; ii < NumAt; ii++) {
            int i1 = (int) (xyz_buff[ii][0] / smallCell[0]);
            int i2 = (int) (xyz_buff[ii][1] / smallCell[1]);
            int i3 = (int) (xyz_buff[ii][2] / smallCell[2]);
            if (cell_atom_map[i1][i2][i3] == null) {
                cell_atom_map[i1][i2][i3] = new Vector();
            }
            logger.debug("ii, i1, i2, i3: " + i1 + ", " + i2 + ", " + i3 + ", "
                    + ii);
            cell_atom_map[i1][i2][i3].addElement(new Integer(ii));
        }

    }

    private boolean[][] ignore;

    public void setIgnore(boolean[][] ignore) {
        this.ignore = ignore;
    }

    public Vector getBonds() {
        return bonds;
    }

    public void doHbond(double[][] xyz) {
        this.xyz = xyz;
        if (!prepareForHbond()) {
            logger.debug("no Hbonds");
            return;
        }
        preProcessForHbond();
        bonds = new Vector();
        if (rbondMax <= 0) {
            return;
        }

        for (int i = 0; i < Hatoms.length; i++) {
            Vector bondInfos = new Vector();
            int i1 = (int) (xyz_buff[i][0] / smallCell[0]);
            int i2 = (int) (xyz_buff[i][1] / smallCell[1]);
            int i3 = (int) (xyz_buff[i][2] / smallCell[2]);
            Vector centerCell = cell_atom_map[i1][i2][i3];
            if (centerCell == null) {
                continue; // このセルには原子はいない
            }
            logger.debug("num. atoms in cell: " + centerCell.size());
            int ind1 = i;

            for (int ii1 = i1 - 1; ii1 <= i1 + 1; ii1++) {
                for (int ii2 = i2 - 1; ii2 <= i2 + 1; ii2++) {
                    for (int ii3 = i3 - 1; ii3 <= i3 + 1; ii3++) {
                        if (ii1 < 0 || ii2 < 0 || ii3 < 0 || ii1 >= ndiv[0]
                                || ii2 >= ndiv[1] || ii3 >= ndiv[2])
                            continue; // 周期境界条件は考慮しない; 差し当たり...
                        Vector vec = cell_atom_map[ii1][ii2][ii3];
                        if (vec == null) { // このセルには原子はいない
                            continue;
                        }
                        for (int ic2 = 0; ic2 < vec.size(); ic2++) {
                            int ind2 = ((Integer) vec.get(ic2)).intValue();
                            if (ind2 < Hatoms.length)
                                continue; // 水素原子

                            if (ignore != null && ignore[map[ind1]][map[ind2]])
                                continue;
                            double rbond2 = Math.pow(covalentRadius[ind2], 2);
                            double r2 = Math.pow(xyz_buff[ind1][0]
                                    - xyz_buff[ind2][0], 2.d)
                                    + Math.pow(xyz_buff[ind1][1]
                                            - xyz_buff[ind2][1], 2.d)
                                    + Math.pow(xyz_buff[ind1][2]
                                            - xyz_buff[ind2][2], 2.d);
                            if (r2 < rbond2 && r2 != 0) {
                                logger.debug("neighboring cell: rbond2, r2 "
                                        + rbond2 + ", " + r2);
                                logger.debug("and their indeces: " + ind1
                                        + ", " + ind2);
                                double[] NearestNeighbor_ = new double[6];
                                NearestNeighbor_[0] = hxyz[ind1][0];
                                NearestNeighbor_[1] = hxyz[ind1][1];
                                NearestNeighbor_[2] = hxyz[ind1][2];
                                NearestNeighbor_[3] = hxyz[ind2][0];
                                NearestNeighbor_[4] = hxyz[ind2][1];
                                NearestNeighbor_[5] = hxyz[ind2][2];
                                BondInfo binf = new BondInfo(ind1, ind2,
                                        NearestNeighbor_);
                                binf.setBondLengthSquared(r2);
                                binf.setCriticalDistanceSquared(rbond2);
                                bondInfos.add(binf);
                            }
                        }
                    }
                }
            }
            if (bondInfos.size() < 2) {
                bondInfos.clear();
                continue;
            }

            BondInfo[] bi = new BondInfo[bondInfos.size()];
            bondInfos.copyInto(bi);
            Arrays.sort(bi);
            double lenmax = mCD.getLenMax();
            logger.debug("bi 0 and 1 " + bi[0].getBondLength() * lenmax + " "
                    + bi[1].getBondLength() * lenmax);
            bonds.addElement(bi[1]);
        }
    }

    private Vector bonds;

    public void doIt(double[][] xyz) {
        this.xyz = xyz;
        preProcess();
        bonds = new Vector();
        if (rbondMax <= 0) {
            return;
        }
        boolean[][][] passed = new boolean[ndiv[0]][ndiv[1]][ndiv[2]];
        for (int i1 = 0; i1 < ndiv[0]; i1++) {
            for (int i2 = 0; i2 < ndiv[1]; i2++) {
                for (int i3 = 0; i3 < ndiv[2]; i3++) {
                    passed[i1][i2][i3] = false;
                }
            }
        }
        for (int i1 = 0; i1 < ndiv[0]; i1++) {
            for (int i2 = 0; i2 < ndiv[1]; i2++) {
                for (int i3 = 0; i3 < ndiv[2]; i3++) {
                    passed[i1][i2][i3] = true;
                    Vector centerCell = cell_atom_map[i1][i2][i3];
                    if (centerCell == null) {
                        continue; // このセルには原子はいない
                    }
                    logger.debug("num. atoms in cell: " + centerCell.size());
                    // 同じセル内は通常のアルゴリズムで.
                    for (int ic1 = 0; ic1 < centerCell.size() - 1; ic1++) {
                        for (int ic2 = ic1 + 1; ic2 < centerCell.size(); ic2++) {
                            int ind1 = ((Integer) centerCell.get(ic1))
                                    .intValue();
                            int ind2 = ((Integer) centerCell.get(ic2))
                                    .intValue();
                            if (ignore != null && ignore[ind1][ind2])
                                continue;
                            double rbond2 = Math.pow(covalentRadius[ind1]
                                    + covalentRadius[ind2], 2);
                            double r2 = Math.pow(xyz_buff[ind1][0]
                                    - xyz_buff[ind2][0], 2.d)
                                    + Math.pow(xyz_buff[ind1][1]
                                            - xyz_buff[ind2][1], 2.d)
                                    + Math.pow(xyz_buff[ind1][2]
                                            - xyz_buff[ind2][2], 2.d);
                            logger.debug("foobar: rbond2, r2 " + rbond2 + ", "
                                    + r2);
                            if (r2 < rbond2) {
                                double[] NearestNeighbor_ = new double[6];
                                NearestNeighbor_[0] = xyz[ind1][0];
                                NearestNeighbor_[1] = xyz[ind1][1];
                                NearestNeighbor_[2] = xyz[ind1][2];
                                NearestNeighbor_[3] = xyz[ind2][0];
                                NearestNeighbor_[4] = xyz[ind2][1];
                                NearestNeighbor_[5] = xyz[ind2][2];
                                BondInfo binf = new BondInfo(ind1, ind2,
                                        NearestNeighbor_);
                                binf.setCriticalDistanceSquared(rbond2);
                                bonds.add(binf);
                            }
                        }
                    }

                    for (int ic1 = 0; ic1 < centerCell.size(); ic1++) {
                        int ind1 = ((Integer) centerCell.get(ic1)).intValue();
                        for (int ii1 = i1 - 1; ii1 <= i1 + 1; ii1++) {
                            for (int ii2 = i2 - 1; ii2 <= i2 + 1; ii2++) {
                                for (int ii3 = i3 - 1; ii3 <= i3 + 1; ii3++) {
                                    if (ii1 < 0 || ii2 < 0 || ii3 < 0
                                            || ii1 >= ndiv[0] || ii2 >= ndiv[1]
                                            || ii3 >= ndiv[2])
                                        continue; // 周期境界条件は考慮しない

                                    if (passed[ii1][ii2][ii3])
                                        continue; // すでに考慮ずみ

                                    Vector vec = cell_atom_map[ii1][ii2][ii3];
                                    if (vec == null) { // このセルには原子はいない
                                        continue;
                                    }
                                    if (ii1 == 0 && ii2 == 0 && ii3 == 0) { // 同じセル内;
                                                                            // 処理済
                                        continue;
                                    }
                                    for (int ic2 = 0; ic2 < vec.size(); ic2++) {
                                        int ind2 = ((Integer) vec.get(ic2))
                                                .intValue();
                                        if (ignore != null
                                                && ignore[ind1][ind2])
                                            continue;
                                        double rbond2 = Math.pow(
                                                covalentRadius[ind1]
                                                        + covalentRadius[ind2],
                                                2);
                                        double r2 = Math.pow(xyz_buff[ind1][0]
                                                - xyz_buff[ind2][0], 2.d)
                                                + Math.pow(xyz_buff[ind1][1]
                                                        - xyz_buff[ind2][1],
                                                        2.d)
                                                + Math.pow(xyz_buff[ind1][2]
                                                        - xyz_buff[ind2][2],
                                                        2.d);
                                        if (r2 < rbond2 && r2 != 0) {
                                            logger.debug("ind1 & 2: " + ind1
                                                    + ", " + ind2);
                                            logger.debug("i1, i2, and i3: "
                                                    + i1 + ", " + i2 + ", "
                                                    + i3);
                                            logger.debug("ii1, ii2, and ii3: "
                                                    + ii1 + ", " + ii2 + ", "
                                                    + ii3);
                                            double[] NearestNeighbor_ = new double[6];
                                            NearestNeighbor_[0] = xyz[ind1][0];
                                            NearestNeighbor_[1] = xyz[ind1][1];
                                            NearestNeighbor_[2] = xyz[ind1][2];
                                            NearestNeighbor_[3] = xyz[ind2][0];
                                            NearestNeighbor_[4] = xyz[ind2][1];
                                            NearestNeighbor_[5] = xyz[ind2][2];
                                            BondInfo binf = new BondInfo(ind1,
                                                    ind2, NearestNeighbor_);
                                            binf.setCriticalDistanceSquared(rbond2);
                                            bonds.addElement(binf);
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        logger.debug("num. bonds: " + bonds.size());
    }

}
