#ifndef _GRID3D_H_
#define _GRID3D_H_

#include <stdlib.h>
#include <iostream>
#include <fstream>

/// written by Christoph Federrath, Sep 2010
/// based on a 2D class written by Philipp Girichidis

class GRID3D
{
  public:
    double * field;
    bool * is_set;

  private:
    long   dimx, dimy, dimz, ntot;
    double xl, xh, yl, yh, zl, zh;
    double Lx, Ly, Lz;
    double minval, maxval;
    double cell_dx, cell_dy, cell_dz;
    static const bool Debug = false;

  //////////////////////////////////////////////////////////////
  // CONSTRUCTOR & SET GLOBAL PROPERTIES
  //////////////////////////////////////////////////////////////

  public: GRID3D()
  {
      field = 0;
      is_set = 0;
      minval = 0.0;
      maxval = 0.0;
      dimx = 0;
      dimy = 0;
      dimz = 0;
      ntot = 0;
      if (Debug) std::cout<<"GRID3D: default constructor called."<<std::endl;
  };
  public: GRID3D(int N)
  {
      ntot = N*N*N;
      field = new double[ntot];
      is_set = new bool[ntot];
      minval = 0.0;
      maxval = 0.0;
      dimx = N;
      dimy = N;
      dimz = N;
      clear();
      if (Debug) std::cout<<"GRID3D: constructor(N) called."<<std::endl;
  };
  public: GRID3D(int Nx, int Ny, int Nz)
  {
      ntot = Nx*Ny*Nz;
      field = new double[ntot];
      is_set = new bool[ntot];
      minval = 0.0;
      maxval = 0.0;
      dimx = Nx;
      dimy = Ny;
      dimz = Nz;
      clear();
      if (Debug) std::cout<<"GRID3D: constructor(Nx,Ny,Nz) called."<<std::endl;
  };

  /// copy constructor (to make a deep copy of a given GRID3D)
  GRID3D(const GRID3D &t)
  {
    ntot = t.ntot;
    dimx = t.dimx;
    dimy = t.dimy;
    dimz = t.dimz;
    minval = t.minval;
    maxval = t.maxval;
    if (t.field != 0 && t.is_set != 0)
    {
	field = new double[ntot];
	is_set = new bool[ntot];
	for (long n=0; n<ntot; n++)
	{
	  field[n] = t.field[n];
	  is_set[n] = t.is_set[n];
	}
    }
    else { field = 0; is_set = 0; }
  };
  
  /// assigment operator
  GRID3D &operator=(const GRID3D &t)
  {
    if (Debug) std::cout<<"GRID3D: operator= called."<<std::endl;
    if (&t == this) return *this;
    else
    {
      ntot = t.ntot;
      dimx = t.dimx;
      dimy = t.dimy;
      dimz = t.dimz;
      minval = t.minval;
      maxval = t.maxval;
      if (field != 0 && is_set != 0) // delete old fields
      {
	delete [] field; delete [] is_set; field = 0; is_set = 0;
      }
      if (t.field != 0 && t.is_set != 0) // assign
      {
	field = new double[ntot];
	is_set = new bool[ntot];
	for (long n=0; n<ntot; n++)
	{
	  field[n] = t.field[n];
	  is_set[n] = t.is_set[n];
	}
      }
      return *this;
    }
  };
  
  public: ~GRID3D()
  {
    if (field != 0) { delete [] field; field = 0; if (Debug) std::cout<<"GRID3D: field destroyed."<<std::endl;}
    if (is_set != 0) { delete [] is_set; is_set = 0; if (Debug) std::cout<<"GRID3D: is_set destroyed."<<std::endl;}
    if (Debug) std::cout<<"GRID3D: destructor called."<<std::endl;
  };

