3D conversion

DLASystem.cpp

@@ -40,7 +40,9 @@ void DLASystem::Reset() {
     // set the grid to zero
     for (int i = 0; i < gridSize; i++) {
         for (int j = 0; j < gridSize; j++) {
-            grid[i][j] = 0;
+            for (int k = 0; k < gridSize; ++k) {
+                grid[i][j][k] = 0;
+            }
         }
     }
 
@@ -49,7 +51,7 @@ void DLASystem::Reset() {
     killCircle = 2.0 * addCircle;
     clusterRadius = 0.0;
     // add a single particle at the origin
-    double pos[] = {0.0, 0.0};
+    double pos[] = {0.0, 0.0, 0.0};
     addParticle(pos);
 }
 
@@ -58,13 +60,13 @@ void DLASystem::Reset() {
 // but this corresponds to the middle of the grid array ie grid[ halfGrid ][ halfGrid ]
 void DLASystem::setGrid(double pos[], int val) {
     int halfGrid = gridSize / 2;
-    grid[(int) (pos[0] + halfGrid)][(int) (pos[1] + halfGrid)] = val;
+    grid[(int) (pos[0] + halfGrid)][(int) (pos[1] + halfGrid)][(int) (pos[2] + halfGrid)] = val;
 }
 
 // read the grid cell for a given position
 int DLASystem::readGrid(double pos[]) {
     int halfGrid = gridSize / 2;
-    return grid[(int) (pos[0] + halfGrid)][(int) (pos[1] + halfGrid)];
+    return grid[(int) (pos[0] + halfGrid)][(int) (pos[1] + halfGrid)][(int) (pos[2] + halfGrid)];
 }
 
 // check if the cluster is big enough and we should stop:
@@ -93,10 +95,14 @@ void DLASystem::addParticle(double pos[]) {
 // if we hit an occupied site then we do nothing except print a message
 // (this should never happen)
 void DLASystem::addParticleOnAddCircle() {
-    double pos[2];
-    double theta = rgen.random01() * 2 * M_PI;
-    pos[0] = ceil(addCircle * cos(theta));
-    pos[1] = ceil(addCircle * sin(theta));
+    double pos[3];
+    double theta = rgen.random01() * M_PI;
+    double phi = rgen.random01() * 2 * M_PI;
+
+    pos[0] = ceil(addCircle * sin(theta) * cos(phi));
+    pos[1] = ceil(addCircle * sin(theta) * sin(phi));
+    pos[2] = ceil(addCircle * cos(theta));
+
     if (readGrid(pos) == 0)
         addParticle(pos);
     else
@@ -111,18 +117,32 @@ void DLASystem::setPosNeighbour(double setpos[], double pos[], int val) {
         case 0:
             setpos[0] = pos[0] + 1.0;
             setpos[1] = pos[1];
+            setpos[2] = pos[2];
             break;
         case 1:
             setpos[0] = pos[0] - 1.0;
             setpos[1] = pos[1];
+            setpos[2] = pos[2];
             break;
         case 2:
             setpos[0] = pos[0];
             setpos[1] = pos[1] + 1.0;
+            setpos[2] = pos[2];
             break;
         case 3:
             setpos[0] = pos[0];
             setpos[1] = pos[1] - 1.0;
+            setpos[2] = pos[2];
+            break;
+        case 4:
+            setpos[0] = pos[0];
+            setpos[1] = pos[1];
+            setpos[2] = pos[2] + 1.0;
+            break;
+        case 5:
+            setpos[0] = pos[0];
+            setpos[1] = pos[1];
+            setpos[2] = pos[2] - 1.0;
             break;
     }
 }
@@ -150,8 +170,8 @@ void DLASystem::updateClusterRadius(double pos[]) {
 
 // make a random move of the last particle in the particleList
 void DLASystem::moveLastParticle() {
-    int rr = rgen.randomInt(4);  // pick a random number in the range 0-3, which direction do we hop?
-    double newpos[2];
+    int rr = rgen.randomInt(6);  // pick a random number in the range 0-3, which direction do we hop?
+    double newpos[3];
 
     Particle *lastP = particleList[numParticles - 1];
 
@@ -163,13 +183,12 @@ void DLASystem::moveLastParticle() {
         particleList.pop_back();  // remove particle from particleList
         numParticles--;
         setParticleInactive();
-    }
-        // check if destination is empty
-    else if (readGrid(newpos) == 0) {
+    } else if (readGrid(newpos) == 0) {
         setGrid(lastP->pos, 0);  // set the old grid site to empty
         // update the position
         particleList[numParticles - 1]->pos[0] = newpos[0];
         particleList[numParticles - 1]->pos[1] = newpos[1];
+        particleList[numParticles - 1]->pos[2] = newpos[2];
         setGrid(lastP->pos, 1);  // set the new grid site to be occupied
 
         // check if we stick
@@ -178,42 +197,41 @@ void DLASystem::moveLastParticle() {
             setParticleInactive();  // make the particle inactive (stuck)
             updateClusterRadius(lastP->pos);  // update the cluster radius, addCircle, etc.
         }
-    } else {
-        // if we get to here then we are trying to move to an occupied site
-        // (this should never happen as long as the sticking probability is 1.0)
-        cerr << "reject " << rr << endl;
-        cerr << lastP->pos[0] << " " << lastP->pos[1] << endl;
-        //cerr << newpos[0] << " " << newpos[1] << " " << (int)newpos[0] << endl;
-        //printOccupied();
     }
 }
 
 // check if the last particle should stick (to a neighbour)
 int DLASystem::checkStick() {
     Particle *lastP = particleList[numParticles - 1];
-    int result = 0;
     // loop over neighbours
-    for (int i = 0; i < 4; i++) {
-        double checkpos[2];
+    for (int i = 0; i < 6; i++) {
+        double checkpos[3];
         setPosNeighbour(checkpos, lastP->pos, i);
-        // if the neighbour is occupied...
-        if (readGrid(checkpos) == 1)
-            result = 1;
+        // if the neighbour is occupied and the particle will stick probabilistically.
+        if (readGrid(checkpos) == 1 && rgen.random01() < stickProbability) {
+            return 1;
+        }
     }
-    return result;
+
+    return 0;
 }
 
 
 // constructor
-DLASystem::DLASystem(int maxParticles, string csvPath) {
+DLASystem::DLASystem(const int maxParticles, const string &csvPath, const double stickProbability) {
     cerr << "creating system, gridSize " << gridSize << endl;
     numParticles = 0;
     endNum = maxParticles;
+    this->stickProbability = stickProbability;
 
     // allocate memory for the grid, remember to free the memory in destructor
-    grid = new int *[gridSize];
+    grid = new int **[gridSize];
     for (int i = 0; i < gridSize; i++) {
-        grid[i] = new int[gridSize];
+        grid[i] = new int*[gridSize];
+
+        for (int j = 0; j < gridSize; ++j) {
+            grid[i][j] = new int[gridSize];
+        }
     }
 
     // reset initial parameters
@@ -232,8 +250,12 @@ DLASystem::~DLASystem() {
     // delete the particles
     clearParticles();
     // delete the grid
-    for (int i = 0; i < gridSize; i++)
+    for (int i = 0; i < gridSize; i++) {
+        for (int j = 0; j < gridSize; ++j) {
+            delete[] grid[i][j];
+        }
         delete[] grid[i];
+    }
     delete[] grid;
 
     if (csv.is_open()) {
@@ -243,9 +265,9 @@ DLASystem::~DLASystem() {
 }
 
 void DLASystem::exportData() {
-    csv << "x,y" << endl;
+    csv << "x,y,z" << endl;
 
     for (auto particle: particleList) {
-        csv << particle->pos[0] << "," << particle->pos[1] << endl;
+        csv << particle->pos[0] << "," << particle->pos[1] << "," << particle->pos[2] << endl;
     }
 }

DLASystem.h

@@ -24,6 +24,11 @@ class DLASystem {
     vector<Particle*> particleList;
     int numParticles;
 
+    /*
+     * The probability that a particle will stick at a given site when adjacent.
+     * */
+    double stickProbability;
+
     // delete particles and clear the particle list
     void clearParticles();
 
@@ -35,7 +40,7 @@ class DLASystem {
   
     // size of grid
     static const int gridSize = 1600;
-    int **grid;  // this will be a 2d array that stores whether each site is occupied
+    int ***grid;  // this will be a 2d array that stores whether each site is occupied
   
     // random number generator, class name is rnd, instance is rgen
     rnd rgen;
@@ -66,7 +71,7 @@ class DLASystem {
     int lastParticleIsActive;
   
     // constructor
-    explicit DLASystem(int maxParticles, string csvPath);
+    explicit DLASystem(int maxParticles, const string &csvPath, double stickProbability);
     // destructor
     ~DLASystem();
   
@@ -88,7 +93,7 @@ class DLASystem {
 
     // return the distance of a given point from the origin
     double distanceFromOrigin(double pos[]) {
-      return sqrt( pos[0]*pos[0] + pos[1]*pos[1] );
+      return sqrt( pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2] );
     }
 
     // set whether there is an active particle in the system or not

Particle.h

@@ -2,7 +2,7 @@
 
 class Particle {
   public:
-    static const int dim = 2;  // we are in two dimensions
+    static const int dim = 3;  // we are in two dimensions
     double *pos;  // pointer to an array of size dim, to store the position
   
     // default constructor

mainDLA.cpp

@@ -3,18 +3,19 @@
 #include "DLASystem.h"
 
 int main(int argc, char **argv) {
-    if (argc != 4) {
-        cerr << "Usage " << argv[0] << " <seed> <maxParticles> <csvPath>" << endl;
+    if (argc != 5) {
+        cerr << "Usage " << argv[0] << " <seed> <maxParticles> <stickProbability> <csvPath>" << endl;
         return 1;
     }
 
     int seed = std::stoi(argv[1]);
     int maxParticles = std::stoi(argv[2]);
-    string csvPath = argv[3];
+    double stickProbability = std::stod(argv[3]);
+    string csvPath = argv[4];
     std::cerr << "Setting seed " << seed << endl;
 
     // create the system
-    auto *sys = new DLASystem(maxParticles,  csvPath);
+    auto *sys = new DLASystem(maxParticles, csvPath, stickProbability);
     sys->setSeed(seed);
     sys->running = true;