【解決方法】表面の原子の被覆率の平均を取得する方法 (C)

正方形の表面にトラップされた 2 つの原子の被覆率の平均を計算する必要があります。 このような平均は、特定のステップ数の後に取得されます。たとえば、モンテカルロの 10 ステップごとに、そのようなカバレッジの平均を取得し、各平均とそれぞれの「y」値をファイルに出力します (y は 1 以下です)。 ）。

私が試したこと：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

//define the dimensions of the grid
#define MAX_X 4
#define MAX_Y 4
//define the iterations 
#define ITERATIONS (MAX_Y*MAX_X)

FILE *data;

//define the states of a cell in grid`
typedef enum  { S_EMPTY, P1_OCCUPIED, P2_OCCUPIED, S_NONE } gstate;

// help generate random coordinate of the grid
int gridrnd(int max)
{
    return (rand() % max);
}


// generates random coordinates of the grid
int generate_coords(int* j, int* i ) 
{
    if (!i || !j)
        return 1;

    *i = gridrnd(MAX_X);
    *j = gridrnd(MAX_Y);

   // printf("(%d,%d)\n\n", *j, *i);
    return 0;
}

//function to initialize the grid as empty
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;
        }
    }
}

//funstion to calculate the mean of the coverages
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; // counters for the number of particle1 and particle2
    int availcells = MAX_X * MAX_Y; //first we initialize with all the cells of the matrix available
    int fullcells = 0;
    int rounds = 0;
    double N = 1.0*sizeof(grid)/sizeof(grid[0][0]); //number of the total sites in the matrix
    float r;
    double cover1 = 0.0, cover2 = 0.0, sumacov = 0.0;
    double average1 = 0.0, average2 = 0.0; //we define the average of the coverages of particle 1 and particle 2 
    double sum1 = 0.0, sum2 = 0.0; //sum of the coverages of both particle 1 and 2. useful to calculate the average of the coverages
    float Y = 0.0;
    double MCSTEPS = 10;

    srand((unsigned)time(0));
    
    // Initialize grid to be S_EMPTY
    grid_init(grid);

for(int iter = 0; iter < ITERATIONS; iter++){

        while (Y <= 1) // repeat till the grid is full
        {
            sum1=0.0;
            sum2=0.0;
                       
            for(int time = 0; time < MCSTEPS; time++ ) {
                    //LOCATE AN ENTRY OF THE MATRIX RANDOMLY
                    generate_coords(&j, &i);

                    //EVALUATE THE CHOOSEN SITE
                    switch (grid[j][i])
                    {
                        case S_EMPTY:
                            //printf("IT'S S_EMPTY, LET'S FILL IT WITH A PARTICLE. FIRST LET'S GENERATE TO DECIDE IFIT WILL BE TRAPPED\n\n");
                            
                            r = rand()/(float)RAND_MAX;
                            
                            if(r <= Y){//The particle #1 is chosen         
                              //printf("r = %lf is less than Y = %lf. We choose the particle #1\n\n", r, Y);
                                grid[j][i] = P1_OCCUPIED;
                                particle1++;
                                availcells--;
                                fullcells++;
                            }
                            else{//The particle #2 is chosen
                              //printf("r = %lf is greater than Y = %lf. We choose the particle #2\n\n", r, Y);
                                grid[j][i] = P2_OCCUPIED;
                                particle2++;
                                availcells--;
                                fullcells++;                
                            }  
                            break;
                            
                        case P1_OCCUPIED:
                            //printf("IT'S OCCUPIED WITH THE PARTICLE #1. PLEASE, GENERATE ANOTHER SITE ON THE SURFACE\n\n");
                            break;

                        case P2_OCCUPIED:
                            //printf("IT'S OCCUPIED WITH THE PARTICLE #2. PLEASE, GENERATE ANOTHER SITE ON THE SURFACE\n\n");
                            break;
                        
                        }  

                    cover1 = (double)particle1/(double)N;
                    cover2 = (double)particle2/(double)N;
                    sumacov = cover1 + cover2;// To verify if the sum of the coverages is effectively equal to 1

                    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);//total of particle1 adsorbed
    printf("#particle2 = %d\n\n", particle2);//total of particle2 adsorbed
    printf("#availcells = %d\n\n",availcells);
    printf("#fullcells = %d\n\n",fullcells); 

    fclose(data);

    return 0;
}

問題は、私のコードがパーティクル #1 を選択することはほとんどないことです。これが、パーティクル #1 または #2 を選択するために乱数「r」を生成する方法に関連しているかどうかはわかりません。 その結果、粒子 1 のカバレッジの平均は 0.000000 になります。

検証のみ：カバレッジとその平均を計算する方法は、各粒子に適していますか？