/*****************************************************************************

    TRAVIS - Trajectory Analyzer and Visualizer

    http://www.travis-analyzer.de/

    Copyright (c) 2009-2020 Martin Brehm
                  2012-2020 Martin Thomas
                  2016-2020 Sascha Gehrke

    Please cite:  J. Chem. Phys. 2020, 152 (16), 164105.         (DOI 10.1063/5.0005078 )
                  J. Chem. Inf. Model. 2011, 51 (8), 2007-2023.  (DOI 10.1021/ci200217w )

    ---------------------------------------------------------------------------

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*****************************************************************************/

// Voro++, a 3D cell-based Voronoi library
//
// Author   : Chris H. Rycroft (LBL / UC Berkeley)
// Email    : chr@alum.mit.edu
// Date     : August 30th 2011

/** \file container_prd.cc
 * \brief Function implementations for the container_periodic_base and
 * related classes. */


 // This must always be the first include directive
#include "config.h"

#include "v_container_prd.h"
#include "globalvar.h"


const char *GetRevisionInfo_v_container_prd(unsigned int len) {
	static char buf[256];
	GET_REVISION_INFO( buf, len );
	return buf;
}


const char *GetSourceVersion_v_container_prd() {
	static char buf[256];
	GET_SOURCE_VERSION( buf );
	return buf;
}


/** The class constructor sets up the geometry of container, initializing the
 * minimum and maximum coordinates in each direction, and setting whether each
 * direction is periodic or not. It divides the container into a rectangular
 * grid of blocks, and allocates memory for each of these for storing particle
 * positions and IDs.
 * \param[in] (bx_) The x coordinate of the first unit vector.
 * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
 * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
 *                            vector.
 * \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
 *                       coordinate directions.
 * \param[in] init_mem_ the initial memory allocation for each block.
 * \param[in] ps_ the number of floating point entries to store for each
 *                particle. */
container_periodic_base::container_periodic_base(double bx_,double bxy_,double by_,
		double bxz_,double byz_,double bz_,int nx_,int ny_,int nz_,int init_mem_,int ps_)
	: unitcell(bx_,bxy_,by_,bxz_,byz_,bz_), voro_base(nx_,ny_,nz_,bx_/nx_,by_/ny_,bz_/nz_),
	ey(int(max_uv_y*ysp+1)), ez(int(max_uv_z*zsp+1)), wy(ny+ey), wz(nz+ez),
	oy(ny+2*ey), oz(nz+2*ez), oxyz(nx*oy*oz), id(new int*[oxyz]), p(new double*[oxyz]),
	co(new int[oxyz]), mem(new int[oxyz]), img(new char[oxyz]), init_mem(init_mem_), ps(ps_) {
	int i,j,k,l;

	// Clear the global arrays
	int *pp=co;while(pp<co+oxyz) *(pp++)=0;
	pp=mem;while(pp<mem+oxyz) *(pp++)=0;
	char *cp=img;while(cp<img+oxyz) *(cp++)=0;

	// Set up memory for the blocks in the primary domain
	for(k=ez;k<wz;k++) for(j=ey;j<wy;j++) for(i=0;i<nx;i++) {
		l=i+nx*(j+oy*k);
		mem[l]=init_mem;
		id[l]=new int[init_mem];
		p[l]=new double[ps*init_mem];
	}
}

/** The container destructor frees the dynamically allocated memory. */
container_periodic_base::~container_periodic_base() {
	int l;
	for(l=0;l<oxyz;l++) if(mem[l]>0) delete [] p[l];
	for(l=0;l<oxyz;l++) if(mem[l]>0) delete [] id[l];
	delete [] p;
	delete [] id;
	delete [] mem;
	delete [] co;
}


/** The class constructor sets up the geometry of container.
 * \param[in] (bx_) The x coordinate of the first unit vector.
 * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
 * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
 *                            vector.
 * \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
 *			    coordinate directions.
 * \param[in] init_mem_ the initial memory allocation for each block. */
container_periodic_poly::container_periodic_poly(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_,
	int nx_,int ny_,int nz_,int init_mem_)
	: container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,4),
	vc(*this,2*nx_+1,2*ey+1,2*ez+1) {ppr=p;}

