/*
!=======================================================================
!
!  PROGRAM  PHASE-Viewer  (PHASE-Viewer 2014.01 ver.3.3.0)
!
!  Created on 2006/12/12 15:40:29
!  AUTHOR(S): KOGA, Junichiro
!  File : BravaisLattice.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.atomcoord.symmetry;

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.util.Arrays;
import java.util.HashMap;
import java.util.Vector;

import Jama.Matrix;
import ciss.phase_viewer.atomcoord.Atom;
import ciss.phase_viewer.atomcoord.AtomCoords;
import ciss.phase_viewer.atomcoord.AtomList;
import ciss.phase_viewer.atomcoord.VolumetricData;
import ciss.phase_viewer.common.ConstParameters;
import ciss.phase_viewer.inputinterface.InputInterface;
import ciss.phase_viewer.inputinterface.InputInterfaceEntryChangeEvent;
import ciss.phase_viewer.inputinterface.InputInterfaceEntryChangeListener;
import ciss.phase_viewer.inputinterface.InputInterfacePrimitiveEntry;
import ciss.phase_viewer.main.PluginLoader;

/**
 * Bravais格子を表すクラスのベースクラス.
 * 
 * @author KOGA, Junichiro
 * 
 */
public abstract class BravaisLattice implements
        InputInterfaceEntryChangeListener {
    /*
     * (non-Javadoc)
     * 
     * @see ciss.phase_viewer.inputinterface.InputInterfaceEntryChangeListener#
     * inputInterfaceEntryChanged
     * (ciss.phase_viewer.inputinterface.InputInterfaceEntryChangeEvent)
     */
    public void inputInterfaceEntryChanged(InputInterfaceEntryChangeEvent e) {
    }

    /*
     * (non-Javadoc)
     * 
     * @see ciss.phase_viewer.inputinterface.InputInterfaceEntryChangeListener#
     * inputInterfaceInitialized()
     */
    public void inputInterfaceInitialized() {
    }

    /*
     * (non-Javadoc)
     * 
     * @see ciss.phase_viewer.inputinterface.InputInterfaceEntryChangeListener#
     * inputInterfaceInitializing()
     */
    public void inputInterfaceInitializing() {
    }

    private static org.apache.log4j.Logger logger = org.apache.log4j.Logger
            .getLogger(BravaisLattice.class.getName());

    public static final String PRIMITIVE = "primitive";

    public static final String RHOMBOHEDRAL = "rhombohedral";

    public static final String HEXAGONAL = "hexagonal";

    public static final String BODY_CENTERED = "bodycentered";

    public static final String BASE_CENTERED = "basecentered";

    public static final String FACE_CENTERED = "facecentered";

    protected String latticeSystem = PRIMITIVE;

    protected String[] latticeSystems = { PRIMITIVE, RHOMBOHEDRAL, HEXAGONAL,
            BODY_CENTERED, BASE_CENTERED, FACE_CENTERED };

    protected double a;

    protected double b;

    protected double c;

    protected double alpha;

    protected double beta;

    protected double gamma;

    protected HashMap operatorMap;

    protected BravaisLattice() {
        operatorMap = new HashMap();
        operatorMap.put(FACE_CENTERED, new double[][] { { 0.5, 0.5, 0 },
                { 0.5, 0, 0.5 }, { 0, 0.5, 0.5 } });
        operatorMap.put(BODY_CENTERED, new double[][] { { 0.5, 0.5, 0.5 } });
        operatorMap.put(BASE_CENTERED, new double[][] { { 0.5, 0.5, 0 } });
        operatorMap.put(RHOMBOHEDRAL, new double[][] {
                { 2.d / 3.d, 1.d / 3.d, 1.d / 3.d },
                { 1.d / 3.d, 2.d / 3.d, 2.d / 3.d } });
    }

    /**
     * lattice systemに応じた演算子の配列を返す.
     * 
     * @return 上記の通り.
     */
    public double[][] getOperators() {
        Object obj = operatorMap.get(latticeSystem);
        if (obj == null)
            return null;

        return (double[][]) obj;
    }

    /**
     * ブラベー格子の名前を取得.
     * 
     * @return ブラベー格子の名前
     */
    public abstract String getName();

    /**
     * 結晶系の名前を取得.
     * 
     * @return 結晶系の名前
     */
    public abstract String getCrystalSystem();

    public String getLatticeSystem() {
        return latticeSystem;
    }

    public void setLatticeSystem(String latticeSystem) {
        this.latticeSystem = latticeSystem;
    }

    /**
     * 基本格子ベクトルを計算, 取得.
     * 
     * @return
     */
    public abstract double[][] getPrimitiveLatticeVector();

    /**
     * 対応する結晶系の可能なlattice systemを返す.
     */
    public abstract String[] getLatticeSystemCandidate();

    /**
     * 基本格子ベクトルからブラベー格子のインスタンスを作成.
     * 
     * @param avec
     *            a軸のベクトル
     * @param bvec
     *            b軸のベクトル
     * @param cvec
     *            c軸のベクトル
     * @return 結果作成されたブラベー格子; うまく行かなかった場合null.
     */
    public static BravaisLattice getBravaisLattice(double[] avec,
            double[] bvec, double[] cvec) {
        BufferedReader breader = null;
        InputStream is = null;
        logger.debug("creating bravais lattice");
        double[] cellparams = BravaisLattice
                .getCellParameters(avec, bvec, cvec);
        logger.debug("a, b, c, alpha, beta, gamma: " + cellparams[0] + ", "
                + cellparams[1] + ", " + cellparams[2] + ", " + cellparams[3]
                + ", " + cellparams[4] + ", " + cellparams[5]);
        try {
            is = BravaisLattice.class
                    .getResourceAsStream("/ciss/phase_viewer/atomcoord/symmetry/bravaislattice.properties");
            breader = new BufferedReader(new InputStreamReader(is));
            String line = "";
            while ((line = breader.readLine()) != null) {
                line = line.trim();
                if (line.startsWith("#")) {
                    continue;
                }
                String[] foo = line.split("=");
                if (foo == null || foo.length < 2) {
                    continue;
                }

                String path = foo[1];
                BravaisLattice lat = (BravaisLattice) PluginLoader.instantiate(
                        path, null);
                lat.a = cellparams[0];
                lat.b = cellparams[1];
                lat.c = cellparams[2];
                lat.alpha = cellparams[3];
                lat.beta = cellparams[4];
                lat.gamma = cellparams[5];
                BravaisLattice ret = lat.tryToCreateBravaisLattice(avec, bvec,
                        cvec);
                if (ret != null)
                    return ret;
            }
        } catch (IOException ioe) {
            ioe.printStackTrace();
        } finally {
            try {
                breader.close();
                is.close();
            } catch (IOException e) {
                e.printStackTrace();
            }
        }

        logger.warn("no matching crystal!!!");
        BravaisLattice lat = new Generic();
        lat.a = cellparams[0];
        lat.b = cellparams[1];
        lat.c = cellparams[2];
        lat.alpha = cellparams[3];
        lat.beta = cellparams[4];
        lat.gamma = cellparams[5];
        return lat;
    }

    /**
     * 引数の条件からブラベー格子のインスタンスを作成.
     * 
     * @param a
     * @param b
     * @param c
     * @param alpha
     * @param beta
     * @param gamma
     * @param lattice_system
     *            nullの場合PRIMITIVE.
     * @return
     */
    public static BravaisLattice getBravaisLattice(double a, double b,
            double c, double alpha, double beta, double gamma,
            String latticeSystem) {
        BufferedReader breader = null;
        InputStream is = null;
        logger.debug("creating bravais lattice");
        try {
            is = BravaisLattice.class
                    .getResourceAsStream("/ciss/phase_viewer/atomcoord/symmetry/bravaislattice.properties");
            breader = new BufferedReader(new InputStreamReader(is));
            String line = "";
            while ((line = breader.readLine()) != null) {
                line = line.trim();
                if (line.startsWith("#")) {
                    continue;
                }
                String[] foo = line.split("=");
                if (foo == null || foo.length < 2) {
                    continue;
                }
                String na = foo[0];
                String path = foo[1];
                BravaisLattice lat = (BravaisLattice) PluginLoader.instantiate(
                        path, null);
                if (lat.isThisForMe(a, b, c, alpha, beta, gamma)) {
                    lat.a = a;
                    lat.b = b;
                    lat.c = c;
                    lat.alpha = alpha;
                    lat.beta = beta;
                    lat.gamma = gamma;
                    lat.latticeSystem = latticeSystem;
                    logger.debug("created bravais lattice : " + lat);
                    return lat;
                }
            }
        } catch (IOException ioe) {
            ioe.printStackTrace();
        } finally {
            try {
                breader.close();
                is.close();
            } catch (IOException e) {
                e.printStackTrace();
            }
        }

        logger.warn("no matching crystal!!!");
        BravaisLattice lat = new Primitive();
        lat.a = a;
        lat.b = b;
        lat.c = c;
        lat.alpha = alpha;
        lat.beta = beta;
        lat.gamma = gamma;
        lat.isBravais(false);
        return lat;
    }

    protected abstract boolean isThisForMe(double a, double b, double c,
            double alpha, double beta, double gamma);

    /**
     * 基本格子ベクトルからブラベー格子作成を試みる. うまくいかない場合null.
     * 
     * @param avec
     * @param bvec
     * @param cvec
     * @return
     */
    protected abstract BravaisLattice tryToCreateBravaisLattice(double[] avec,
            double[] bvec, double[] cvec);

    private double epsilon = 1.0e-3;

    protected boolean epsilonEquals(double a, double b) {
        if (Math.abs(a - b) < epsilon)
            return true;
        return false;
    }

    public String toString() {
        return getName();
    }

    /**
     * 指定のInputInterfaceから, ある座標データがPrimitiveなのかBravaisなのかを判定するメソッド.
     * 
     * @param inputInterface
     *            このInputInterfaceから判定する.
     * @return Bravaisの場合真.
     */
    public static boolean isBravais(InputInterface inputInterface) {
        InputInterfacePrimitiveEntry uni = inputInterface
                .getInputInterfacePrimitiveEntry("structure.unit_cell_type");
        if (uni.getValue().equalsIgnoreCase("primitive")) {
            return false;
        }

        return true;
    }

    /**
     * 指定のInputInterfaceから, 格子の型を判定して結果を返す.
     * structure.symmetry.crystal_structureとstructure
     * .symmetry.tspace.lattice_systemの二種類があるけど, tspaceの方を優先する? べきかどうかは分からない.
     * 
     * @param inputInterface
     *            このInputInterfaceオブジェクトから判定する.
     * @return 格子の型; このクラスのstatic 文字列のいずれかが返される.
     */
    public static String getLatticeSystemFrom(InputInterface inputInterface) {
        InputInterfacePrimitiveEntry crystal_symmetry = inputInterface
                .getInputInterfacePrimitiveEntry("structure.symmetry.crystal_structure");
        if (crystal_symmetry.getValue().equalsIgnoreCase("fcc")
                || crystal_symmetry.getValue().equalsIgnoreCase("diamond"))
            return FACE_CENTERED;
        if (crystal_symmetry.getValue().equalsIgnoreCase("bcc"))
            return BODY_CENTERED;

        InputInterfacePrimitiveEntry lattice_system = inputInterface
                .getInputInterfacePrimitiveEntry("structure.symmetry.tspace.lattice_system");
        String val = lattice_system.getValue();
        if (val.equalsIgnoreCase("rhombohedral")
                || val.equalsIgnoreCase("trigonal")
                || val.equalsIgnoreCase("r") || val.equalsIgnoreCase("t")
                || val.equalsIgnoreCase("-1"))
            return RHOMBOHEDRAL;
        if (val.equalsIgnoreCase("hexagonal") || val.equalsIgnoreCase("h")
                || val.equalsIgnoreCase("1"))
            return HEXAGONAL;
        if (val.equalsIgnoreCase("primitive") || val.equalsIgnoreCase("simple")
                || val.equalsIgnoreCase("p") || val.equalsIgnoreCase("s"))
            return PRIMITIVE;
        if (val.equalsIgnoreCase("facecentered") || val.equalsIgnoreCase("f")
                || val.equalsIgnoreCase("2"))
            return FACE_CENTERED;
        if (val.equalsIgnoreCase("bodycentered") || val.equalsIgnoreCase("b")
                || val.equalsIgnoreCase("3"))
            return BODY_CENTERED;
        if (val.equalsIgnoreCase("bottomcentered")
                || val.equalsIgnoreCase("basecentered")
                || val.equalsIgnoreCase("onefacecentered")
                || val.equalsIgnoreCase("bot") || val.equalsIgnoreCase("ba")
                || val.equalsIgnoreCase("o") || val.equalsIgnoreCase("4"))
            return BASE_CENTERED;
        return PRIMITIVE;
    }

    /**
     * 基本格子ベクトルから, a,b,c,alpha,beta,gammaをもとめて返す.
     * 
     * @param avec
     * @param bvec
     * @param cvec
     * @return
     */
    public static double[] getCellParameters(double[] avec, double[] bvec,
            double[] cvec) {
        double[] ret = new double[6];

        double ar = 0.;
        double br = 0.;
        double cr = 0.;
        for (int i = 0; i < 3; i++) {
            ar += Math.pow(avec[i], 2);
            br += Math.pow(bvec[i], 2);
            cr += Math.pow(cvec[i], 2);
        }
        ar = Math.sqrt(ar);
        br = Math.sqrt(br);
        cr = Math.sqrt(cr);

        double cda = 0.;
        double cdb = 0.;
        double adb = 0.;

        for (int i = 0; i < 3; i++) {
            cda += cvec[i] * avec[i];
            cdb += cvec[i] * bvec[i];
            adb += avec[i] * bvec[i];
        }

        double cosbeta = 0.;
        double cosalph = 0.;
        double cosgamm = 0.;

        double beta = 0.;
        double alph = 0.;
        double gamm = 0.;

        try {
            cosbeta = cda / (cr * ar);
            cosalph = cdb / (cr * br);
            cosgamm = adb / (ar * br);

            beta = Math.acos(cosbeta) * ConstParameters.RAD_TO_ANG;
            alph = Math.acos(cosalph) * ConstParameters.RAD_TO_ANG;
            gamm = Math.acos(cosgamm) * ConstParameters.RAD_TO_ANG;
        } catch (Exception e) {
            return null;
        }

        ret[0] = ar;
        ret[1] = br;
        ret[2] = cr;
        ret[3] = alph;
        ret[4] = beta;
        ret[5] = gamm;
        return ret;
    }

    /**
     * ブラベー格子の格子定数から, セルベクトルを取得する.
     * 
     * @return
     */
    public double[][] getGenericPrimitiveVector() {
        double[] cellLength = new double[] { a, b, c };
        double alpha = this.alpha * ConstParameters.ANG_TO_RAD;
        double beta = this.beta * ConstParameters.ANG_TO_RAD;
        double gamma = this.gamma * ConstParameters.ANG_TO_RAD;
        double temp = 0.;
        double[][] cellVec = new double[3][3];

        logger.debug("cellLength: " + cellLength[0] + " " + cellLength[1] + " "
                + cellLength[2]);
        for (int i = 0; i < 3; ++i)
            for (int j = 0; j < 3; ++j)
                cellVec[i][j] = 0.;

        cellVec[0][0] = 1.;
        cellVec[1][0] = Math.cos(gamma);
        cellVec[2][0] = Math.cos(beta);
        cellVec[1][1] = Math.sin(gamma);
        temp = (Math.cos(alpha) - Math.cos(beta) * Math.cos(gamma))
                / Math.sin(gamma);
        cellVec[2][1] = temp;
        cellVec[2][2] = Math
                .sqrt(Math.sin(beta) * Math.sin(beta) - temp * temp);
        for (int i = 0; i < 3; ++i)
            for (int j = 0; j < 3; ++j)
                cellVec[j][i] *= cellLength[j];

        return cellVec;
    }

    /**
     * 副格子点を追加する. 追加された原子は, 保存の対象になることはない.
     * 
     * @param coords
     *            副格子点を追加したい原子配置
     */
    public void addSubLatticeAtoms(AtomCoords coords) {
        if (latticeSystem.equals(PRIMITIVE) || latticeSystem.equals(HEXAGONAL))
            return;

        AtomList atomList = coords.getAtomList();
        double[][] cellvec = getGenericPrimitiveVector();
        double[][] operations = getOperators();

        if (operations == null)
            return;

        int undoCount = 0;
        int numat = atomList.size();
        for (int i = 0; i < numat; i++) {
            Atom atom = atomList.getAtomAt(i);
            if (atom.isBoundaryAtom() || atom.isSubLatticeAtom())
                continue;
            double[] pos = atom.getDouble();
            for (int j = 0; j < operations.length; j++) {
                double[] newpos = new double[3];
                for (int l = 0; l < 3; l++) {
                    newpos[l] = pos[l];
                    for (int k = 0; k < 3; k++) {
                        newpos[l] += operations[j][k] * cellvec[k][l];
                    }
                }
                Atom newAtom = new Atom(atom.getElementName(), newpos,
                        atom.getAuxil());
                newAtom.isSubLatticeAtom(true);
                atomList.addAtom(newAtom);
                undoCount++;
            }
        }
        atomList.setUndoCount(undoCount);
    }

    /**
     * lattice systemを推測するメソッド. ただし, 現在PRIMITVEかHEXAGONALの場合以外で, かつ結晶系にあるlattice
     * systemの場合はそちらを優先する.
     * 
     * @param coords
     */
    public void estimateLatticeSystem(AtomCoords coords) {
        String currlat = new String(latticeSystem);
        String[] latsys = getLatticeSystemCandidate();
        if (!(currlat.equalsIgnoreCase(PRIMITIVE) || currlat.equals(HEXAGONAL))) {
            for (int i = 0; i < latsys.length; i++)
                if (latsys[i].equals(currlat))
                    return;
        }
        logger.debug("estimating lattice system");

        for (int i = 0; i < latsys.length; i++) {
            Vector tmpvec = new Vector();
            if (doLatticeSystem(coords, latsys[i], tmpvec)) {
                this.latticeSystem = latsys[i];
                logger.debug("new lattice system: " + this.latticeSystem);
                return;
            }
        }
    }

    /**
     * 副格子点に相当すると思われる原子を削除する(PHASEの入力を作成する際などに利用する) 同時に, lattice systemを推定する.
     */
    public void stripSubLatticeAtoms(AtomCoords coords) {
        String currlat = new String(latticeSystem);
        if (stripSubLatticeAtoms(coords, currlat))
            return;

        String[] latsys = getLatticeSystemCandidate();
        for (int sys = 0; sys < latsys.length; sys++) {
            if (latsys[sys].equals(currlat))
                continue;
            if (stripSubLatticeAtoms(coords, latsys[sys])) {
                latticeSystem = latsys[sys];
                logger.debug("new lattice system: " + latticeSystem);
                return;
            }
        }

    }

    /**
     * InputInterfaceの設定を行う.
     * 
     * @param inputInterface
     *            対象のInputInterfaceオブジェクト
     */
    public void saveToInputInterface(InputInterface inputInterface) {
        logger.debug("isBravais(): " + isBravais);
        InputInterfacePrimitiveEntry brava = inputInterface
                .getInputInterfacePrimitiveEntry("structure.unit_cell_type");
        if (isBravais())
            brava.setValue("bravais");
        else
            brava.setValue("primitive");
        inputInterface.replaceEntry(brava, this, true);

        InputInterfacePrimitiveEntry lats = inputInterface
                .getInputInterfacePrimitiveEntry("structure.symmetry.tspace.lattice_system");
        lats.setValue(latticeSystem);
        inputInterface.replaceEntry(lats, this, true);
    }

    private double[][] pm = { { 1, 1, 1 }, { -1, 1, 1 }, { 1, -1, 1 },
            { 1, 1, -1 }, { -1, -1, 1 }, { -1, 1, -1 }, { 1, -1, -1 },
            { -1, -1, -1 } };

    /**
     * 副格子点に相当すると思われる原子を削除する(PHASEの入力を作成する際などに利用する)
     * 
     * @param coords
     *            対象となる座標データ
     * @param latsys
     *            latticesystemの値.
     * @return stripされた場合true.
     */
    public boolean stripSubLatticeAtoms(AtomCoords coords, String latsys) {
        Vector remvec = new Vector();
        boolean strip = doLatticeSystem(coords, latsys, remvec);

        int remcount = 0;
        AtomList li = coords.getAtomList();
        if (strip) {
            Vector indeces = new Vector();
            for (int i = 0; i < remvec.size(); i++) {
                int index = li.indexOf(remvec.get(i));
                indeces.addElement(new Integer(index));
            }
            Integer[] ints = new Integer[indeces.size()];
            indeces.copyInto(ints);
            Arrays.sort(ints);
            for (int i = ints.length - 1; i >= 0; i--) {
                int index = ints[i].intValue();
                logger.debug("index: " + index);
                if (index < li.getNumAt())
                    li.removeAtomAt(index);
                remcount++;
            }
            logger.debug("num. at: " + li.getNumAt());
            li.setUndoCount(remcount);
        }

        return strip;
    }

    private boolean doLatticeSystem(AtomCoords coords, String latsys,
            Vector remvec) {
        AtomList li = coords.getAtomList();
        AtomList copy = li.getCopy();

        double[][] cellvec = getGenericPrimitiveVector();
        Object obj = operatorMap.get(latsys);
        if (obj == null)
            return false;

        double[][] operators = (double[][]) obj;
        iloop: for (int i = 0; i < copy.getNumAt() - 1; i++) {
            Atom at = copy.getAtomAt(i);
            if (at.isVirtualAtom())
                continue iloop;
            jloop: for (int j = copy.getNumAt() - 1; j >= i + 1; j--) {
                Atom compat = copy.getAtomAt(j);
                if (compat.isVirtualAtom())
                    continue jloop;
                for (int iope = 0; iope < operators.length; iope++) {
                    for (int ipm = 0; ipm < pm.length; ipm++) {
                        double[] tmpop = new double[3];
                        for (int l = 0; l < 3; l++)
                            tmpop[l] = pm[ipm][l] * operators[iope][l];

                        double[] newpos = new double[3];
                        double[] pos = compat.getDouble();
                        for (int l = 0; l < 3; l++) {
                            newpos[l] = pos[l];
                            for (int k = 0; k < 3; k++) {
                                newpos[l] += tmpop[k] * cellvec[k][l];
                            }
                        }

                        Atom at2 = new Atom(compat.getElementName(), newpos);
                        if (at2.epsilonEquals(at)) {
                            if (copy.removeAtom(compat, false))
                                remvec.add(compat);
                        }
                    }
                }
            }
        }

        int numat = copy.getRealNumAt();
        int numorg = numat * (operators.length + 1);
        logger.debug("numat " + numat + ", numorg " + numorg + ", "
                + li.getRealNumAt());
        boolean strip = numorg == li.getRealNumAt();
        return strip;
    }

    private boolean isBravais = true;

    /**
     * 本当にブラベー格子か否かを設定.
     * 
     * @param isBravais
     *            ブラベーならtrue.
     */
    public void isBravais(boolean isBravais) {
        this.isBravais = isBravais;
    }

    /**
     * 本当にブラベー格子ならtrue.
     * 
     * @return 上記の通り
     */
    public boolean isBravais() {
        return this.isBravais;
    }

    /**
     * Bravaisセル用のVolumetricDataからPrimitiveセル用のVolumetricDataを作成して返す.
     * Bravaisは現Bravaisセル, Primitiveは対応するPrimitiveセル. ただし, latticeSystem
     * がPRIMITIVEかHEXAGONALの時は結果をそのまま返す.
     * 
     * @param origData
     * @return
     */
    public VolumetricData getPrimitiveVolumetricDataFromBravais(
            VolumetricData origData) {
        if (latticeSystem.equals(PRIMITIVE) || latticeSystem.equals(HEXAGONAL))
            return origData;

        origData.stripPBC();
        double[][] primlat = getPrimitiveLatticeVector();

        int[] origNdiv = origData.getNumDiv();
        float[][] origDelta = origData.getDelta();
        float[] origDensity = origData.getDensity();
        float[] origOrigin = origData.getOrigin();

        int[] newNdiv = new int[3];
        float[][] newDelta = new float[3][3];
        float[] newOrigin = new float[3];

        for (int i = 0; i < 3; i++)
            newNdiv[i] = origNdiv[i];

        for (int i = 0; i < 3; i++)
            newOrigin[i] = origOrigin[i];

        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                newDelta[i][j] = (float) (primlat[i][j] / (double) newNdiv[i]);

        logger.debug("origNdiv: " + origNdiv[0] + ", " + origNdiv[1] + ", "
                + origNdiv[2]);
        logger.debug("newNdiv: " + newNdiv[0] + ", " + newNdiv[1] + ", "
                + newNdiv[2]);

        logger.debug("origDelta1: " + origDelta[0][0] + ", " + origDelta[0][1]
                + ", " + origDelta[0][2]);
        logger.debug("origDelta2: " + origDelta[1][0] + ", " + origDelta[1][1]
                + ", " + origDelta[1][2]);
        logger.debug("origDelta3: " + origDelta[2][0] + ", " + origDelta[2][1]
                + ", " + origDelta[2][2]);

        logger.debug("newDelta1: " + newDelta[0][0] + ", " + newDelta[0][1]
                + ", " + newDelta[0][2]);
        logger.debug("newDelta2: " + newDelta[1][0] + ", " + newDelta[1][1]
                + ", " + newDelta[1][2]);
        logger.debug("newDelta3: " + newDelta[2][0] + ", " + newDelta[2][1]
                + ", " + newDelta[2][2]);

        float[] newDensity = new float[newNdiv[0] * newNdiv[1] * newNdiv[2]];

        double[][] tmpPoint = new double[3][1];

        double[][] tmpprimlat = new double[3][3];
        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                tmpprimlat[j][i] = primlat[i][j];
        Matrix primMatrix = new Matrix(tmpprimlat);
        Matrix pointMatrix;
        int[] tmpind = new int[3];
        for (int i = 0; i < origNdiv[2]; i++) {
            for (int j = 0; j < origNdiv[1]; j++) {
                for (int k = 0; k < origNdiv[0]; k++) {
                    for (int l = 0; l < 3; l++)
                        tmpPoint[l][0] = i * origDelta[2][l] + j
                                * origDelta[1][l] + k * origDelta[0][l];
                    pointMatrix = new Matrix(tmpPoint);
                    Matrix keisu = primMatrix.solve(pointMatrix);
                    for (int l = 0; l < 3; l++) {
                        double foo = keisu.get(l, 0)
                                - Math.floor(keisu.get(l, 0));
                        tmpind[l] = (int) Math.round(foo * newNdiv[l]);
                        if (tmpind[l] >= newNdiv[l])
                            tmpind[l] = 0;
                    }

                    int newtmpind = tmpind[0] * newNdiv[2] * newNdiv[1]
                            + tmpind[1] * newNdiv[2] + tmpind[2];
                    int oldtmpind = k * origNdiv[2] * origNdiv[1] + j
                            * origNdiv[2] + i;

                    if (newtmpind >= newDensity.length) {
                        System.out.println("ind " + tmpind[0] + ", "
                                + tmpind[1] + ", " + tmpind[2] + " 1dind: "
                                + newtmpind);
                        continue;
                    }

                    newDensity[newtmpind] = origDensity[oldtmpind];
                }
            }
        }

        origData.doPBC(origDensity, origNdiv);

        VolumetricData vdata = new VolumetricData(newDensity, newNdiv,
                newDelta, newOrigin);
        return vdata;
    }

    /**
     * Primitiveセル用のVolumetricDataからBravaisセル用のVolumetricDataに変換して返す.
     * Bravaisは現Bravaisセル, Primitiveは対応するPrimitiveセルとする. ただし, latticeSystem
     * がPRIMITIVEかHEXAGONALの時は結果をそのまま返す.
     * 
     * @param origData
     *            Primitiveセル用のVolumetricData.
     * @return Bravaisセル用のVolumetricData(origDataと同じ参照)
     */
    public VolumetricData bravaisVolumetricData2Primitive(
            VolumetricData origData) {
        if (latticeSystem.equals(PRIMITIVE) || latticeSystem.equals(HEXAGONAL))
            return origData;

        origData.stripPBC();

        double[][] primlat = getPrimitiveLatticeVector();
        double[][] bravlat = getGenericPrimitiveVector();

        int[] origNdiv = origData.getNumDiv();
        float[][] origDelta = origData.getDelta();
        float[] origDensity = origData.getDensity();
        float[] origOrigin = origData.getOrigin();

        int[] newNdiv = new int[3];
        float[][] newDelta = new float[3][3];
        float[] newOrigin = new float[3];

        for (int i = 0; i < 3; i++)
            // newNdiv[i] = (int) Math.round((double) origNdiv[i]
            // * VectorOperations.norm(bravlat[i])
            // / VectorOperations.norm(primlat[i]));
            newNdiv[i] = origNdiv[i];

        for (int i = 0; i < 3; i++) {
            newOrigin[i] = origOrigin[i];
        }

        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                newDelta[i][j] = (float) (bravlat[i][j] / (double) (newNdiv[i]));

        logger.debug("origNdiv: " + origNdiv[0] + ", " + origNdiv[1] + ", "
                + origNdiv[2]);
        logger.debug("newNdiv: " + newNdiv[0] + ", " + newNdiv[1] + ", "
                + newNdiv[2]);

        logger.debug("origDelta1: " + origDelta[0][0] + ", " + origDelta[0][1]
                + ", " + origDelta[0][2]);
        logger.debug("origDelta2: " + origDelta[1][0] + ", " + origDelta[1][1]
                + ", " + origDelta[1][2]);
        logger.debug("origDelta3: " + origDelta[2][0] + ", " + origDelta[2][1]
                + ", " + origDelta[2][2]);

        logger.debug("newDelta1: " + newDelta[0][0] + ", " + newDelta[0][1]
                + ", " + newDelta[0][2]);
        logger.debug("newDelta2: " + newDelta[1][0] + ", " + newDelta[1][1]
                + ", " + newDelta[1][2]);
        logger.debug("newDelta3: " + newDelta[2][0] + ", " + newDelta[2][1]
                + ", " + newDelta[2][2]);

        // return origData;

        float[] newDensity = new float[newNdiv[0] * newNdiv[1] * newNdiv[2]];

        logger.debug("origDensitySize: " + origDensity.length);
        logger.debug("newDensitySize: " + newDensity.length);

        double[][] tmpPoint = new double[3][1];

        double[][] tmpprimlat = new double[3][3];
        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                tmpprimlat[j][i] = primlat[i][j];
        Matrix primMatrix = new Matrix(tmpprimlat);
        Matrix pointMatrix;
        int[] tmpind = new int[3];
        for (int i = 0; i < newNdiv[2]; i++) {
            for (int j = 0; j < newNdiv[1]; j++) {
                for (int k = 0; k < newNdiv[0]; k++) {
                    for (int l = 0; l < 3; l++)
                        tmpPoint[l][0] = i * newDelta[2][l] + j
                                * newDelta[1][l] + k * newDelta[0][l];
                    pointMatrix = new Matrix(tmpPoint);
                    Matrix keisu = primMatrix.solve(pointMatrix);
                    for (int l = 0; l < 3; l++) {
                        double foo = keisu.get(l, 0)
                                - Math.floor(keisu.get(l, 0));
                        tmpind[l] = (int) Math.round(foo * origNdiv[l]);
                        if (tmpind[l] >= origNdiv[l])
                            tmpind[l] = 0;
                    }

                    int oldtmpind = tmpind[0] * origNdiv[2] * origNdiv[1]
                            + tmpind[1] * origNdiv[2] + tmpind[2];
                    int newtmpind = k * newNdiv[2] * newNdiv[1] + j
                            * newNdiv[2] + i;

                    if (oldtmpind >= origDensity.length) {
                        System.out.println("ind " + tmpind[0] + ", "
                                + tmpind[1] + ", " + tmpind[2] + " 1dind: "
                                + oldtmpind);
                    }

                    newDensity[newtmpind] = origDensity[oldtmpind];
                }
            }
        }

        // origData.doPBC(origDensity, origNdiv);
        origData.setDelta(newDelta);
        origData.doPBC(newDensity, newNdiv);

        origData.createBuffer(newDensity, newNdiv);
        return origData;
        // VolumetricData vdata = new VolumetricData(newDensity, newNdiv,
        // newDelta, newOrigin);
        // return vdata;
    }

}
