branches/dynamic/examples/extra/cut_region.cc

   1 // Single Voronoi cell example code
   2 //
   3 // Author   : Chris H. Rycroft (LBL / UC Berkeley)
   4 // Email    : chr@alum.mit.edu
   5 // Date     : July 1st 2008
   6
   7 #include "cell.cc"
   8
   9 const double pi=3.1415926535897932384626433832795;
  10
  11 // This constant sets the tolerance in the bisection search algorithm
  12 const double tolwidth=1e-7;
  13
  14 // This constant determines the density of points to test
  15 const double theta_step=pi/200;
  16
  17 int main() {
  18         double x,y,z,r,rmin,rmax;
  19         double theta,phi,phi_step;
  20         voronoicell v;
  21         ofstream os;
  22         os.open("cell_cut_region.gnu",ofstream::out|ofstream::trunc);
  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("cell.gnu",0,0,0);
  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                         os << x << " " << y << " " << z << endl;
  68                 }
  69         }
  70
  71         os.close();
  72 }