branches/2d/examples/extra/intersect.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++_2d.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 phi_step=pi/400;
  17
  18 int main() {
  19         double x,y,r,rmin,rmax;
  20         double phi;
  21         voronoicell_2d v;
  22         FILE *fp=safe_fopen("intersect_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(-1,1,-1,1);
  27         v.plane(1,1,1);
  28
  29         // Output the cell in gnuplot format
  30         v.draw_gnuplot(0,0,"intersect_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(phi=phi_step*0.5;phi<2*pi;phi+=phi_step) {
  37
  38                 // Calculate a direction to look along
  39                 x=cos(phi);
  40                 y=sin(phi);
  41
  42                 // Now carry out a bisection search. Here, we initialize a
  43                 // minimum and a maximum guess for the distance along this
  44                 // vector. Keep multiplying rmax by two until the plane no
  45                 // longer makes a cut.
  46                 rmin=0;rmax=1;
  47                 while (v.plane_intersects(x,y,rmax)) rmax*=2;
  48
  49                 // We now know that the distance is somewhere between rmin and
  50                 // rmax. Test the point halfway in-between these two. If it
  51                 // intersects, then move rmin to this point; otherwise, move
  52                 // rmax there. At each stage the bounding interval is divided
  53                 // by two. Exit when the width of the interval is smaller than
  54                 // the tolerance.
  55                 while (rmax-rmin>tolwidth) {
  56                         r=(rmax+rmin)*0.5;
  57                         if (v.plane_intersects(x,y,r)) rmin=r;
  58                         else rmax=r;
  59                 }
  60
  61                 // Output this point to file
  62                 r=(rmax+rmin)*0.5;
  63                 x*=r;y*=r;
  64                 fprintf(fp,"%g %g\n",x,y);
  65         }
  66
  67         fclose(fp);
  68 }