/*
!=======================================================================
!
!  PROGRAM  PHASE-Viewer  (PHASE-Viewer 2014.01 ver.3.3.0)
!
!  Created on 2005/11/10, 15:24
!  AUTHOR(S): KOGA, Junichiro
!  File : VolumetricData.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;

import java.util.Vector;

import org.apache.log4j.Logger;

import Jama.Matrix;
import ciss.phase_viewer.acviewer.colormap.BaseColorMap;
import ciss.phase_viewer.acviewer.scenegraphelements.abinitmp.MarchingCube;
import ciss.phase_viewer.common.BinaryOperator;

/**
 * 「電荷密度情報」を格納するクラス. GaussianCubeファイルの内容を表すことを目的とした クラスだが, 不完全であることに注意.
 * 
 * @author
 */
public class VolumetricData implements BinaryOperator {

    public static int CHARGE_DENSITY = 0;
    /** 電荷密度を表すボリュームデータの場合setModeでこの値を設定する. */
    public static int MOLECULAR_ORBITAL = 1;

    /** 分子軌道を表すボリュームデータの場合setModeでこの値を設定する. */

    /*
     * (non-Javadoc)
     * 
     * @see
     * ciss.phase_viewer.common.BinaryOperator#add(ciss.phase_viewer.common.
     * BinaryOperator)
     */
    public void add(BinaryOperator bo) {
        if (!(bo instanceof VolumetricData))
            return;
        VolumetricData vdata = (VolumetricData) bo;
        if (vdata.density.length != this.density.length)
            return;

        for (int i = 0; i < this.density.length; i++)
            this.density[i] += vdata.density[i];
    }

    /**
     * 実装しない.
     */
    public void divide(BinaryOperator bo) {
    }

    /**
     * 実装しない.
     */
    public void multiply(BinaryOperator bo) {
    }

    /**
     * @see ciss.phase_viewer.common.BinaryOperator#subtract(ciss.phase_viewer.common.BinaryOperator)
     */
    public void subtract(BinaryOperator bo) {
        if (!(bo instanceof VolumetricData))
            return;
        VolumetricData vdata = (VolumetricData) bo;
        if (vdata.density.length != this.density.length)
            return;

        for (int i = 0; i < this.density.length; i++)
            this.density[i] -= vdata.density[i];
    }

    private int mode = CHARGE_DENSITY;

    /**
     * ボリュームデータの"mode"を設定;
     * 
     * @param mode
     *            設定したいmode, CHARGE_DENSITYかMOLECULAR_ORBITAL
     */
    public void setMode(int mode) {
        this.mode = mode;
        this.interpolationScheme = VolumetricData.LINEAR;
    }

    /**
     * このボリュームデータの"mode"を取得.
     * 
     * @return このボリュームデータのmode; CHARGE_DENSITY かMOLECULAR_ORBITAL.
     */
    public int getMode() {
        return this.mode;
    }

    private Vector molecularOrbitalIndeces = null;

    /**
     * 分子軌道のインデックスを足す. 単なる整数だけど...
     * 
     * @param obj
     *            分子軌道のインデックス
     */
    public void addMolecularOrbitalIndex(Object obj) {
        if (molecularOrbitalIndeces == null) {
            molecularOrbitalIndeces = new Vector();
            molecularOrbitalIndeces.add(obj);
        }
    }

    /**
     * 分子軌道のインデックスを削除.
     * 
     * @param obj
     *            削除したいインデックス
     */
    public void removeMolecularOrbitalIndex(Object obj) {
        if (molecularOrbitalIndeces != null)
            molecularOrbitalIndeces.remove(obj);
    }

    public Object[] getMolecularOrbitalIndeces() {
        if (molecularOrbitalIndeces == null)
            return null;

        Object[] ob = new Object[molecularOrbitalIndeces.size()];
        molecularOrbitalIndeces.copyInto(ob);
        return ob;
    }