/** Put a particle into the correct region of the container.
 * \param[in] n the numerical ID of the inserted particle.
 * \param[in] (x,y,z) the position vector of the inserted particle.
 * \param[in] r the radius of the particle. */
void container_periodic_poly::put(int n,double x,double y,double z,double r) {
	int ijk;
	put_locate_block(ijk,x,y,z);
	id[ijk][co[ijk]]=n;
	double *pp=p[ijk]+4*co[ijk]++;
	*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
	if(max_radius<r) max_radius=r;
}

/** Put a particle into the correct region of the container.
 * \param[in] n the numerical ID of the inserted particle.
 * \param[in] (x,y,z) the position vector of the inserted particle.
 * \param[in] r the radius of the particle.
 * \param[out] (ai,aj,ak) the periodic image displacement that the particle is
 * 			  in, with (0,0,0) corresponding to the primary domain.
 */
void container_periodic_poly::put(int n,double x,double y,double z,double r,int &ai,int &aj,int &ak) {
	int ijk;
	put_locate_block(ijk,x,y,z,ai,aj,ak);
	id[ijk][co[ijk]]=n;
	double *pp=p[ijk]+4*co[ijk]++;
	*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
	if(max_radius<r) max_radius=r;
}


/** Takes a particle position vector and computes the region index into which
 * it should be stored. If the container is periodic, then the routine also
 * maps the particle position to ensure it is in the primary domain. If the
 * container is not periodic, the routine bails out.
 * \param[out] ijk the region index.
 * \param[in,out] (x,y,z) the particle position, remapped into the primary
 *                        domain if necessary.
 * \return True if the particle can be successfully placed into the container,
 * false otherwise. */
void container_periodic_base::put_locate_block(int &ijk,double &x,double &y,double &z) {

	// Remap particle in the z direction if necessary
	int k=step_int(z*zsp);
	if(k<0||k>=nz) {
		int ak=step_div(k,nz);
		z-=ak*bz;y-=ak*byz;x-=ak*bxz;k-=ak*nz;
	}

	// Remap particle in the y direction if necessary
	int j=step_int(y*ysp);
	if(j<0||j>=ny) {
		int aj=step_div(j,ny);
		y-=aj*by;x-=aj*bxy;j-=aj*ny;
	}

	// Remap particle in the x direction if necessary
	ijk=step_int(x*xsp);
	if(ijk<0||ijk>=nx) {
		int ai=step_div(ijk,nx);
		x-=ai*bx;ijk-=ai*nx;
	}

	// Compute the block index and check memory allocation
	j+=ey;k+=ez;
	ijk+=nx*(j+oy*k);
	if(co[ijk]==mem[ijk]) add_particle_memory(ijk);
}

/** Takes a particle position vector and computes the region index into which
 * it should be stored. If the container is periodic, then the routine also
 * maps the particle position to ensure it is in the primary domain. If the
 * container is not periodic, the routine bails out.
 * \param[out] ijk the region index.
 * \param[in,out] (x,y,z) the particle position, remapped into the primary
 *                        domain if necessary.
 * \param[out] (ai,aj,ak) the periodic image displacement that the particle is
 *                        in, with (0,0,0) corresponding to the primary domain.
 * \return True if the particle can be successfully placed into the container,
 * false otherwise. */
void container_periodic_base::put_locate_block(int &ijk,double &x,double &y,double &z,int &ai,int &aj,int &ak) {

	// Remap particle in the z direction if necessary
	int k=step_int(z*zsp);
	if(k<0||k>=nz) {
		ak=step_div(k,nz);
		z-=ak*bz;y-=ak*byz;x-=ak*bxz;k-=ak*nz;
	} else ak=0;

	// Remap particle in the y direction if necessary
	int j=step_int(y*ysp);
	if(j<0||j>=ny) {
		aj=step_div(j,ny);
		y-=aj*by;x-=aj*bxy;j-=aj*ny;
	} else aj=0;

	// Remap particle in the x direction if necessary
	ijk=step_int(x*xsp);
	if(ijk<0||ijk>=nx) {
		ai=step_div(ijk,nx);
		x-=ai*bx;ijk-=ai*nx;
	} else ai=0;

	// Compute the block index and check memory allocation
	j+=ey;k+=ez;
	ijk+=nx*(j+oy*k);
	if(co[ijk]==mem[ijk]) add_particle_memory(ijk);
}

