I need to calculate the mean of the coverages of two trapped atoms in a square surface. Such averages are taken after certain number of steps, for example, every 10 monte carlo steps I take the average for such coverages, then I print in a file each average and its respective "y" value (while y is less or equal than 1).

**What I have tried:**
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define MAX_X 4
#define MAX_Y 4
#define ITERATIONS (MAX_Y*MAX_X)
FILE *data;
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;
}
}
}
double calculate_mean(double sum, double count){
double average=0.0;
average=sum/count;
return average;
}
int main(){
data=fopen("data.txt","w");
int i = 0, j = 0;
gstate grid[MAX_Y][MAX_X];
int particle1 = 0, particle2 = 0;
int availcells = MAX_X * MAX_Y;
int fullcells = 0;
int rounds = 0;
double N = 1.0*sizeof(grid)/sizeof(grid[0][0]);
float r;
double cover1 = 0.0, cover2 = 0.0, sumacov = 0.0;
double average1 = 0.0, average2 = 0.0;
double sum1 = 0.0, sum2 = 0.0;
float Y = 0.0;
double MCSTEPS = 10;
srand((unsigned)time(0));
grid_init(grid);
for(int iter = 0; iter < ITERATIONS; iter++){
while (Y <= 1)
{
sum1=0.0;
sum2=0.0;
for(int time = 0; time < MCSTEPS; time++ ) {
generate_coords(&j, &i);
switch (grid[j][i])
{
case S_EMPTY:
r = rand()/(float)RAND_MAX;
if(r <= Y){
grid[j][i] = P1_OCCUPIED;
particle1++;
availcells--;
fullcells++;
}
else{
grid[j][i] = P2_OCCUPIED;
particle2++;
availcells--;
fullcells++;
}
break;
case P1_OCCUPIED:
break;
case P2_OCCUPIED:
break;
}
cover1 = (double)particle1/(double)N;
cover2 = (double)particle2/(double)N;
sumacov = cover1 + cover2;
sum1+=cover1;
sum2+=cover2;
}
average1=calculate_mean(sum1,(double)MCSTEPS);
average2=calculate_mean(sum2, (double)MCSTEPS);
Y = Y + 1/((float)(MCSTEPS*ITERATIONS));
fprintf(data,"%f %f %f\n", average1, average2, Y);
printf("%f %f %f\n", average1, average2, Y);
}
}
printf("The process took %d rounds\n\n", rounds);
printf("#particle1 = %d\n\n", particle1);
printf("#particle2 = %d\n\n", particle2);
printf("#availcells = %d\n\n",availcells);
printf("#fullcells = %d\n\n",fullcells);
fclose(data);
return 0;
}

The problem is that my code almost never selects the particle #1 and I don't know if this is related to the way I'm generating the random number "r" to select between particle #1 or #2. IN consequence the average of the coverage for the particle1 gives 0.000000.

Only to verify: Is the way I'm calculating the coverages and its averages appropiate for each particle?

In the mentioned paper the authors simulate the interactions between both particles. The idea is to plot Y versus both coverages and obtain the same results they obtained before. But first I'm trying to "fill" that surface with atoms, before adding interactions between them. However, I have this problem while trying to choose between particles #1 and #2.

I can't think of any other way to incorporate "Y" in the simulation to be able to choose between these two atoms. I don't have a lot of programming experience but this simulation is not supposed to be difficult, and I don't know what to do.