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:
-
Start with the random collision of a particle on a square lattice.
-
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").
-
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.
-
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.