/** Takes a position vector and remaps it into the primary domain.
 * \param[out] (ai,aj,ak) the periodic image displacement that the vector is in,
 *                        with (0,0,0) corresponding to the primary domain.
 * \param[out] (ci,cj,ck) the index of the block that the position vector is
 *                        within, once it has been remapped.
 * \param[in,out] (x,y,z) the position vector to consider, which is remapped
 *                        into the primary domain during the routine.
 * \param[out] ijk the block index that the vector is within. */
inline void container_periodic_base::remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk) {

	// Remap particle in the z direction if necessary
	ck=step_int(z*zsp);
	if(ck<0||ck>=nz) {
		ak=step_div(ck,nz);
		z-=ak*bz;y-=ak*byz;x-=ak*bxz;ck-=ak*nz;
	} else ak=0;

	// Remap particle in the y direction if necessary
	cj=step_int(y*ysp);
	if(cj<0||cj>=ny) {
		aj=step_div(cj,ny);
		y-=aj*by;x-=aj*bxy;cj-=aj*ny;
	} else aj=0;

	// Remap particle in the x direction if necessary
	ci=step_int(x*xsp);
	if(ci<0||ci>=nx) {
		ai=step_div(ci,nx);
		x-=ai*bx;ci-=ai*nx;
	} else ai=0;

	cj+=ey;ck+=ez;
	ijk=ci+nx*(cj+oy*ck);
}


/** Takes a vector and finds the particle whose Voronoi cell contains that
 * vector. Additional wall classes are not considered by this routine.
 * \param[in] (x,y,z) the vector to test.
 * \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell
 *                        contains the vector. If the container is periodic,
 *                        this may point to a particle in a periodic image of
 *                        the primary domain.
 * \param[out] pid the ID of the particle.
 * \return True if a particle was found. If the container has no particles,
 * then the search will not find a Voronoi cell and false is returned. */
bool container_periodic_poly::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) {
	int ai,aj,ak,ci,cj,ck,ijk;
	particle_record w;
	double mrs;

	// Remap the vector into the primary domain and then search for the
	// Voronoi cell that it is within
	remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
	vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);

	if(w.ijk!=-1) {

		// Assemble the position vector of the particle to be returned,
		// applying a periodic remapping if necessary
		ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);
		rx=p[w.ijk][4*w.l]+ak*bxz+aj*bxy+ai*bx;
		ry=p[w.ijk][4*w.l+1]+ak*byz+aj*by;
		rz=p[w.ijk][4*w.l+2]+ak*bz;
		pid=id[w.ijk][w.l];
		return true;
	}
	return false;
}

/** Increase memory for a particular region.
 * \param[in] i the index of the region to reallocate. */
void container_periodic_base::add_particle_memory(int i) {

	// Handle the case when no memory has been allocated for this block
	if(mem[i]==0) {
		mem[i]=init_mem;
		id[i]=new int[init_mem];
		p[i]=new double[ps*init_mem];
		return;
	}

	// Otherwise, double the memory allocation for this block. Carry out a
	// check on the memory allocation size, and print a status message if
	// requested.
	int l,nmem(mem[i]<<1);
	if(nmem>max_particle_memory)
		voro_fatal_error("Absolute maximum memory allocation exceeded",VOROPP_MEMORY_ERROR);
#if VOROPP_VERBOSE >=3
	if ( !g_bVoroSilent )
		mprintf("Particle memory in region %d scaled up to %d\n",i,nmem);
	if (nmem > g_iVoroMemory)
		g_iVoroMemory = nmem;
#endif

	// Allocate new memory and copy in the contents of the old arrays
	int *idp=new int[nmem];
	for(l=0;l<co[i];l++) idp[l]=id[i][l];
	double *pp=new double[ps*nmem];
	for(l=0;l<ps*co[i];l++) pp[l]=p[i][l];

	// Update pointers and delete old arrays
	mem[i]=nmem;
	delete [] id[i];id[i]=idp;
	delete [] p[i];p[i]=pp;
}


/** Import a list of particles from an open file stream into the container.
 * Entries of five numbers (Particle ID, x position, y position, z position,
 * radius) are searched for. If the file cannot be successfully read, then the
 * routine causes a fatal error.
 * \param[in] fp the file handle to read from. */