  public: inline void set_bnds(double low, double high)
  {
    set_bnds(low, high, low, high, low, high);
  }
  public: inline void set_bnds(double xlow, double xhigh, double ylow, double yhigh, double zlow, double zhigh)
  {
    if(xhigh-xlow <= 0.0)
    {
	std::cout << " xlow > xhigh !! swapping these values !!" << std::endl;
	double temp = xlow; xlow = xhigh; xhigh = temp;
    }
    if(yhigh-ylow <= 0.0)
    {
	std::cout << " ylow > yhigh !! swapping these values !!" << std::endl;
	double temp = ylow; ylow = yhigh; yhigh = temp;
    }
    if(zhigh-zlow <= 0.0)
    {
	std::cout << " zlow > zhigh !! swapping these values !!" << std::endl;
	double temp = zlow; zlow = zhigh; zhigh = temp;
    }
    xl = xlow; xh = xhigh;
    yl = ylow; yh = yhigh;
    zl = zlow; zh = zhigh;
    Lx = xh-xl;
    Ly = yh-yl;
    Lz = zh-zl;
    cell_dx = (xhigh-xlow)/(double)dimx;
    cell_dy = (yhigh-ylow)/(double)dimy;
    cell_dz = (zhigh-zlow)/(double)dimz;
  }


  //////////////////////////////////////////////////////////////
  // GLOBAL OBJECT INFORMATION
  //////////////////////////////////////////////////////////////

  public: inline int get_dimx() { return dimx; }
  public: inline int get_dimy() { return dimy; }
  public: inline int get_dimz() { return dimz; }
  public: inline long get_ntot() { return ntot; }

  public: inline double get_xmin() { return xl; }
  public: inline double get_xmax() { return xh; }
  public: inline double get_ymin() { return yl; }
  public: inline double get_ymax() { return yh; }
  public: inline double get_zmin() { return zl; }
  public: inline double get_zmax() { return zh; }

  public: inline double get_dx() { return cell_dx; }
  public: inline double get_dy() { return cell_dy; }
  public: inline double get_dz() { return cell_dz; }

  public: inline double get_Lx() { return Lx; }
  public: inline double get_Ly() { return Ly; }
  public: inline double get_Lz() { return Lz; }

  public: inline double get_dV() { return cell_dx*cell_dy*cell_dz; }

  public: inline void print_bnds(std::ostream & out = std::cout)
  {
	out << " grid bounds"
    	<< " [" << xl << ":" << xh << "]"
		<< " [" << yl << ":" << yh << "]"
    	<< " [" << zl << ":" << zh << "]" << std::endl;
  }
  public: inline void print_minmax(std::ostream & out = std::cout)
  {
    double MIN, MAX;
    minmax(MIN, MAX);
    out << " min value : " << minval << std::endl;
    out << " max value : " << maxval << std::endl;
  }
  public: inline void print_sum(std::ostream & out = std::cout)
  {
    double val = sum();
    out << " sum    : " << val << std::endl;
  }
  public: inline void print_sum_vw(std::ostream & out = std::cout)
  {
    double val_vw = sum_vw();
    out << " sum_vw : " << val_vw << std::endl;
  }
  
  public: double first_set()
  {
    // find first data point that is set and use it as starting point for min and max
    long counter = 0;
    while(counter < ntot)
    {
	if (is_set[counter])
	    return field[counter];
	counter++;
    }
    // no field is set so far!!
    std::cout << " NO DATA POINT IS SET IN YOUR FIELD, ABORT!!" << std::endl;
    exit(1);
    return 0.0;
  }

  public: double min()
  {
    double MIN = first_set();
    for (long n=0; n<ntot; n++)
	if (is_set[n] && (field[n] < MIN)) MIN = field[n];
    minval = MIN;
    return minval;
  }
  public: double max()
  {
    double MAX = first_set();
    for (long n=0; n<ntot; n++)
	if (is_set[n] && (field[n] > MAX)) MAX = field[n];
    maxval = MAX;
    return maxval;
  }
  public: void minmax(double &tmin, double &tmax)
  {
    tmin = first_set();
    tmax = tmin;
    for (long n=0; n<ntot; n++)
    {
	  if(is_set[n] && (field[n] < tmin)) tmin = field[n];
	  if(is_set[n] && (field[n] > tmax)) tmax = field[n];
    }
    minval = tmin;
    maxval = tmax;
  }
  public: double sum()
  {
    double ret = 0.;
    for (long n=0; n<ntot; n++)
    {
	  if(is_set[n]) ret += field[n];
    }
    return ret;
  }
  public: double sum_vw()
  {
    double ret = sum()*Lx*Ly*Lz/(double)dimx/(double)dimy/(double)dimz;
    return ret;
  }


