vicsek.cpp

#include <iostream>
#include <cmath>
#include <random>
#include <vector>
#include <tuple>
#include <valarray>

#define PI 3.14159265359
using cell_type = std::vector<std::vector<std::vector<int>>>;
using cell_idx = std::tuple<int, int>;

// Viscsek model struct
struct Vicsek{
    int L;
    double rho;
    double v0;
    double R;
    double eta;
    int n_particle;
    std::valarray<double> xs;
    std::valarray<double> ys;
    std::valarray<double> phases;
};

void mod(std::valarray<double>& xs, double L)
{
    for (auto& x: xs)
    {
        x = std::fmod(x+L, L);
    }
}

double orderparameter(Vicsek& m)
{
    auto rx = std::cos(m.phases).sum() / m.n_particle;
    auto ry = std::sin(m.phases).sum() / m.n_particle;
    return std::sqrt(rx*rx + ry*ry);
}

cell_type CellIdx(std::valarray<double> xs, std::valarray<double> ys, int L)
{
    cell_type cellidx(L, std::vector<std::vector<int>>(L, std::vector<int>()));
    for (int i=0; i<xs.size(); i++)
    {
        int x = (int)std::floor(xs[i]);
        int y = (int)std::floor(ys[i]);
        cellidx[y][x].push_back(i);
    }
    return cellidx;
}

std::vector<cell_idx> nearest_cell(double x, double y, int L)
{
    int center_x = (int)std::floor(x);
    int center_y = (int)std::floor(y);
    std::vector<cell_idx> cells;
    for (int k=L-1; k<L+2; k++)
    {
        int idx_x = (center_x + k) % L;
        for (int l=L-1; l<L+2; l++)
        {
            int idx_y = (center_y + l) % L;
            cell_idx cell = std::make_tuple(idx_x, idx_y);
            cells.push_back(cell);
        }
    }
    return cells;
}

// initialize Vicsek model struct
Vicsek init_vicsek(int L, double rho, double v0, double R, double eta){
    // calculate number of particles
    int n_particle = (int)std::floor(L * L * rho);
    // set random device
    std::random_device seed_gen;
    std::mt19937 engine(seed_gen());
    std::uniform_real_distribution<> dist(0.0, 1.0);
    // initialize position/direction of particles
    std::valarray<double> xs(n_particle);
    std::valarray<double> ys(n_particle);
    std::valarray<double> phases(n_particle);
    for (int i = 0; i < n_particle; i++)
    {
        xs[i] = L * dist(engine);
        ys[i] = L * dist(engine);
        phases[i] = 2.0 * PI * dist(engine);
    }
    Vicsek m = {L, rho, v0, R, eta, n_particle, xs, ys, phases};
    return m;
}

void update_phase(Vicsek& m, std::valarray<double> noises){
    auto coss = std::cos(m.phases);
    auto sins = std::sin(m.phases);
    auto cellidx = CellIdx(m.xs, m.ys, m.L);
    double dir_x, dir_y;
    for (int i = 0; i < m.n_particle; i++)
    {
        dir_x = 0.0;
        dir_y = 0.0;
        double xi = m.xs[i];
        double yi = m.ys[i];
        auto cells = nearest_cell(xi, yi, m.L);
        for (auto& cell: cells)
        {
            auto [idx_x, idx_y] = cell;
            for (auto j: cellidx[idx_y][idx_x])
            {
                double dist = std::pow(std::remainder(xi-m.xs[j], m.L), 2) + std::pow(std::remainder(yi-m.ys[j], m.L), 2);
                if (dist <= m.R*m.R)
                {
                    dir_x += coss[j];
                    dir_y += sins[j];
                }
            }
        }
        m.phases[i] = std::atan2(dir_y, dir_x) + noises[i];
    }
}

void update_position(Vicsek& m){
    m.xs += m.v0 * std::cos(m.phases);
    m.ys += m.v0 * std::sin(m.phases);
    mod(m.xs, m.L);
    mod(m.ys, m.L);
}

std::valarray<double> run(Vicsek& m, int n_step)
{
    // set random device
    std::random_device seed_gen;
    std::mt19937 engine(seed_gen());
    std::uniform_real_distribution<> dist(-0.5*m.eta, 0.5*m.eta);

    std::valarray<double> rs(n_step);
    std::valarray<double> noises(m.n_particle);
    // run!!
    for (int i=0; i<n_step; i++)
    {
        // record orderparameter
        rs[i] = orderparameter(m);
        // generate noise
        for (int j=0; j<m.n_particle; j++)
        {
            noises[j] = dist(engine);
        }
        // update phase and position
        update_phase(m, noises);
        update_position(m);
    }

    return rs;
}

int main(void)
{
    // initialize vicsek model
    int L = 512;
    double rho = 1.0/8.0;
    double v0 = 0.5;
    double R = 1.0;
    double eta = 0.55;
    auto m = init_vicsek(L, rho, v0, R, eta);

    auto rs = run(m, 1000);
    for (auto& r: rs)
    {
        std::cout << r << std::endl;
    }
}