    private Logger logger = Logger.getLogger(VolumetricData.class.getName());

    private float[] density;

    private float[] density_buff;

    private float[] sortedDensity;

    private int[] ndiv;

    private int[] ndiv_buff;

    private float[][] delta;

    private float[][] delta_buff;

    private float[] origin = new float[] { 0, 0, 0 };

    private int Ntot;

    private float minVal = 10000.f;

    private float maxVal = -10000.f;

    private float midVal = 0.0f;

    public static int LOG = MarchingCube.LOG;

    public static int LINEAR = MarchingCube.LINEAR;

    private int interpolationScheme = LOG;

    /** この値と等しい値を見つけたら無視 */
    public static float meaninglessValue = 1000f;

    private boolean pbc = true;

    /**
     * 電荷密度を表すに必要な情報を渡す.
     * 
     * @param density
     *            空間上の各点での電荷密度の値そのもの.
     * @param ndiv
     *            a軸, b軸, c軸の分割数
     * @param delta
     *            a軸, b軸, c軸の「変化量」
     * @param origin
     *            原点を指定する配列.
     */
    public VolumetricData(float[] density, int[] ndiv, float[][] delta,
            float[] origin) {
        if (!pbc) {
            this.density = density;
            this.ndiv = ndiv;
            this.delta = delta;
            this.origin = origin;
            Ntot = ndiv[0] * ndiv[1] * ndiv[2];
            this.init(density, ndiv);
        } else {
            this.delta = delta;
            this.origin = origin;
            this.density = density;
            this.ndiv = ndiv;
            this.Ntot = ndiv[0] * ndiv[1] * ndiv[2];
            doPBC(density, ndiv);
            this.init(density, ndiv);
        }
    }

    /**
     * このVolumetricDataの状態を保存しておく.
     * 
     */
    public void save() {
        density_buff = new float[density.length];
        for (int i = 0; i < density.length; i++)
            density_buff[i] = density[i];

        ndiv_buff = new int[3];
        for (int i = 0; i < ndiv.length; i++)
            ndiv_buff[i] = ndiv[i];

        delta_buff = new float[3][3];
        for (int i = 0; i < delta.length; i++)
            for (int j = 0; j < delta[i].length; j++)
                delta_buff[i][j] = delta[i][j];

    }

    /**
     * save()メソッドを呼ぶ直前の状態に戻す.
     */
    public void restore() {
        this.density = this.density_buff;
        this.ndiv = this.ndiv_buff;
        this.delta = this.delta_buff;
        // init(density, ndiv);

        logger.debug("restored data ");
        logger.debug("ndiv: " + ndiv[0] + ", " + ndiv[1] + ", " + ndiv[2]);
        for (int i = 0; i < 3; i++)
            logger.debug("delta " + i + ": " + delta[i][0] + ", " + delta[i][1]
                    + ", " + delta[i][2]);
    }

    /**
     * 周期境界条件より付加している端のデータを削除する.
     * 
     */
    public void stripPBC() {
        if (!pbc)
            return;
        int[] ndiv_tmp = new int[3];
        for (int i = 0; i < 3; i++)
            ndiv_tmp[i] = ndiv[i] - 1;
        int Ntot_tmp = ndiv_tmp[0] * ndiv_tmp[1] * ndiv_tmp[2];
        float[] density_tmp = new float[Ntot_tmp];
        for (int i = 0; i < ndiv_tmp[0]; i++) {
            for (int j = 0; j < ndiv_tmp[1]; j++) {
                for (int k = 0; k < ndiv_tmp[2]; k++) {
                    int index = i * ndiv[1] * ndiv[2] + j * ndiv[2] + k;
                    int newind = i * ndiv_tmp[1] * ndiv_tmp[2] + j
                            * ndiv_tmp[2] + k;
                    density_tmp[newind] = this.density[index];
                }
            }
        }
        this.density = density_tmp;
        this.ndiv = ndiv_tmp;
        this.Ntot = Ntot_tmp;
    }