  //////////////////////////////////////////////////////////////
  // GET SINGLE FIELD INFORMATION
  //////////////////////////////////////////////////////////////

  public: inline double getX(int i)
  {
    return xl + ((double)i+0.5) * Lx/(double)dimx;
  }
  public: inline double getY(int i)
  {
    return yl + ((double)i+0.5) * Ly/(double)dimy;
  }
  public: inline double getZ(int i)
  {
    return zl + ((double)i+0.5) * Lz/(double)dimz;
  }
  public: inline double get_val(int ix, int iy, int iz)
  {
    if((ix>=0) && (ix<dimx) && (iy>=0) && (iy<dimy) && (iz>=0) && (iz<dimz))
    {
      long n = index(ix,iy,iz);
      return field[n];
    }
    else return 0.0;
  }
  /// returns the index of the cell in which (x,y,z) is located
  private: inline void get_bins(double x, double y, double z, int &ix, int &iy, int &iz)
  {
    double xmxl = x-xl;
    double ymyl = y-yl;
    double zmzl = z-zl;
    if (xmxl >= 0) ix = (int)(xmxl/Lx*(double)dimx);
    else           ix = (int)(xmxl/Lx*(double)dimx)-1;
    if (ymyl >= 0) iy = (int)(ymyl/Ly*(double)dimy);
    else           iy = (int)(ymyl/Ly*(double)dimy)-1;
    if (zmzl >= 0) iz = (int)(zmzl/Lz*(double)dimz);
    else           iz = (int)(zmzl/Lz*(double)dimz)-1;
  }
  private: inline long index(int ix, int iy, int iz)
  {
    return iz*dimx*dimy + iy*dimx + ix;
  }
  private: inline double cell_center_x(int i)
  {
    return xl + ((double)i + 0.5)*cell_dx;
  }
  private: inline double cell_center_y(int i)
  {
    return yl + ((double)i + 0.5)*cell_dy;
  }
  private: inline double cell_center_z(int i)
  {
    return zl + ((double)i + 0.5)*cell_dz;
  }
  private: inline double get_cell_xlow(int ix)
  {
    return xl + ((double)ix)*cell_dx;
  }
  private: inline double get_cell_xhigh(int ix)
  {
    return xl + ((double)(ix+1))*cell_dx;
  }
  private: inline double get_cell_ylow(int iy)
  {
    return yl + ((double)iy)*cell_dy;
  }
  private: inline double get_cell_yhigh(int iy)
  {
    return yl + ((double)(iy+1))*cell_dy;
  }
  private: inline double get_cell_zlow(int iz)
  {
    return zl + ((double)iz)*cell_dz;
  }
  private: inline double get_cell_zhigh(int iz)
  {
    return zl + ((double)(iz+1))*cell_dz;
  }
  private: inline void cell_corner_x(int i, double &xlow, double &xhigh)
  {
    xlow  = xl + ((double)i)*cell_dx;
    xhigh = xl + ((double)(i+1))*cell_dx;
  }
  private: inline void cell_corner_y(int i, double &ylow, double &yhigh)
  {
    ylow  = yl + ((double)i)*cell_dy;
    yhigh = yl + ((double)(i+1))*cell_dy;
  }
  private: inline void cell_corner_z(int i, double &zlow, double &zhigh)
  {
    zlow  = zl + ((double)i)*cell_dz;
    zhigh = zl + ((double)(i+1))*cell_dz;
  }
  private: inline void cell_corners(int ix, int iy, int iz,
				    double &xlow, double &xhigh,
				    double &ylow, double &yhigh,
				    double &zlow, double &zhigh)
  {
    xlow  = xl + ((double)ix)*cell_dx;
    xhigh = xl + ((double)(ix+1))*cell_dx;
    ylow  = yl + ((double)iy)*cell_dy;
    yhigh = yl + ((double)(iy+1))*cell_dy;
    zlow  = zl + ((double)iz)*cell_dz;
    zhigh = zl + ((double)(iz+1))*cell_dz;
  }


