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