H! I need to calculate the mean of the coverages of two types of particles that could be trapped on a square surface. They are particle #1 and particle #2. Such averages are taken after certain number of steps, for example, every 10 Monte Carlo steps, and considering a number of iterations.

The coverages mentioned before are defined as (number of particles trapped in the surface)/(number of total cells).

For illustrative purposes, the following are the steps that my simulation follows:

  1. Start with the random collision of a particle on a square lattice.

  2. We chose particle #1 with a given probability Y and particle #2 with probability 1 — Y. ("Y is determined by the mole fraction of the particle #1 in the gas phase").

  3. If the chosen particle is particle #1, we choose a site on the lattice at random. If that site is full the trial ends. Otherwise, particle #1 is trapped.

  4. If the chosen particle is particle #2, we choose a site on the lattice at random. If that site is full the trial ends. Otherwise, particle #2 is trapped.

This is the code I have so far:

#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)
#define Y 0.55
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;
        }
    }
}




int main(){
    data=fopen("data.txt","w");
    int i = 0, j = 0;
    gstate grid[MAX_Y][MAX_X];
    double particle1 = 0, particle2 = 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, sumacov = 0.0;
    double average1 = 0.0, average2 = 0.0; //we define the average of the coverages of particle 1 and particle 2 
    double num_1[MCSTEPS*ITERATIONS];
    double num_2[MCSTEPS*ITERATIONS];
    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
    double MCSTEPS = 5;

    srand((unsigned)time(0));

    // Initialize grid to be S_EMPTY
    grid_init(grid);

sum1=0.0;
sum2=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);
                                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 = particle1/N;
            cover2 = particle2/N;

            num_1[time] += cover1;
            num_2[time] += cover2;

            sum_1 += num_1[time];
            sum_2 += num_2[time];
        }    
    }  

    average1 = sum_1/(MCSTEPS*ITERATIONS);
    average2 = sum_2/(MCSTEPS*ITERATIONS);

    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); 
    return 0;
}

I would like to calculate the mean and standard deviation of such coverages in separate functions since I will to add more interactions between the particles later on and to do not have parts of my code "spread out".my

The problem is that I'm not sure how to write properly those functions due to the content that I already have inside the for loops.

Note that such averages and standard deviations should be calculated at the end of my loops.

usman_37 commented: The formula which is used in this program is mean = average of the numbers<avariance = (summation( ( Xi – average of numbers) * ( Xi – average of numb +0

for me it generates error here
verage2 = sum_2/(MCSTEPS*ITERATIONS);

commented: There are two main ways to calculate standard deviation: population standard deviation and sample standard deviationhtIf you collect data from all mem +0

You have to declare your variables first. MCSTEPS as a constant expression, and before trying to use it. sum_1, sum_2 need declaring also. You can log your output and variables separately, or however you want, and read the file(s) later.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <iostream>
#include <fstream>
using namespace std;

//define the dimensions of the grid
constexpr auto MAX_X = 4;
constexpr auto MAX_Y = 4;
//define the iterations 
#define ITERATIONS (MAX_Y*MAX_X)
#define Y 0.55
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;
        }
    }
}



int main() {
    ofstream log;
    log.open("log.txt", ios::app);
    FILE* fp;
    errno_t err;

    err = fopen_s(&fp, "data.txt", "w");
    if (err == 0)
    {
        printf("The file 'data.txt' was opened\n");
    }
    else
    {
        printf("The file 'data.txt' was not opened\n");
    }
    constexpr double MCSTEPS = 5.0;
    double sum_1 = 0;
    double sum_2 = 0;
    int i = 0, j = 0;
    gstate grid[MAX_Y][MAX_X];
    double particle1 = 0, particle2 = 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, sumacov = 0.0;
    double average1 = 0.0, average2 = 0.0; //we define the average of the coverages of particle 1 and particle 2 
    double num_1[static_cast<int>(MCSTEPS) * ITERATIONS] = {};
    double num_2[static_cast<int>(MCSTEPS) * ITERATIONS] = {};
    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

    srand((unsigned)time(0));

    // Initialize grid to be S_EMPTY
    grid_init(grid);

    sum1 = 0.0;
    sum2 = 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);
                    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 = particle1 / N;
            cover2 = particle2 / N;

            num_1[time] += cover1;
            num_2[time] += cover2;

            sum_1 += num_1[time];
            sum_2 += num_2[time];
        }
    }

    average1 = sum_1 / (MCSTEPS * ITERATIONS);
    average2 = sum_2 / (MCSTEPS * ITERATIONS);
    rounds = ITERATIONS;

    std::cout << "The process took " << rounds << std::endl; // ITERATIONS?
    std::cout << "#particle1 = " << particle1 << std::endl;
    std::cout << "#particle2 = " << particle2 << std::endl;
    std::cout << "#availcells = " << availcells << std::endl;
    std::cout << "#fullcells = " << fullcells << std::endl;

    std::cout << "#cover1 = " << cover1 << std::endl;
    std::cout << "#cover2 = " << cover2 << std::endl;

    std::cout << "#sum1 = " << sum_1 << std::endl;
    std::cout << "#sum2 = " << sum_2 << std::endl; 

    std::cout << "#average1 = " << average1 << std::endl;
    std::cout << "#average2 = " << average2 << std::endl;

    log << "The process took " << rounds << std::endl; // ITERATIONS?
    log << "#particle1 = " << particle1 << std::endl;
    log << "#particle2 = " << particle2 << std::endl;
    log << "#availcells = " << availcells << std::endl;
    log << "#fullcells = " << fullcells << std::endl;

    log << "#cover1 = " << cover1 << std::endl;
    log << "#cover2 = " << cover2 << std::endl;

    log << "#sum1 = " << sum_1 << std::endl;
    log << "#sum2 = " << sum_2 << std::endl;

    log << "#average1 = " << average1 << std::endl;
    log << "#average2 = " << average2 << std::endl;
    log.close();
    return 0;
}



The file 'data.txt' was opened
The process took 16
#particle1 = 10
#particle2 = 6
#availcells = 0
#fullcells = 16
#cover1 = 0.625
#cover2 = 0.375
#sum1 = 334.875
#sum2 = 199.188
#average1 = 4.18594
#average2 = 2.48984
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.