  /// this subroutine adds cell-centered data values given at (x,y,z) to the appropriate cell of the GRID3D
  /// taking into account all overlapping fractions of the original data
  /// This routine does not use interpolation (it is thus conserving the integral over all datavalues)
  public: inline void add_coord_fields(double x, double y, double z, double dx, double dy, double dz, double val)
  {
    // get cell at the center x,y
    double xlow  = x-0.5*dx;
    double xhigh = x+0.5*dx;
    double ylow  = y-0.5*dy;
    double yhigh = y+0.5*dy;
    double zlow  = z-0.5*dz;
    double zhigh = z+0.5*dz;
    double cell_vol = cell_dx*cell_dy*cell_dz;
    double data_vol = dx*dy*dz;
    int ix0, ix1, iy0, iy1, iz0, iz1;
    get_bins(xlow, ylow, zlow, ix0, iy0, iz0);
    get_bins(xhigh, yhigh, zhigh, ix1, iy1, iz1);

    // check whether data cell fits entirely in grid cell
    // equal to ix0==ix1 and iy0==iy1 and iz0==iz1
    if((ix0==ix1) && (iy0==iy1) && (iz0==iz1))
    {
	if (Debug) std::cout << "data cell fits entirely in grid cell" << std::endl;
	add(ix0, iy0, iz0, val/cell_vol);
    }
    else
    {
	// get specific value for one cell
	// N_cells for dx*dy*dz = dx*dy*dz/(cell_dx*cell_dy*cell_dz)
	double spec_val = val/data_vol;

	// data cell occupies more than one cell in x direction
	if((ix0!=ix1) && (iy0==iy1) && (iz0==iz1))
	{
	    if (Debug) std::cout << "data cell occupies more than one cell in x direction" << std::endl;
	    double frac;
	    // left cell
	    frac = (get_cell_xhigh(ix0)-xlow)*dy*dz/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // middle cells
	    for(int xx=ix0+1; xx<ix1; xx++)
	    {
		frac = cell_dx*dy*dz/cell_vol;
		add(xx,iy0,iz0,frac*spec_val);
	    }
	    // right cell
	    frac = (xhigh-get_cell_xlow(ix1))*dy*dz/cell_vol;
	    add(ix1,iy0,iz0,frac*spec_val);
	}

	// data cell occupies more than one cell in y direction
	if((ix0==ix1) && (iy0!=iy1) && (iz0==iz1))
	{
	    if (Debug) std::cout << "data cell occupies more than one cell in y direction" << std::endl;
	    double frac;
	    // bottom cell
	    frac = dx*(get_cell_yhigh(iy0)-ylow)*dz/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // middle cells
	    for (int yy=iy0+1; yy<iy1; yy++)
	    {
		frac = dx*cell_dy*dz/cell_vol;
		add(ix0,yy,iz0,frac*spec_val);
	    }
	    // top cell
	    frac = dx*(yhigh-get_cell_ylow(iy1))*dz/cell_vol;
	    add(ix0,iy1,iz0,frac*spec_val);
	}

	// data cell occupies more than one cell in z direction
	if((ix0==ix1) && (iy0==iy1) && (iz0!=iz1))
	{
	    if (Debug) std::cout << "data cell occupies more than one cell in z direction" << std::endl;
	    double frac;
	    // front cell
	    frac = dx*dy*(get_cell_zhigh(iz0)-zlow)/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // middle cells
	    for (int zz=iz0+1; zz<iz1; zz++)
	    {
		frac = dx*dy*cell_dz/cell_vol;
		add(ix0,iy0,zz,frac*spec_val);
	    }
	    // back cell
	    frac = dx*dy*(zhigh-get_cell_zlow(iz1))/cell_vol;
	    add(ix0,iy0,iz1,frac*spec_val);
	}

	// data cell occupies more than one cell in x and y direction
	if((ix0!=ix1) && (iy0!=iy1) && (iz0==iz1))
	{
	    if (Debug) std::cout << "data cell occupies more than one cell in x and y direction" << std::endl;
	    double frac;
	    double dxl = get_cell_xhigh(ix0)-xlow;
	    double dyl = get_cell_yhigh(iy0)-ylow;
	    double dxh = xhigh-get_cell_xlow(ix1);
	    double dyh = yhigh-get_cell_ylow(iy1);
	    // left lower cell
	    frac = dxl*dyl*dz/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // right lower cell
	    frac = dxh*dyl*dz/cell_vol;
	    add(ix1,iy0,iz0,frac*spec_val);
	    // left upper cell
	    frac = dxl*dyh*dz/cell_vol;
	    add(ix0,iy1,iz0,frac*spec_val);
	    // right upper cell
	    frac = dxh*dyh*dz/cell_vol;
	    add(ix1,iy1,iz0,frac*spec_val);
	    // intermediate border cells along x
	    for(int xx=ix0+1; xx<ix1; xx++)
	    {
		// lower
		frac = cell_dx*dyl*dz/cell_vol;
		add(xx,iy0,iz0,frac*spec_val);
		// upper
		frac = cell_dx*dyh*dz/cell_vol;
		add(xx,iy1,iz0,frac*spec_val);
	    }
	    // intermediate border cells along y
	    for(int yy=iy0+1; yy<iy1; yy++)
	    {
		// left
		frac = dxl*cell_dy*dz/cell_vol;
		add(ix0,yy,iz0,frac*spec_val);
		// right
		frac = dxh*cell_dy*dz/cell_vol;
		add(ix1,yy,iz0,frac*spec_val);
	    }
	    // middle cells
	    for(int yy=iy0+1; yy<iy1; yy++)
	      for(int xx=ix0+1; xx<ix1; xx++)
	      {
		  frac = cell_dx*cell_dy*dz/cell_vol;
		  add(xx,yy,iz0,frac*spec_val);
	      }
	}

	// data cell occupies more than one cell in x and z direction
	if((ix0!=ix1) && (iy0==iy1) && (iz0!=iz1))
	{
	    if (Debug) std::cout << "data cell occupies more than one cell in x and z direction" << std::endl;
	    double frac;
	    double dxl = get_cell_xhigh(ix0)-xlow;
	    double dzl = get_cell_zhigh(iz0)-zlow;
	    double dxh = xhigh-get_cell_xlow(ix1);
	    double dzh = zhigh-get_cell_zlow(iz1);
	    // left front cell
	    frac = dxl*dy*dzl/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // right front cell
	    frac = dxh*dy*dzl/cell_vol;
	    add(ix1,iy0,iz0,frac*spec_val);
	    // left back cell
	    frac = dxl*dy*dzh/cell_vol;
	    add(ix0,iy0,iz1,frac*spec_val);
	    // right back cell
	    frac = dxh*dy*dzh/cell_vol;
	    add(ix1,iy0,iz1,frac*spec_val);
	    // intermediate border cells along x
	    for(int xx=ix0+1; xx<ix1; xx++)
	    {
		// front
		frac = cell_dx*dy*dzl/cell_vol;
		add(xx,iy0,iz0,frac*spec_val);
		// back
		frac = cell_dx*dy*dzh/cell_vol;
		add(xx,iy0,iz1,frac*spec_val);
	    }
	    // intermediate border cells along z
	    for(int zz=iz0+1; zz<iz1; zz++)
	    {
		// left
		frac = dxl*dy*cell_dz/cell_vol;
		add(ix0,iy0,zz,frac*spec_val);
		// right
		frac = dxh*dy*cell_dz/cell_vol;
		add(ix1,iy0,zz,frac*spec_val);
	    }
	    // middle cells
	    for(int zz=iz0+1; zz<iz1; zz++)
	      for(int xx=ix0+1; xx<ix1; xx++)
	      {
		  frac = cell_dx*dy*cell_dz/cell_vol;
		  add(xx,iy0,zz,frac*spec_val);
	      }
	}

	// data cell occupies more than one cell in y and z direction
	if((ix0==ix1) && (iy0!=iy1) && (iz0!=iz1))
	{
	    if (Debug) std::cout << "data cell occupies more than one cell in y and z direction" << std::endl;
	    double frac;
	    double dyl = get_cell_yhigh(iy0)-ylow;
	    double dzl = get_cell_zhigh(iz0)-zlow;
	    double dyh = yhigh-get_cell_ylow(iy1);
	    double dzh = zhigh-get_cell_zlow(iz1);
	    // lower front cell
	    frac = dx*dyl*dzl/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // upper front cell
	    frac = dx*dyh*dzl/cell_vol;
	    add(ix0,iy1,iz0,frac*spec_val);
	    // lower back cell
	    frac = dx*dyl*dzh/cell_vol;
	    add(ix0,iy0,iz1,frac*spec_val);
	    // upper back cell
	    frac = dx*dyh*dzh/cell_vol;
	    add(ix0,iy1,iz1,frac*spec_val);
	    // intermediate border cells along y
	    for(int yy=iy0+1; yy<iy1; yy++)
	    {
		// front
		frac = dx*cell_dy*dzl/cell_vol;
		add(ix0,yy,iz0,frac*spec_val);
		// back
		frac = dx*cell_dy*dzh/cell_vol;
		add(ix0,yy,iz1,frac*spec_val);
	    }
	    // intermediate border cells along z
	    for(int zz=iz0+1; zz<iz1; zz++)
	    {
		// lower
		frac = dx*dyl*cell_dz/cell_vol;
		add(ix0,iy0,zz,frac*spec_val);
		// upper
		frac = dx*dyh*cell_dz/cell_vol;
		add(ix0,iy1,zz,frac*spec_val);
	    }
	    // middle cells
	    for(int zz=iz0+1; zz<iz1; zz++)
	      for(int yy=iy0+1; yy<iy1; yy++)
	      {
		  frac = dx*cell_dy*cell_dz/cell_vol;
		  add(ix0,yy,zz,frac*spec_val);
	      }
	}

	// data cell occupies a larger volume
	if((ix0!=ix1) && (iy0!=iy1) && (iz0!=iz1))
	{
	    if (Debug) std::cout << "data cell occupies a larger volume, ix0,ix1,iy0,iy1,iz0,iz1: " <<ix0<<" "<<ix1<<" "<<iy0<<" "<<iy1<<" "<<iz0<<" "<<iz1<<std::endl;
	    // do all eight corner cells separately
	    double frac;
	    double dxl = get_cell_xhigh(ix0)-xlow;
	    double dyl = get_cell_yhigh(iy0)-ylow;
	    double dzl = get_cell_zhigh(iz0)-zlow;
	    double dxh = xhigh-get_cell_xlow(ix1);
	    double dyh = yhigh-get_cell_ylow(iy1);
	    double dzh = zhigh-get_cell_zlow(iz1);
	    // left lower front
	    frac = dxl*dyl*dzl/cell_vol;
	    add(ix0,iy0,iz0,frac*spec_val);
	    // right lower front
	    frac = dxh*dyl*dzl/cell_vol;
	    add(ix1,iy0,iz0,frac*spec_val);
	    // left upper front
	    frac = dxl*dyh*dzl/cell_vol;
	    add(ix0,iy1,iz0,frac*spec_val);
	    // right upper front
	    frac = dxh*dyh*dzl/cell_vol;
	    add(ix1,iy1,iz0,frac*spec_val);
	    // left lower back
	    frac = dxl*dyl*dzh/cell_vol;
	    add(ix0,iy0,iz1,frac*spec_val);
	    // right lower back
	    frac = dxh*dyl*dzh/cell_vol;
	    add(ix1,iy0,iz1,frac*spec_val);
	    // left upper back
	    frac = dxl*dyh*dzh/cell_vol;
	    add(ix0,iy1,iz1,frac*spec_val);
	    // right upper back
	    frac = dxh*dyh*dzh/cell_vol;
	    add(ix1,iy1,iz1,frac*spec_val);
	    // intermediate border cells along x
	    for(int xx=ix0+1; xx<ix1; xx++)
	    {
		// lower front
		frac = cell_dx*dyl*dzl/cell_vol;
		add(xx,iy0,iz0,frac*spec_val);
		// upper front
		frac = cell_dx*dyh*dzl/cell_vol;
		add(xx,iy1,iz0,frac*spec_val);
		// lower back
		frac = cell_dx*dyl*dzh/cell_vol;
		add(xx,iy0,iz1,frac*spec_val);
		// upper back
		frac = cell_dx*dyh*dzh/cell_vol;
		add(xx,iy1,iz1,frac*spec_val);
	    }
	    // intermediate border cells along y
	    for(int yy=iy0+1; yy<iy1; yy++)
	    {
		// left front
		frac = dxl*cell_dy*dzl/cell_vol;
		add(ix0,yy,iz0,frac*spec_val);
		// right front
		frac = dxh*cell_dy*dzl/cell_vol;
		add(ix1,yy,iz0,frac*spec_val);
		// left back
		frac = dxl*cell_dy*dzh/cell_vol;
		add(ix0,yy,iz1,frac*spec_val);
		// right back
		frac = dxh*cell_dy*dzh/cell_vol;
		add(ix1,yy,iz1,frac*spec_val);
	    }
	    // intermediate border cells along z
	    for(int zz=iz0+1; zz<iz1; zz++)
	    {
		// left lower
		frac = dxl*dyl*cell_dz/cell_vol;
		add(ix0,iy0,zz,frac*spec_val);
		// right lower
		frac = dxh*dyl*cell_dz/cell_vol;
		add(ix1,iy0,zz,frac*spec_val);
		// left upper
		frac = dxl*dyh*cell_dz/cell_vol;
		add(ix0,iy1,zz,frac*spec_val);
		// right upper
		frac = dxh*dyh*cell_dz/cell_vol;
		add(ix1,iy1,zz,frac*spec_val);
	    }
	    // faces xy
	    for(int xx=ix0+1; xx<ix1; xx++)
	      for(int yy=iy0+1; yy<iy1; yy++)
	      {
		  // front
		  frac = cell_dx*cell_dy*dzl/cell_vol;
		  add(xx,yy,iz0,frac*spec_val);
		  // back
		  frac = cell_dx*cell_dy*dzh/cell_vol;
		  add(xx,yy,iz1,frac*spec_val);
	      }
	    // faces xz
	    for(int xx=ix0+1; xx<ix1; xx++)
	      for(int zz=iz0+1; zz<iz1; zz++)
	      {
		  // lower
		  frac = cell_dx*dyl*cell_dz/cell_vol;
		  add(xx,iy0,zz,frac*spec_val);
		  // upper
		  frac = cell_dx*dyh*cell_dz/cell_vol;
		  add(xx,iy1,zz,frac*spec_val);
	      }
	    // faces yz
	    for(int yy=iy0+1; yy<iy1; yy++)
	      for(int zz=iz0+1; zz<iz1; zz++)
	      {
		  // left
		  frac = dxl*cell_dy*cell_dz/cell_vol;
		  add(ix0,yy,zz,frac*spec_val);
		  // right
		  frac = dxh*cell_dy*cell_dz/cell_vol;
		  add(ix1,yy,zz,frac*spec_val);
	      }
	    // central cells with fraction 1 (data cell covers entire grid cell)
	    for(int zz=iz0+1; zz<iz1; zz++)
	      for(int yy=iy0+1; yy<iy1; yy++)
		for(int xx=ix0+1; xx<ix1; xx++)
		{
		    if (Debug) std::cout << "xx,yy,zz: " << xx << " " << yy << " " << zz << std::endl;
		    add(xx,yy,zz,spec_val);
		}
	} // data cell occupies a larger volume
    }
  }