    public void doPBC(float[] den, int[] ndi) {
        this.ndiv = new int[3];
        ndiv[0] = ndi[0] + 1;
        ndiv[1] = ndi[1] + 1;
        ndiv[2] = ndi[2] + 1;
        Ntot = ndiv[0] * ndiv[1] * ndiv[2];
        this.density = new float[Ntot];
        for (int i = 0; i < ndi[0]; i++) {
            for (int j = 0; j < ndi[1]; j++) {
                for (int k = 0; k < ndi[2]; k++) {
                    int index = i * ndi[1] * ndi[2] + j * ndi[2] + k;
                    int newind = i * ndiv[1] * ndiv[2] + j * ndiv[2] + k;
                    density[newind] = den[index];
                }
            }
        }
        for (int i = 0; i < ndiv[0]; i++) {
            for (int j = 0; j < ndiv[1]; j++) {
                int ind = i * ndiv[1] * ndiv[2] + j * ndiv[2];
                density[ind + ndiv[2] - 1] = density[ind];
            }
        }
        for (int i = 0; i < ndiv[0]; i++) {
            for (int j = 0; j < ndiv[2]; j++) {
                int ind = i * ndiv[1] * ndiv[2] + j;
                density[ind + ndiv[2] * (ndiv[1] - 1)] = density[ind];
            }
        }
        for (int i = 0; i < ndiv[1]; i++) {
            for (int j = 0; j < ndiv[2]; j++) {
                int ind = i * ndiv[2] + j;
                density[ndiv[2] * ndiv[1] * (ndiv[0] - 1) + ind] = density[ind];
            }
        }
    }

    private float[] doPBC(float[] den) {
        int[] ndi = new int[3];
        for (int i = 0; i < ndi.length; i++)
            ndi[i] = ndiv[i] - 1;

        float[] tmpden = new float[Ntot];

        for (int i = 0; i < ndi[0]; i++) {
            for (int j = 0; j < ndi[1]; j++) {
                for (int k = 0; k < ndi[2]; k++) {
                    int index = i * ndi[1] * ndi[2] + j * ndi[2] + k;
                    int newind = i * ndiv[1] * ndiv[2] + j * ndiv[2] + k;
                    tmpden[newind] = den[index];
                }
            }
        }
        for (int i = 0; i < ndiv[0]; i++) {
            for (int j = 0; j < ndiv[1]; j++) {
                int ind = i * ndiv[1] * ndiv[2] + j * ndiv[2];
                tmpden[ind + ndiv[2] - 1] = tmpden[ind];
            }
        }
        for (int i = 0; i < ndiv[0]; i++) {
            for (int j = 0; j < ndiv[2]; j++) {
                int ind = i * ndiv[1] * ndiv[2] + j;
                tmpden[ind + ndiv[2] * (ndiv[1] - 1)] = tmpden[ind];
            }
        }
        for (int i = 0; i < ndiv[1]; i++) {
            for (int j = 0; j < ndiv[2]; j++) {
                int ind = i * ndiv[2] + j;
                tmpden[ndiv[2] * ndiv[1] * (ndiv[0] - 1) + ind] = tmpden[ind];
            }
        }
        return tmpden;
    }

    private int nshift0 = 0;

    private int nshift1 = 0;

    private int nshift2 = 0;

