我需要计算两种可能被困在正方形表面上的粒子的覆盖率的平均值，它们是粒子#1和粒子#2。这些平均值是在一定数量的步骤后获得的，例如，每10个蒙特卡罗步骤，并考虑多次迭代。前面提到的覆盖率定义为
（表面捕获的颗粒数）/（总细胞数）。
模拟的步骤如下：
1.从一个粒子在正方形晶格上的随机碰撞开始。
1.我们选择具有给定概率Y的粒子#1和具有概率1 - Y的粒子#2（“Y由粒子#1在气相中的摩尔分数决定”）。
1.如果选择的粒子是1号粒子，我们随机选择一个位置，如果该位置满了，试验结束，否则，我们检查它的4个最近的邻居，如果其中一个邻居中有2号粒子，则生成3号粒子，否则，1号粒子仍被困在表面上。
1.如果选择的粒子是2号粒子，我们随机选择一个位置，如果该位置满了，试验结束，否则，我们检查它的4个最近的邻居，如果1号粒子在它的一个邻居中，则生成3号粒子，否则，2号粒子仍被困在表面上。
这是我的代码：

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

//define the dimensions of the grid
#define MAX_X 5
#define MAX_Y 5
//define the iterations 
#define ITERATIONS (MAX_Y*MAX_X)
#define Y 0.44

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

double Square( double value )
{
    return value * value;
}

double CalcMean( double number[], int steps)
{
    double sum = 0;
    double pMean = 0;

     // calculate mean

    for (int k = 0; k < steps; k++) {
         sum += number[k];
    }
    
    pMean = sum / (double) steps;   // save mean value

    return pMean;
}  
// returns standard deviation

double CalcStdDev( double number[], int steps, double * pMean )
{
    double sum = 0;
    double stddev = 0;
    
    // calculate standard deviation

    sum = 0;   // NOTE: your code is missing this statement.

    for (int k = 0; k < steps; k++) {
        sum += Square( number[k] - *pMean );
    }
    
    stddev = Square( sum / (double) steps );

    return stddev;
}

//Function that locates the four nearest neighbors of the chosen site, considering periodic boundary conditions
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;
    }
}

    //Function that prints individually the "right" neighbor of the chosen site
    gstate rightneighbor(int x, int 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;   
    }


int main(){
    int i = 0, j = 0;
    int rx, ry, lx, ly, ux, uy, dx, dy;
    gstate grid[MAX_Y][MAX_X];
    double particle1 = 0, particle2 = 0, particle3 = 0; // counters for the number of particle1 and particle2
    double availcells = MAX_X * MAX_Y; //first we initialize with all the cells of the matrix available
    double 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, created3 = 0.0, sumacov = 0.0, covera = 0.0;
    double average1 = 0.0, average2 = 0.0, average3 = 0.0, averagea = 0.0; //we define the average of the coverages of particle 1 and particle 2   
    double std1 = 0.0, std2 = 0.0, std3=0.0;
    int MCSTEPS = 100;
    double num_1[MCSTEPS];
    double num_2[MCSTEPS];
    double num_3[MCSTEPS];
    double num_a[MCSTEPS];
    double sum_1=0.0, sum_2=0.0, sum_3=0.0, sum_a=0.0; //sum of the coverages of both particle 1 and 2. useful to calculate the average of the coverages

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

    sum_1=0.0;
    sum_2=0.0; 
    sum_3=0.0;
    sum_a=0.0;                  
    for(int time = 0; time < MCSTEPS; time++ ) {
        for(int iter = 0; iter < ITERATIONS; iter++){
            //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);
                        rightneighbor(i, j, &rx, &ry);

                        if(grid[ry][rx]==P2_OCCUPIED){
                            grid[ry][rx] = S_EMPTY;
                            particle3++;
                            particle2--;
                            fullcells--;
                            availcells++;
                        } 
                      
                        else{ 
                            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);

                        rightneighbor(i, j, &rx, &ry);

                        if(grid[ry][rx]==P1_OCCUPIED){
                            grid[ry][rx] = S_EMPTY;
                            particle3++;
                            particle1--;
                            fullcells--;
                            availcells++;
                        } 
                    
                        else{ 
                            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 = particle1/N;
            cover2 = particle2/N;
            created3 = particle3/N;
            covera = availcells/N;
        }    
        printf("cover1 = %lf    cover2 = %lf    sumacov = %lf     created3 = %lf    covera = %lf\n\n", cover1,cover2,cover1+cover2, created3, covera);
        
        num_1[time] += cover1;
        num_2[time] += cover2;
        num_3[time] += created3;
        num_a[time] += covera;
    } 

    average1 = CalcMean(num_1, MCSTEPS);
    std1 = CalcStdDev( num_1, MCSTEPS, &average1);
    average2 = CalcMean(num_2, MCSTEPS);
    std2 = CalcStdDev( num_2, MCSTEPS, &average2);
    average3 = CalcMean(num_3, MCSTEPS);
    std3 = CalcStdDev( num_3, MCSTEPS, &average3);

    printf("#particle1 = %lf\n\n", particle1);//total of particle1 adsorbed
    printf("#particle2 = %lf\n\n", particle2);//total of particle2 adsorbed
    printf("#particle3 = %lf\n\n", particle3);//total of particle3 adsorbed
    printf("#availcells = %lf\n\n",availcells);
    printf("#fullcells = %lf\n\n",fullcells); 
    printf("average1 = %lf      average2 = %lf       average3 = %lf      std1 = %lf         std2 = %lf        std3 = %lf      Y = %lf\n\n",average1, average2,average3, std1, std2, std3, Y);
    return 0;
}

问题是我得到的平均结果是
average1 = -nan, average2 = -nan, average3 = -nan, std1 = -nan, std2 = -nan, std3 = -nan
但我不明白为什么我会得到这样的结果“南”是什么意思？
先谢谢你。