branches/exact/examples/extra/cut_region.cc

   1 // Cell cutting region example code
   2 //
   3 // Author   : Chris H. Rycroft (LBL / UC Berkeley)
   4 // Email    : chr@alum.mit.edu
   5 // Date     : August 30th 2011
   6
   7 #include "voro++.hh"
   8 using namespace voro;
   9
  10 const double pi=3.1415926535897932384626433832795;
  11
  12 // This constant sets the tolerance in the bisection search algorithm
  13 const double tolwidth=1e-7;
  14
  15 // This constant determines the density of points to test
  16 const double theta_step=pi/200;
  17
  18 int main() {
  19         double x,y,z,r,rmin,rmax;
  20         double theta,phi,phi_step;
  21         voronoicell v;
  22         FILE *fp=safe_fopen("cell_cut_region.gnu","w");
  23
  24         // Initialize the Voronoi cell to be an octahedron and make a single
  25         // plane cut to add some variation
  26         v.init_octahedron(1);
  27         v.plane(0.4,0.3,1,0.1);
  28
  29         // Output the cell in gnuplot format
  30         v.draw_gnuplot(0,0,0,"cell.gnu");
  31
  32         // Now test over direction vectors from the center of the sphere. For
  33         // each vector, carry out a search to find the maximum distance along
  34         // that vector that a plane will intersect with cell, and save it to
  35         // the output file.
  36         for(theta=theta_step*0.5;theta<pi;theta+=theta_step) {
  37                 phi_step=2*pi/(int(2*pi*sin(theta)/theta_step)+1);
  38                 for(phi=phi_step*0.5;phi<2*pi;phi+=phi_step) {
  39
  40                         // Calculate a direction to look along
  41                         x=sin(theta)*cos(phi);
  42                         y=sin(theta)*sin(phi);
  43                         z=cos(theta);
  44
  45                         // Now carry out a bisection search. Here, we initialize
  46                         // a minimum and a maximum guess for the distance
  47                         // along this vector. Keep multiplying rmax by two until
  48                         // the plane no longer makes a cut.
  49                         rmin=0;rmax=1;
  50                         while (v.plane_intersects(x,y,z,rmax)) rmax*=2;
  51
  52                         // We now know that the distance is somewhere between
  53                         // rmin and rmax. Test the point halfway in-between
  54                         // these two. If it intersects, then move rmin to this
  55                         // point; otherwise, move rmax there. At each stage the
  56                         // bounding interval is divided by two. Exit when the
  57                         // width of the interval is smaller than the tolerance.
  58                         while (rmax-rmin>tolwidth) {
  59                                 r=(rmax+rmin)*0.5;
  60                                 if (v.plane_intersects(x,y,z,r)) rmin=r;
  61                                 else rmax=r;
  62                         }
  63
  64                         // Output this point to file
  65                         r=(rmax+rmin)*0.5;
  66                         x*=r;y*=r;z*=r;
  67                         fprintf(fp,"%g %g %g\n",x,y,z);
  68                 }
  69         }
  70
  71         fclose(fp);
  72 }