branches/exact/examples/degenerate/degenerate2.cc

   1 // Degenerate vertex 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 // The total number of points to create as degenerate vertices
  13 const int points=100;
  14
  15 // The number of planes that will be cut around each point
  16 const int n=64;
  17 const double step=2*pi/n;
  18
  19 // The angle (in radians) of the cutting planes from horizontal
  20 const double theta=0.04;
  21
  22 // This function returns a random floating point number between 0 and 1
  23 double rnd() {return double(rand())/RAND_MAX;}
  24
  25 int main() {
  26         double x,y,z,rsq,r,phi;
  27         voronoicell v;
  28         int n=0;
  29
  30         // Initialize the Voronoi cell to be a cube of side length 2, centered on
  31         // the origin
  32         v.init(-1,1,-1,1,-1,1);
  33
  34         // Plane cutting
  35         while(n<points) {
  36
  37                 // Choose a random point
  38                 x=2*rnd()-1;
  39                 y=2*rnd()-1;
  40                 z=2*rnd()-1;
  41
  42                 // Skip it if it's outside the unit sphere or too close to the
  43                 // origin
  44                 rsq=x*x+y*y+z*z;
  45                 if(rsq<0.01||rsq>1) continue;
  46
  47                 // Rescale the point so that it has modulus 1, and then apply
  48                 // plane cuts around this point
  49                 r=1/sqrt(rsq);x*=r;y*=r;z*=r;
  50                 rsq=sqrt(x*x+y*y);r=z/rsq;
  51                 for(phi=rnd()*step;phi<2*pi;phi+=step)
  52                         v.plane(x*cos(theta)+sin(theta)*(-y*cos(phi)/rsq-x*r*sin(phi)),
  53                                 y*cos(theta)+sin(theta)*(x*cos(phi)/rsq-y*r*sin(phi)),
  54                                 z*cos(theta)+sin(theta)*rsq*sin(phi),1);
  55                 n++;
  56         }
  57
  58         // Output the Voronoi cell to a file in Gnuplot format
  59         v.draw_gnuplot(0,0,0,"degenerate2.gnu");
  60
  61         // Optional POV-Ray output
  62         v.draw_pov(0,0,0,"degenerate2_v.pov");
  63         v.draw_pov_mesh(0,0,0,"degenerate2_m.pov");
  64 }