void container_periodic_poly::import(FILE *fp) {
	int i,j;
	double x,y,z,r;
	while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(i,x,y,z,r);
	if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}

/** Outputs the a list of all the container regions along with the number of
 * particles stored within each. */
void container_periodic_base::region_count() {
	int i,j,k,*cop=co;
	for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)
		mprintf("Region (%d,%d,%d): %d particles\n",i,j,k,*(cop++));
}


/** Clears a container of particles, also clearing resetting the maximum radius
 * to zero. */
void container_periodic_poly::clear() {
	for(int *cop=co;cop<co+nxyz;cop++) *cop=0;
	max_radius=0;
}


/** Computes all the Voronoi cells and saves customized
 * information about them.
 * \param[in] format the custom output string to use.
 * \param[in] fp a file handle to write to. */
void container_periodic_poly::print_custom(const char *format,FILE *fp) {
	c_loop_all_periodic vl(*this);
	print_custom(vl,format,fp);
}


/** Computes all the Voronoi cells and saves customized
 * information about them
 * \param[in] format the custom output string to use.
 * \param[in] filename the name of the file to write to. */
void container_periodic_poly::print_custom(const char *format,const char *filename) {
	FILE *fp=safe_fopen(filename,"w");
	print_custom(format,fp);
	fclose(fp);
}


/** Computes all of the Voronoi cells in the container, but does nothing
 * with the output. It is useful for measuring the pure computation time
 * of the Voronoi algorithm, without any additional calculations such as
 * volume evaluation or cell output. */
void container_periodic_poly::compute_all_cells() {
	voronoicell_neighbor c;
	c_loop_all_periodic vl(*this);
	if(vl.start()) do compute_cell(c,vl);while(vl.inc());
}


/** Calculates all of the Voronoi cells and sums their volumes. In most cases
 * without walls, the sum of the Voronoi cell volumes should equal the volume
 * of the container to numerical precision.
 * \return The sum of all of the computed Voronoi volumes. */
double container_periodic_poly::sum_cell_volumes() {
	voronoicell_neighbor c;
	double vol=0;
	c_loop_all_periodic vl(*this);
	if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc());
	return vol;
}

/** This routine creates all periodic images of the particles. It is meant for
 * diagnostic purposes only, since usually periodic images are dynamically
 * created in when they are referenced. */
void container_periodic_base::create_all_images() {
	int i,j,k;
	for(k=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++) create_periodic_image(i,j,k);
}

/** Checks that the particles within each block lie within that block's bounds.
 * This is useful for diagnosing problems with periodic image computation. */
void container_periodic_base::check_compartmentalized() {
	int c,l,i,j,k;
	double mix,miy,miz,max,may,maz,*pp;
	for(k=l=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++,l++) if(mem[l]>0) {

		// Compute the block's bounds, adding in a small tolerance
		mix=i*boxx-tolerance;max=mix+boxx+tolerance;
		miy=(j-ey)*boxy-tolerance;may=miy+boxy+tolerance;
		miz=(k-ez)*boxz-tolerance;maz=miz+boxz+tolerance;

		// Print entries for any particles that lie outside the block's
		// bounds
		for(pp=p[l],c=0;c<co[l];c++,pp+=ps) if(*pp<mix||*pp>max||pp[1]<miy||pp[1]>may||pp[2]<miz||pp[2]>maz)
			mprintf("%d %d %d %d %f %f %f %f %f %f %f %f %f\n",
			       id[l][c],i,j,k,*pp,pp[1],pp[2],mix,max,miy,may,miz,maz);
	}
}

/** Creates particles within an image block that is aligned with the primary
 * domain in the z axis. In this case, the image block may be comprised of
 * particles from two primary blocks. The routine considers these two primary
 * blocks, and adds the needed particles to the image. The remaining particles
 * from the primary blocks are also filled into the neighboring images.
 * \param[in] (di,dj,dk) the index of the block to consider. The z index must
 *			 satisfy ez<=dk<wz. */