    /**
     * このデータを"並進"を設定する. 実際にデータが並進されるのは各種getterを呼んだタイミングで行われる
     * 
     * @param n1
     *            一つ目の方向にどれくらいずらすか指定
     * @param n2
     *            二つ目の方向にどれくらいずらすか指定
     * @param n3
     *            三つ目の方向にどれくらいずらすか指定
     */
    public void setShift(int n1, int n2, int n3) {
        if (n1 > ndiv[0] || n2 > ndiv[1] || n3 > ndiv[2]) {
            logger.error("can't shift that much!");
            return;
        }
        this.nshift0 = n1;
        this.nshift1 = n2;
        this.nshift2 = n3;
        // nshift0 += n1;
        // nshift1 += n2;
        // nshift2 += n3;
        // if ( nshift0 < 0 ) nshift0 = 0;
        // if ( nshift0 >= ndiv[0] ) nshift0=ndiv[0]-1;
        // if ( nshift1 < 0 ) nshift1 = 0;
        // if ( nshift1 >= ndiv[1] ) nshift1=ndiv[1]-1;
        // if ( nshift2 < 0 ) nshift2 = 0;
        // if ( nshift2 >= ndiv[2] ) nshift2=ndiv[2]-1;
    }

    /**
     * このデータを"並進"を設定する. ただし, セルが(実際に)どれくらいシフトしたか, によって 決める.
     * 実際にデータが並進されるのは各種getterを呼んだタイミングで行われる
     * 
     * @param d1
     *            一つ目の方向にどれくらいずらすか指定
     * @param d2
     *            二つ目の方向にどれくらいずらすか指定
     * @param d3
     *            三つ目の方向にどれくらいずらすか指定
     */
    public void setShiftOld(double d1, double d2, double d3) {
        double cell0 = (double) (delta[0][0] * ndiv[0] + delta[0][1] * ndiv[1] + delta[0][2]
                * ndiv[2]);
        double cell1 = (double) (delta[1][0] * ndiv[0] + delta[1][1] * ndiv[1] + delta[1][2]
                * ndiv[2]);
        double cell2 = (double) (delta[2][0] * ndiv[0] + delta[2][1] * ndiv[1] + delta[2][2]
                * ndiv[2]);
        double[][] cell = new double[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                cell[i][j] = delta[i][j] * ndiv[j];
            }
        }
        double[][] shift = new double[3][1];
        shift[0][0] = d1;
        shift[1][0] = d2;
        shift[2][0] = d3;
        Matrix A = new Matrix(cell);
        Matrix B = new Matrix(shift);
        Matrix N = A.solve(B);
        nshift0 = (int) Math.round(N.get(0, 0) * ndiv[0]);
        nshift1 = (int) Math.round(N.get(1, 0) * ndiv[1]);
        nshift2 = (int) Math.round(N.get(2, 0) * ndiv[2]);
        logger.debug("set shift : " + nshift0 + ", " + nshift1 + ", " + nshift2
                + " from float: " + d1 + ", " + d2 + ", " + d3);
    }

    public void setShift(double d1, double d2, double d3) {
        logger.debug("d1, d2, d3: " + d1 + ", " + d2 + ", " + d3);
        nshift0 = (int) Math.round(d1 * ndiv[0]);
        nshift1 = (int) Math.round(d2 * ndiv[1]);
        nshift2 = (int) Math.round(d3 * ndiv[2]);
        logger.debug("nshift: " + nshift0 + ", " + nshift1 + ", " + nshift2);
    }

    public void setDiffShift(double d1, double d2, double d3) {
        nshift0 += (int) Math.round(d1 * ndiv[0]);
        nshift1 += (int) Math.round(d2 * ndiv[1]);
        nshift2 += (int) Math.round(d3 * ndiv[2]);
    }

    /**
     * このデータを"拡大"する.
     * 
     * @param n1
     *            一つ目の方向へ何倍するか指定
     * @param n2
     *            二つ目の方向に何倍するか指定
     * @param n3
     *            三つ目の方向に何倍するか指定
     */
    public void enlarge(int n1, int n2, int n3) {
        try {
            int NtotOld = ndiv_buff[0] * ndiv_buff[1] * ndiv_buff[2];
            ndiv = new int[ndiv.length];
            ndiv[0] = ndiv_buff[0] * n1;
            ndiv[1] = ndiv_buff[1] * n2;
            ndiv[2] = ndiv_buff[2] * n3;
            Ntot = ndiv[0] * ndiv[1] * ndiv[2];

            float[] density_tmp = new float[ndiv[2] * ndiv_buff[1]
                    * ndiv_buff[0]];
            int ntmp1 = 0;
            int ntmp2 = 0;
            int ntmp3 = 0;
            int nbase = 0;
            for (int i = 0; i < ndiv_buff[0] * ndiv_buff[1]; i++) {
                for (int j = 0; j < ndiv_buff[2]; j++) {
                    int nold = i * ndiv_buff[2] + j;
                    for (int k = 0; k < n3; k++) {
                        int tmp = k * ndiv_buff[2];
                        int nnew = i * ndiv[2] + j + tmp;
                        density_tmp[nnew] = density_buff[nold];
                    }
                }
            }

            float[] density_tmp2 = new float[ndiv[2] * ndiv[1] * ndiv_buff[0]];
            for (int i = 0; i < ndiv_buff[0]; i++) {
                for (int j = 0; j < n3 * ndiv_buff[2] * ndiv_buff[1]; j++) {
                    int nold = i * n3 * ndiv_buff[2] * ndiv_buff[1] + j;
                    for (int k = 0; k < n2; k++) {
                        int tmp = k * ndiv_buff[1] * ndiv[2];
                        int nnew = i * ndiv[1] * ndiv[2] + j + tmp;
                        density_tmp2[nnew] = density_tmp[nold];
                    }
                }
            }

            if (!pbc) {
                density = new float[Ntot];
                for (int i = 0; i < density_tmp2.length; i++) {
                    for (int j = 0; j < n1; j++) {
                        int nnew = j * density_tmp2.length + i;
                        density[nnew] = density_tmp2[i];
                    }
                }
            } else {
                float[] den = new float[Ntot];
                int[] ndi = new int[3];
                ndi[0] = ndiv[0];
                ndi[1] = ndiv[1];
                ndi[2] = ndiv[2];
                for (int i = 0; i < density_tmp2.length; i++) {
                    for (int j = 0; j < n1; j++) {
                        int nnew = j * density_tmp2.length + i;
                        den[nnew] = density_tmp2[i];
                    }
                }
                doPBC(den, ndi);
            }
        } catch (Exception exc) {
            exc.printStackTrace();
        }
    }

    /**
     * 補間の形式をセット. デフォルトはLOG.
     * 
     * @param 補間の形式
     */
    public void setInterpolationScheme(int interpolationScheme) {
        this.interpolationScheme = interpolationScheme;
    }

    /**
     * 指定の値よりも大きな値を持つ, 一番最初のインデックスを返す.
     * 
     * @param value
     *            キーとなる値
     * @return valueよりも大きな値を持つ, 一番最初のインデックス.
     */
    public int binarySearch(float value) {
        java.util.Arrays.sort(sortedDensity);
        int ret = java.util.Arrays.binarySearch(sortedDensity, value);
        if (ret < 0) {
            ret *= -1;
        }
        return ret;
    }

    /**
     * 密度の配列を返す. ただし, 拡大した場合においても拡大前の密度の情報を返す.
     * 
     * @return 密度の配列; ただし拡大している場合拡大前の密度.
     */
    public float[] getOrigDensity() {
        return density_buff;
    }

    /**
     * 原点を取得
     * 
     * @return このデータの原点
     */
    public float[] getOrigin() {
        return origin;
    }

    // private void init() {
    // if ( density != null ) {
    // for ( int i=0 ; i<density.length ; i++ ) {
    // if ( density[i] == meaninglessValue ) {
    // continue;
    // }
    // if ( density[i] < minVal ) {
    // minVal = density[i];
    // } else if ( density[i] > maxVal ) {
    // maxVal = density[i];
    // }
    // }
    // sortedDensity = new float[density.length];
    // density_buff = new float[density.length];
    // for ( int i=0 ; i<sortedDensity.length ; i++ ) {
    // sortedDensity[i] = density[i];
    // density_buff[i] = density[i];
    // }
    // midVal = (float)
    // BaseColorMap.getValueFromLogScaledValue((double)minVal,(double)maxVal,0.5)
    // ;
    // }
    // if ( ndiv != null ) {
    // ndiv_buff = new int[ndiv.length];
    // for ( int i=0 ; i<ndiv_buff.length ; i++ ) {
    // ndiv_buff[i] = ndiv[i];
    // }
    // }
    // }

    public void createBuffer(float[] den, int[] ndi) {
        init(den, ndi);
    }

    private void init(float[] den, int[] ndi) {
        if (den != null) {
            for (int i = 0; i < den.length; i++) {
                if (den[i] == meaninglessValue) {
                    continue;
                }
                if (den[i] < minVal) {
                    minVal = den[i];
                } else if (den[i] > maxVal) {
                    maxVal = den[i];
                }
            }
            sortedDensity = new float[den.length];
            density_buff = new float[den.length];
            for (int i = 0; i < sortedDensity.length; i++) {
                sortedDensity[i] = den[i];
                density_buff[i] = den[i];
            }
            if (interpolationScheme == LOG) {
                midVal = (float) BaseColorMap.getValueFromLogScaledValue(
                        (double) minVal, (double) maxVal, 0.5);
            } else if (interpolationScheme == LINEAR) {
                midVal = (float) BaseColorMap.getValueFromLinearScaledValue(
                        (double) minVal, (double) maxVal, 0.5);
            }
        }
        if (ndiv != null) {
            ndiv_buff = new int[ndi.length];
            for (int i = 0; i < ndiv_buff.length; i++) {
                ndiv_buff[i] = ndi[i];
            }
        }
        if (delta != null) {
            delta_buff = new float[3][3];
            for (int i = 0; i < delta.length; i++)
                for (int j = 0; j < delta[i].length; j++)
                    delta_buff[i][j] = delta[i][j];
        }
    }

    /**
     * 密度の最小値を取得.
     * 
     * @return 密度の最小値
     */
    public float getMinVal() {
        return minVal;
    }

    /**
     * 密度の最大値を取得.
     * 
     * @return 密度の最大値
     */
    public float getMaxVal() {
        return maxVal;
    }

    /**
     * 格子点の分割数を返す
     * 
     * @return 格子点の分割数配列
     */
    public int[] getNumDiv() {
        return ndiv;
    }

    /**
     * 全格子の数を返す. 変なデータの場合は-1を返す.
     * 
     * @return 格子点の数. 「変なデータ」の場合-1.
     */
    public int getNtot() {
        if (ndiv == null || ndiv.length < 3) {
            return -1;
        }
        return ndiv[0] * ndiv[1] * ndiv[2];
    }

    /**
     * 密度の中点?を取得
     * 
     * @return 密度の中点
     */
    public float getMidVal() {
        return midVal;
    }

    private float[] getShiftedDensity() {
        stripPBC();

        float[] den_tmp = new float[density.length];
        for (int i = 0; i < den_tmp.length; i++)
            den_tmp[i] = density[i];

        for (int i = 0; i < ndiv[2]; i++) {
            for (int j = 0; j < ndiv[1]; j++) {
                for (int k = 0; k < ndiv[0]; k++) {
                    int itmp = i + nshift2;
                    int jtmp = j + nshift1;
                    int ktmp = k + nshift0;
                    if (itmp >= ndiv[2]) {
                        itmp -= ndiv[2];
                    } else if (itmp < 0) {
                        itmp += ndiv[2];
                    }
                    if (jtmp >= ndiv[1]) {
                        jtmp -= ndiv[1];
                    } else if (jtmp < 0) {
                        jtmp += ndiv[1];
                    }
                    if (ktmp >= ndiv[0]) {
                        ktmp -= ndiv[0];
                    } else if (ktmp < 0) {
                        ktmp += ndiv[0];
                    }

                    int ind1 = i + j * ndiv[2] + k * ndiv[2] * ndiv[1];
                    int ind2 = itmp + jtmp * ndiv[2] + ktmp * ndiv[2] * ndiv[1];
                    den_tmp[ind2] = density[ind1];
                }
            }
        }
        doPBC(this.density, this.ndiv);
        den_tmp = doPBC(den_tmp);
        return den_tmp;
    }

    /**
     * 空間上の各点における電荷密度の値の配列を返す.
     * 
     * @return 電荷密度の値を並べた配列.
     */
    public float[] getDensity() {
        if (nshift0 == 0 && nshift1 == 0 && nshift2 == 0) {
            return this.density;
        } else {
            return getShiftedDensity();
        }
    }

    /**
     * i,j,k -> coordinateみたいな配列を返す
     * 
     * @return 第一次元:a軸のインデックス, 第二次元:b軸のインデックス, 第三次元:c軸のインデックス, 第四次元: x,y,zという配列.
     */
    public float[][][][] getCoordinates_explicitIndeces() {
        float[][][][] ret = new float[ndiv[2]][ndiv[1]][ndiv[0]][3];
        for (int i = 0; i < ndiv[2]; i++) {
            for (int j = 0; j < ndiv[1]; j++) {
                for (int k = 0; k < ndiv[0]; k++) {
                    ret[i][j][k][0] = k * delta[0][0] + j * delta[1][0] + i
                            * delta[2][0];
                    ret[i][j][k][1] = k * delta[0][1] + j * delta[1][1] + i
                            * delta[2][1];
                    ret[i][j][k][2] = k * delta[0][2] + j * delta[1][2] + i
                            * delta[2][2];
                }
            }
        }
        return ret;
    }

    /**
     * i,j,k -> coordinateみたいな配列を返す スケールも行う.
     * 
     * @return 第一次元:a軸のインデックス, 第二次元:b軸のインデックス, 第三次元:c軸のインデックス, 第四次元: x,y,zという配列.
     */
    public float[][][][] getCoordinates_explicitIndeces(float lenmax,
            float[] jusin) {
        float[][][][] ret = new float[ndiv[2]][ndiv[1]][ndiv[0]][3];
        for (int i = 0; i < ndiv[2]; i++) {
            for (int j = 0; j < ndiv[1]; j++) {
                for (int k = 0; k < ndiv[0]; k++) {
                    ret[i][j][k][0] = (k * delta[0][0] + j * delta[1][0] + i
                            * delta[2][0] - jusin[0])
                            / lenmax;
                    ret[i][j][k][1] = (k * delta[0][1] + j * delta[1][1] + i
                            * delta[2][1] - jusin[1])
                            / lenmax;
                    ret[i][j][k][2] = (k * delta[0][2] + j * delta[1][2] + i
                            * delta[2][2] - jusin[2])
                            / lenmax;
                }
            }
        }
        return ret;
    }

    /**
     * i,j,k->値, みたいな配列を返す.
     * 
     * @return 第一次元:a軸のインデックス, 第二次元:b軸のインデックス, 第三次元:c軸のインデックスという配列を返す
     */
    public float[][][] getValue_explicitIndeces() {
        float ret[][][] = new float[ndiv[2]][ndiv[1]][ndiv[0]];
        if (nshift0 == 0 && nshift1 == 0 && nshift2 == 0) {
            for (int i = 0; i < ndiv[2]; i++) {
                for (int j = 0; j < ndiv[1]; j++) {
                    for (int k = 0; k < ndiv[0]; k++) {
                        int n = i + j * ndiv[2] + k * ndiv[2] * ndiv[1];
                        ret[i][j][k] = density[n];
                    }
                }
            }
        } else {
            float[] den = getShiftedDensity();
            for (int i = 0; i < ndiv[2]; i++) {
                for (int j = 0; j < ndiv[1]; j++) {
                    for (int k = 0; k < ndiv[0]; k++) {
                        int n = i + j * ndiv[2] + k * ndiv[2] * ndiv[1];
                        ret[i][j][k] = den[n];
                    }
                }
            }
        }
        return ret;
    }

    /**
     * 「格子点の座標」を返す
     * 
     * @return 格子点の座標; 第一次元がx,y,zのラベル, 第二次元が対応する格子点配列.
     */
    public float[][] getCoordinates() {
        float[][] ret = new float[3][Ntot];
        for (int i = 0; i < ndiv[2]; i++) {
            for (int j = 0; j < ndiv[1]; j++) {
                for (int k = 0; k < ndiv[0]; k++) {
                    int n = i + j * ndiv[2] + k * ndiv[2] * ndiv[1];
                    ret[0][n] = k * delta[0][0] + j * delta[1][0] + i
                            * delta[2][0];
                    ret[1][n] = k * delta[0][1] + j * delta[1][1] + i
                            * delta[2][1];
                    ret[2][n] = k * delta[0][2] + j * delta[1][2] + i
                            * delta[2][2];
                }
            }
        }
        // Vector tmp = new Vector();
        // for ( int i=0 ; i<Ntot ; i++ ) {
        // tmp.addElement(getCoordinatesFrom(i));
        // }
        // if ( tmp.size() == 0 ) {
        // return null;
        // }
        // float [][] ret = new float[3][tmp.size()];
        // for ( int i=0 ; i<tmp.size() ; i++ ) {
        // float [] foo = (float[]) tmp.elementAt(i);
        // for ( int j=0 ; j<3 ; j++ ) {
        // ret[j][i] = foo[j];
        // }
        // }
        return ret;
    }

    /**
     * 一次元で表されたインデックスより, 三次元のインデックスを返す.
     * 
     * @param index1d
     *            一次元で表されたインデックス.
     * @return 三次元で表されたインデックス.
     */
    private int[] get3DIndexFrom1DIndex(int index1d) {
        int nz, ny, nx;
        nz = index1d % ndiv[2];
        ny = ((index1d - nz) / ndiv[2]) % ndiv[1];
        nx = ((index1d - nz) / ndiv[2] - ny) / ndiv[1];
        return new int[] { nx, ny, nz };
    }

    /**
     * 三次元のインデックス配列より座標を計算.
     * 
     * @param index
     *            三次元のインデックス
     */
    private float[] getCoordinatesFrom(int[] index) {
        float[] ret = new float[3];
        for (int i = 0; i < 3; i++) {
            ret[i] = 0.f;
            for (int j = 0; j < 3; j++) {
                ret[i] += index[j] * delta[j][i];
            }
            ret[i] += origin[i];
        }
        return ret;
    }

    /**
     * 一次元のインデックス配列より座標を計算.
     * 
     * @param index
     *            一次元のインデックス
     */
    private float[] getCoordinatesFrom(int index) {
        return getCoordinatesFrom(get3DIndexFrom1DIndex(index));
    }

    /**
     * 補間の形式を取得. デフォルトはLOG.
     * 
     * @return 補間の形式
     */
    public int getInterpolationScheme() {
        return interpolationScheme;
    }

    public String toString() {
        String ret = "min: " + String.valueOf(minVal) + ", max: "
                + String.valueOf(maxVal);
        return ret;
    }

    /**
     * a軸, b軸, c軸の変化分を取得
     * 
     * @return 上記の通り
     */
    public float[][] getDelta() {
        return this.delta;
    }

    /**
     * a軸, b軸, c軸の変化分を設定.
     * 
     * @param delta
     *            新しい変化分
     */
    public void setDelta(float[][] delta) {
        this.delta = delta;
    }

    public void setDensity(float[] density) {
        this.density = density;
    }

}
