#include "voro++.hh"
using namespace voro;
const double h=0.05;
const double hcube=h*h*h;
const int particles=20;
double rnd() {return double(rand())/RAND_MAX;}
int main() {
int i;
double x,y,z,r,rx,ry,rz;
container con(0,1,0,1,0,1,5,5,5,false,false,false,8);
for(i=0;i<particles;i++) {
x=rnd();
y=rnd();
z=rnd();
con.put(i,x,y,z);
}
con.draw_particles("find_voro_cell_p.gnu");
FILE *f1=safe_fopen("find_voro_cell.vec","w");
for(x=0.5*h;x<1;x+=h) for(y=0.5*h;y<1;y+=h) {
if(con.find_voronoi_cell(x,y,0.5,rx,ry,rz,i))
fprintf(f1,"%g %g %g %g %g %g %g\n",x,y,0.5,rx-x,ry-y,rz-0.5,
sqrt((rx-x)*(rx-x)+(ry-y)*(ry-y)+(rz-0.5)*(rz-0.5)));
else fprintf(stderr,"# find_voronoi_cell error for %g %g 0.5\n",x,y);
}
fclose(f1);
int samp_v[particles];
for(i=0;i<particles;i++) samp_v[i]=0;
for(z=0.5*h;z<1;z+=h) for(y=0.5*h;y<1;y+=h) for(x=0.5*h;x<1;x+=h) {
if(con.find_voronoi_cell(x,y,z,rx,ry,rz,i)) samp_v[i]++;
else fprintf(stderr,"# find_voronoi_cell error for %g %g %g\n",x,y,z);
}
f1=safe_fopen("find_voro_cell.vol","w");
FILE *f2=safe_fopen("find_voro_cell_v.gnu","w");
c_loop_all cla(con);
voronoicell c;
if(cla.start()) do if (con.compute_cell(c,cla)) {
cla.pos(i,x,y,z,r);
fprintf(f1,"%d %g %g %g %g %g\n",i,x,y,z,c.volume(),samp_v[i]*hcube);
c.draw_gnuplot(x,y,z,f2);
} while (cla.inc());
fclose(f1);
fclose(f2);
}