void container_periodic_base::create_side_image(int di,int dj,int dk) {
	int l,dijk=di+nx*(dj+oy*dk),odijk,ima=step_div(dj-ey,ny);
	int qua=di+step_int(-ima*bxy*xsp),quadiv=step_div(qua,nx);
	int fi=qua-quadiv*nx,fijk=fi+nx*(dj-ima*ny+oy*dk);
	double dis=ima*bxy+quadiv*bx,switchx=di*boxx-ima*bxy-quadiv*bx,adis;

	// Left image computation
	if((img[dijk]&1)==0) {
		if(di>0) {
			odijk=dijk-1;adis=dis;
		} else {
			odijk=dijk+nx-1;adis=dis+bx;
		}
		img[odijk]|=2;
		for(l=0;l<co[fijk];l++) {
			if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,dis,by*ima,0);
			else put_image(odijk,fijk,l,adis,by*ima,0);
		}
	}

	// Right image computation
	if((img[dijk]&2)==0) {
		if(fi==nx-1) {
			fijk+=1-nx;switchx+=(1-nx)*boxx;dis+=bx;
		} else {
			fijk++;switchx+=boxx;
		}
		if(di==nx-1) {
			odijk=dijk-nx+1;adis=dis-bx;
		} else {
			odijk=dijk+1;adis=dis;
		}
		img[odijk]|=1;
		for(l=0;l<co[fijk];l++) {
			if(p[fijk][ps*l]<switchx) put_image(dijk,fijk,l,dis,by*ima,0);
			else put_image(odijk,fijk,l,adis,by*ima,0);
		}
	}

	// All contributions to the block now added, so set both two bits of
	// the image information
	img[dijk]=3;
}

/** Creates particles within an image block that is not aligned with the
 * primary domain in the z axis. In this case, the image block may be comprised
 * of particles from four primary blocks. The routine considers these four
 * primary blocks, and adds the needed particles to the image. The remaining
 * particles from the primary blocks are also filled into the neighboring
 * images.
 * \param[in] (di,dj,dk) the index of the block to consider. The z index must
 *			 satisfy dk<ez or dk>=wz. */