  private: inline void add(int ix, int iy, int iz, double val)
  {
    if ((ix >= 0) && (ix < dimx) && (iy >= 0) && (iy < dimy) && (iz >= 0) && (iz < dimz))
    {
	long n = index(ix,iy,iz);
	field[n] += val;
	is_set[n] = true;
    }
  }

  //////////////////////////////////////////////////////////////
  // MANIPULATE ENTIRE FIELD
  //////////////////////////////////////////////////////////////

  public: inline void clear()
  {
    for (long n=0; n<ntot; n++)
    {
	field[n] = 0.0;
	is_set[n] = false;
    }
  }

  public: inline void init(double val)
  {
    for (long n=0; n<ntot; n++)
    {
	field[n] = val;
	is_set[n] = true;
    }
  }

  //////////////////////////////////////////////////////////////
  // WRITE OUTPUT / DATA / PICTURE
  //////////////////////////////////////////////////////////////
  
  public: void print_field()
  {
    for (int k=0; k<dimz; k++)
    {
      for (int j=0; j<dimy; j++)
      {
	std::cout<<std::endl;
	for (int i=0; i<dimx; i++)
	{
	  long n = index(i,j,k);
	  std::cout<<" (i,j,k)=("<<i<<","<<j<<","<<k<<")f="<<field[n];
	}
      }
    }
    std::cout<<std::endl;
  }
  

};
#endif