compb-dla-model/DLASystem.cpp

//
//  DLASystem.cpp
//

#include "DLASystem.h"

// this function gets called every step,
//   if there is an active particle then it gets moved,
//   if not then add a particle
void DLASystem::update() {
    if (lastParticleIsActive == 1) {
        moveLastParticle();
    } else if (numParticles < endNum) {
        addParticleOnAddCircle();
        setParticleActive();
    } else {
        this->running = false;
    }
}


void DLASystem::clearParticles() {
    // delete particles and the particle list
    for (int i = 0; i < numParticles; i++) {
        delete particleList[i];
    }
    particleList.clear();
    numParticles = 0;
}

// remove any existing particles and setup initial condition
void DLASystem::Reset() {
    // stop running
    this->running = false;

    clearParticles();

    lastParticleIsActive = 0;

    // set the grid to zero
    for (int i = 0; i < gridSize; i++) {
        for (int j = 0; j < gridSize; j++) {
            grid[i][j] = 0;
        }
    }

    // setup initial condition and parameters
    addCircle = 10;
    killCircle = 2.0 * addCircle;
    clusterRadius = 0.0;
    // add a single particle at the origin
    double pos[] = {0.0, 0.0};
    addParticle(pos);
}

// set the value of a grid cell for a particular position
// note the position has the initial particle at (0,0)
// 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;
}

// 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)];
}

// check if the cluster is big enough and we should stop:
// to be safe, we need the killCircle to be at least 2 less than the edge of the grid
int DLASystem::checkStop() {
    if (killCircle + 2 >= gridSize / 2) {
        this->running = false;
        cerr << "STOP" << endl;
        return 1;
    } else return 0;
}

// add a particle to the system at a specific position
void DLASystem::addParticle(double pos[]) {
    // create a new particle
    Particle *p = new Particle(pos);
    // push_back means "add this to the end of the list"
    particleList.push_back(p);
    numParticles++;

    // pos coordinates should be -gridSize/2 < x < gridSize/2
    setGrid(pos, 1);
}

// add a particle to the system at a random position on the addCircle
// 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));
    if (readGrid(pos) == 0)
        addParticle(pos);
    else
        cerr << "FAIL " << pos[0] << " " << pos[1] << endl;
}

// send back the position of a neighbour of a given grid cell
// NOTE: there is no check that the neighbour is inside the grid,
// this has to be done separately...
void DLASystem::setPosNeighbour(double setpos[], double pos[], int val) {
    switch (val) {
        case 0:
            setpos[0] = pos[0] + 1.0;
            setpos[1] = pos[1];
            break;
        case 1:
            setpos[0] = pos[0] - 1.0;
            setpos[1] = pos[1];
            break;
        case 2:
            setpos[0] = pos[0];
            setpos[1] = pos[1] + 1.0;
            break;
        case 3:
            setpos[0] = pos[0];
            setpos[1] = pos[1] - 1.0;
            break;
    }
}

// when we add a particle to the cluster, we should update the cluster radius
// and the sizes of the addCircle and the killCircle
void DLASystem::updateClusterRadius(double pos[]) {

    double rr = distanceFromOrigin(pos);
    if (rr > clusterRadius) {
        clusterRadius = rr;
        // this is how big addCircle is supposed to be:
        //   either 20% more than cluster radius, or at least 5 bigger.
        double check = clusterRadius * addRatio;
        if (check < clusterRadius + 5)
            check = clusterRadius + 5;
        // if it is smaller then update everything...
        if (addCircle < check) {
            addCircle = check;
            killCircle = killRatio * addCircle;
        }
        checkStop();
    }
}

// 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];

    Particle *lastP = particleList[numParticles - 1];

    setPosNeighbour(newpos, lastP->pos, rr);

    if (distanceFromOrigin(newpos) > killCircle) {
        //cerr << "#deleting particle" << endl;
        setGrid(lastP->pos, 0);
        particleList.pop_back();  // remove particle from particleList
        numParticles--;
        setParticleInactive();
    } 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];
        setGrid(lastP->pos, 1);  // set the new grid site to be occupied

        // check if we stick
        if (checkStick()) {
            //cerr << "stick" << endl;
            setParticleInactive();  // make the particle inactive (stuck)
            updateClusterRadius(lastP->pos);  // update the cluster radius, addCircle, etc.
        }
    }
}

// check if the last particle should stick (to a neighbour)
int DLASystem::checkStick() {
    Particle *lastP = particleList[numParticles - 1];
    // loop over neighbours
    for (int i = 0; i < 4; i++) {
        double checkpos[2];
        setPosNeighbour(checkpos, lastP->pos, i);
        // if the neighbour is occupied and the particle will stick probabilistically.
        if (readGrid(checkpos) == 1 && rgen.random01() < stickProbability) {
            return 1;
        }
    }

    return 0;
}


// constructor
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];
    for (int i = 0; i < gridSize; i++) {
        grid[i] = new int[gridSize];
    }

    // reset initial parameters
    Reset();

    addRatio = 1.2;   // how much bigger the addCircle should be, compared to cluster radius
    killRatio = 1.7;   // how much bigger is the killCircle, compared to the addCircle

    csv.open(csvPath);
}

// destructor
DLASystem::~DLASystem() {
    // strictly we should not print inside the destructor but never mind...
    cerr << "deleting system" << endl;
    // delete the particles
    clearParticles();
    // delete the grid
    for (int i = 0; i < gridSize; i++)
        delete[] grid[i];
    delete[] grid;

    if (csv.is_open()) {
        csv.flush();
        csv.close();
    }
}

void DLASystem::exportData() {
    csv << "x,y" << endl;

    for (auto particle: particleList) {
        csv << particle->pos[0] << "," << particle->pos[1] << endl;
    }
}