void container_periodic_base::create_vertical_image(int di,int dj,int dk) {
	int l,dijk=di+nx*(dj+oy*dk),dijkl,dijkr,ima=step_div(dk-ez,nz);
	int qj=dj+step_int(-ima*byz*ysp),qjdiv=step_div(qj-ey,ny);
	int qi=di+step_int((-ima*bxz-qjdiv*bxy)*xsp),qidiv=step_div(qi,nx);
	int fi=qi-qidiv*nx,fj=qj-qjdiv*ny,fijk=fi+nx*(fj+oy*(dk-ima*nz)),fijk2;
	double disy=ima*byz+qjdiv*by,switchy=(dj-ey)*boxy-ima*byz-qjdiv*by;
	double disx=ima*bxz+qjdiv*bxy+qidiv*bx,switchx=di*boxx-ima*bxz-qjdiv*bxy-qidiv*bx;
	double switchx2,disxl,disxr,disx2,disxr2;

	if(di==0) {dijkl=dijk+nx-1;disxl=disx+bx;}
	else {dijkl=dijk-1;disxl=disx;}

	if(di==nx-1) {dijkr=dijk-nx+1;disxr=disx-bx;}
	else {dijkr=dijk+1;disxr=disx;}

	// Down-left image computation
	bool y_exist=dj!=0;
	if((img[dijk]&1)==0) {
		img[dijkl]|=2;
		if(y_exist) {
			img[dijkl-nx]|=8;
			img[dijk-nx]|=4;
		}
		for(l=0;l<co[fijk];l++) {
			if(p[fijk][ps*l+1]>switchy) {
				if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,disx,disy,bz*ima);
				else put_image(dijkl,fijk,l,disxl,disy,bz*ima);
			} else {
				if(!y_exist) continue;
				if(p[fijk][ps*l]>switchx) put_image(dijk-nx,fijk,l,disx,disy,bz*ima);
				else put_image(dijkl-nx,fijk,l,disxl,disy,bz*ima);
			}
		}
	}

	// Down-right image computation
	if((img[dijk]&2)==0) {
		if(fi==nx-1) {
			fijk2=fijk+1-nx;switchx2=switchx+(1-nx)*boxx;disx2=disx+bx;disxr2=disxr+bx;
		} else {
			fijk2=fijk+1;switchx2=switchx+boxx;disx2=disx;disxr2=disxr;
		}
		img[dijkr]|=1;
		if(y_exist) {
			img[dijkr-nx]|=4;
			img[dijk-nx]|=8;
		}
		for(l=0;l<co[fijk2];l++) {
			if(p[fijk2][ps*l+1]>switchy) {
				if(p[fijk2][ps*l]>switchx2) put_image(dijkr,fijk2,l,disxr2,disy,bz*ima);
				else put_image(dijk,fijk2,l,disx2,disy,bz*ima);
			} else {
				if(!y_exist) continue;
				if(p[fijk2][ps*l]>switchx2) put_image(dijkr-nx,fijk2,l,disxr2,disy,bz*ima);
				else put_image(dijk-nx,fijk2,l,disx2,disy,bz*ima);
			}
		}
	}

	// Recomputation of some intermediate quantities for boundary cases
	if(fj==wy-1) {
		fijk+=nx*(1-ny)-fi;
		switchy+=(1-ny)*boxy;
		disy+=by;
		qi=di+step_int(-(ima*bxz+(qjdiv+1)*bxy)*xsp);
		int dqidiv=step_div(qi,nx)-qidiv;qidiv+=dqidiv;
		fi=qi-qidiv*nx;
		fijk+=fi;
		disx+=bxy+bx*dqidiv;
		disxl+=bxy+bx*dqidiv;
		disxr+=bxy+bx*dqidiv;
		switchx-=bxy+bx*dqidiv;
	} else {
		fijk+=nx;switchy+=boxy;
	}

	// Up-left image computation
	y_exist=dj!=oy-1;
	if((img[dijk]&4)==0) {
		img[dijkl]|=8;
		if(y_exist) {
			img[dijkl+nx]|=2;
			img[dijk+nx]|=1;
		}
		for(l=0;l<co[fijk];l++) {
			if(p[fijk][ps*l+1]>switchy) {
				if(!y_exist) continue;
				if(p[fijk][ps*l]>switchx) put_image(dijk+nx,fijk,l,disx,disy,bz*ima);
				else put_image(dijkl+nx,fijk,l,disxl,disy,bz*ima);
			} else {
				if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,disx,disy,bz*ima);
				else put_image(dijkl,fijk,l,disxl,disy,bz*ima);
			}
		}
	}

	// Up-right image computation
	if((img[dijk]&8)==0) {
		if(fi==nx-1) {
			fijk2=fijk+1-nx;switchx2=switchx+(1-nx)*boxx;disx2=disx+bx;disxr2=disxr+bx;
		} else {
			fijk2=fijk+1;switchx2=switchx+boxx;disx2=disx;disxr2=disxr;
		}
		img[dijkr]|=4;
		if(y_exist) {
			img[dijkr+nx]|=1;
			img[dijk+nx]|=2;
		}
		for(l=0;l<co[fijk2];l++) {
			if(p[fijk2][ps*l+1]>switchy) {
				if(!y_exist) continue;
				if(p[fijk2][ps*l]>switchx2) put_image(dijkr+nx,fijk2,l,disxr2,disy,bz*ima);
				else put_image(dijk+nx,fijk2,l,disx2,disy,bz*ima);
			} else {
				if(p[fijk2][ps*l]>switchx2) put_image(dijkr,fijk2,l,disxr2,disy,bz*ima);
				else put_image(dijk,fijk2,l,disx2,disy,bz*ima);
			}
		}
	}

	// All contributions to the block now added, so set all four bits of
	// the image information
	img[dijk]=15;
}

/** Copies a particle position from the primary domain into an image block.
 * \param[in] reg the block index within the primary domain that the particle
 *                is within.
 * \param[in] fijk the index of the image block.
 * \param[in] l the index of the particle entry within the primary block.
 * \param[in] (dx,dy,dz) the displacement vector to add to the particle. */
void container_periodic_base::put_image(int reg,int fijk,int l,double dx,double dy,double dz) {
	if(co[reg]==mem[reg]) add_particle_memory(reg);
	double *p1=p[reg]+ps*co[reg],*p2=p[fijk]+ps*l;
	*(p1++)=*(p2++)+dx;
	*(p1++)=*(p2++)+dy;
	*p1=*p2+dz;
	if(ps==4) *(++p1)=*(++p2);
	id[reg][co[reg]++]=id[fijk][l];
}