|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
|
|
#define N0 100000
|
|
#define LAMBDA_A 0.05
|
|
#define LAMBDA_B 0.02
|
|
#define DT 0.5
|
|
#define TMAX 200.0
|
|
|
|
void aggiorna_atomi(int *N, double p);
|
|
void simula_passo(int *NA, int *NB, double pa, double pb);
|
|
|
|
int main() {
|
|
srand48(1209382);
|
|
|
|
int NA = N0;
|
|
int NB = 0;
|
|
|
|
double pa = LAMBDA_A * DT;
|
|
double pb = LAMBDA_B * DT;
|
|
|
|
int max_NB = 0;
|
|
double max_t;
|
|
|
|
FILE *out = fopen("catena.dat", "w");
|
|
if(out == NULL) {
|
|
printf("Errore nell'apertura del file 'catena.dat'\n");
|
|
exit(1);
|
|
}
|
|
for(double t = 0.; t < TMAX; t += DT) {
|
|
simula_passo(&NA, &NB, pa, pb);
|
|
double NA_th = N0 * exp(-LAMBDA_A * t);
|
|
double NB_th = N0 * LAMBDA_A / (LAMBDA_B - LAMBDA_A) * (exp(-LAMBDA_A * t) - exp(-LAMBDA_B * t));
|
|
|
|
fprintf(out, "%.2f %d %.8f %d %.8f\n", t, NA, NA_th, NB, NB_th);
|
|
|
|
if(NB > max_NB) {
|
|
max_NB = NB;
|
|
max_t = t;
|
|
}
|
|
}
|
|
fclose(out);
|
|
|
|
/* ho calcolato il massimo teorico facendo la derivata dell'eq (2), ma va bene
|
|
anche calcolare il massimo numericamente, come si fa per i valori di simulazione */
|
|
double max_t_th = log(LAMBDA_A / LAMBDA_B) / (LAMBDA_A - LAMBDA_B);
|
|
double max_NB_th = N0 * LAMBDA_A / (LAMBDA_B - LAMBDA_A) * (exp(-LAMBDA_A * max_t_th) - exp(-LAMBDA_B * max_t_th));
|
|
printf("I valori massimi numerici e teorici sono: %d (al tempo %lf) e %lf (al tempo %lf)\n", max_NB, max_t, max_NB_th, max_t_th);
|
|
|
|
return 0;
|
|
}
|
|
|
|
void aggiorna_atomi(int *N, double p) {
|
|
if(drand48() < p) *N += 1;
|
|
}
|
|
|
|
void simula_passo(int *NA, int *NB, double pa, double pb){
|
|
int i;
|
|
int diventa_B = 0;
|
|
int diventa_C = 0;
|
|
for(i = 0; i < *NA; i++) {
|
|
aggiorna_atomi(&diventa_B, pa);
|
|
}
|
|
for(i = 0; i < *NB; i++) {
|
|
aggiorna_atomi(&diventa_C, pb);
|
|
}
|
|
|
|
*NA -= diventa_B;
|
|
*NB += diventa_B - diventa_C;
|
|
}
|