I'm simulating the interactions between two kind of particles in a square surface in C.
The steps of my simulations are the following:
1) We start with an empty surface and select at random a site of it.
2) We choose randomly between particle #1 and particle #2. To do so, we define a random number "r" so that if it's less than a given value "Y", particle #1 is chosen. Otherwise, we choose particle #2.
3) If particle #1 is chosen, we look up among the four nearest neighbors of the site. If among the nearest neighbors of such a site, any of them is occupied by particle #2, they react to form particle #3. Otherwise, particle #1 remains adsorbed.
4) If particle #2 is chosen, we look up among the four nearest neighbors of the site. If among the nearest neighbors of such a site, any of them is occupied by particle #1, they react to form particle #3. Otherwise, particle #2 remains adsorbed.
This is the code I have written so far:
What I have tried:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define MAX_X 3
#define MAX_Y 3
#define Y 0.55
typedef enum { S_EMPTY, P1_OCCUPIED, P2_OCCUPIED, S_NONE } gstate;
int gridrnd(int max)
{
return (rand() % max);
}
int generate_coords(int* j, int* i )
{
if (!i || !j)
return 1;
*i = gridrnd(MAX_X);
*j = gridrnd(MAX_Y);
return 0;
}
void grid_init(gstate grid[MAX_Y][MAX_X])
{
for (int j = 0; j < MAX_Y; j++) {
for (int i = 0; i < MAX_X; i++) {
grid[j][i] = S_EMPTY;
}
}
}
int get_limited_coord(int coord, int coord_max){
if (coord >= 0 && coord < coord_max) {
return coord;
} else if (coord >= coord_max) {
return coord - coord_max;
} else {
return coord + coord_max;
}
}
gstate rightneighbor(int *x, int *y, gstate grid[MAX_X][MAX_Y], int *rx, int *ry){
*ry = get_limited_coord(*y, MAX_Y);
*rx = get_limited_coord(*x+1, MAX_X);
printf("right neighbor = (%d,%d)\n\n",*ry,*rx);
return 0;
}
gstate leftneighbor(int *x, int *y, gstate grid[MAX_X][MAX_Y], int *lx, int *ly){
*ly = get_limited_coord(*y , MAX_Y);
*lx = get_limited_coord(*x-1, MAX_X);
printf("left neighbor = (%d,%d)\n\n",*ly,*lx);
return 0;
}
gstate upneighbor(int *x, int *y, gstate grid[MAX_X][MAX_Y], int *ux, int *uy){;
*uy = get_limited_coord(*y-1, MAX_Y);
*ux = get_limited_coord(*x , MAX_X);
printf("up neighbor = (%d,%d)\n\n",*uy,*ux);
return 0;
}
gstate downneighbor(int *x, int *y, gstate grid[MAX_X][MAX_Y], int *dx, int *dy){
*dy = get_limited_coord(*y+1, MAX_Y);
*dx = get_limited_coord(*x , MAX_X);
printf("down neighbor = (%d,%d)\n\n",*dy,*dx);
return 0;
}
typedef struct SimulInfo
{
gstate grid[MAX_Y][MAX_X];
int particle1;
int particle2;
int particle3;
int availcells;
int fullcells;
} SimulInfo;
int adsorb_particle1(SimulInfo * psi,int x, int y){
psi->grid[y][x] = P1_OCCUPIED;
psi->particle1++;
psi->availcells--;
psi->fullcells++;
return 0;
}
int adsorb_particle2(SimulInfo * psi,int x, int y){
psi->grid[y][x] = P2_OCCUPIED;
psi->particle2++;
psi->availcells--;
psi->fullcells++;
return 0;
}
int form_particle3(SimulInfo * psi, int x, int y){
psi->grid[y][x] = S_EMPTY;
psi->particle3++;
psi->availcells++;
psi->fullcells--;
return 0;
}
int main(){
int i = 0, j = 0;
gstate grid[MAX_Y][MAX_X];
int particle[3]={0};
int particle1 = 0, particle2 = 0, particle3 = 0;
int availcells = MAX_X * MAX_Y;
int fullcells = 0;
int rounds = 0;
int rx, ry, lx, ly, ux, uy, dx, dy;
double N = 1.0*sizeof(grid)/sizeof(grid[0][0]);
double r=0, r1=0, r2=0;
SimulInfo psi;
srand((unsigned)time(0));
grid_init(grid);
while(rounds<10){
generate_coords(&j, &i);
switch (grid[j][i])
{
case S_EMPTY:
r = rand()/(float)RAND_MAX;
printf("r = %lf\n",r);
if(r<=Y){
printf("r = %lf is less than Y = %lf. We choose the particle #1\n\n", r, Y);
r1 = rand()/(float)RAND_MAX;
printf("r1 = %lf\n",r1);
rightneighbor(&i, &j, grid, &rx, &ry);
leftneighbor(&i, &j, grid, &lx, &ly);
upneighbor(&i, &j, grid, &ux, &uy);
downneighbor(&i, &j, grid, &dx, &dy);
if(r1<=0.25){
if(psi.grid[ry][rx]==P2_OCCUPIED){
form_particle3(&psi, rx, ry);
}
else{
adsorb_particle1(&psi, rx, ry);
}
}
else if(r1<=0.50){
if(psi.grid[ly][lx]==P2_OCCUPIED){
form_particle3(&psi, lx, ly);
}
else{
adsorb_particle1(&psi, lx, ly);
}
}
else if(r1<=0.75){
if(psi.grid[uy][ux]==P2_OCCUPIED){
form_particle3(&psi, ux, uy);
}
else{
adsorb_particle1(&psi, ux, uy);
}
}
else if(r1<=1.00){
if(psi.grid[dy][dx]==P2_OCCUPIED){
form_particle3(&psi, dx, dy);
}
else{
adsorb_particle1(&psi, dx, dy);
}
}
}
else{
printf("r = %lf is greater than Y = %lf. We choose the particle #2\n\n", r, Y);
r2 = rand()/(float)RAND_MAX;
printf("r2 = %lf\n",r2);
rightneighbor(&i, &j, grid, &rx, &ry);
leftneighbor(&i, &j, grid, &lx, &ly);
upneighbor(&i, &j, grid, &ux, &uy);
downneighbor(&i, &j, grid, &dx, &dy);
if(r2<=0.25){
if(psi.grid[ry][rx]==P1_OCCUPIED){
form_particle3(&psi, rx, ry);
}
else{
adsorb_particle2(&psi, rx, ry);
}
}
else if(r2<=0.50){
if(psi.grid[ly][lx]==P1_OCCUPIED){
form_particle3(&psi, lx, ly);
}
else{
adsorb_particle2(&psi, lx, ly);
}
}
else if(r2<=0.75){
if(psi.grid[uy][ux]==P1_OCCUPIED){
form_particle3(&psi, ux, uy);
}
else{
adsorb_particle2(&psi, ux, uy);
}
}
else if(r2<=1.00){
if(psi.grid[dy][dx]==P1_OCCUPIED){
form_particle3(&psi, dx, dy);
}
else{
adsorb_particle2(&psi, dx, dy);
}
}
}
break;
case P1_OCCUPIED:
break;
case P2_OCCUPIED:
break;
}
rounds++;
}
printf("The process took %d rounds\n\n", rounds);
printf("#particle1 = %d\n\n", psi.particle1);
printf("#particle2 = %d\n\n", psi.particle2);
printf("#particle3 = %d\n\n", psi.particle3);
printf("#availcells = %d\n\n", psi.availcells);
printf("#fullcells = %d\n\n", psi.fullcells);
return 0;
}
As you can see, it is large and I'm not sure how to make it easier to maintain, debug, etc, mainly because I need to add additional particles and interactions in the future. Could you give me some suggestions to rewrite my code in a more condensed way?