#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
const double pi = 3.1415926; // value of pi
double a = 1; // inner radius of shield
double b = 1.5; // outer radius of shield
double tau = 0.2; // mean free time for v=1
double p_a = 0.05; // absorption probability
double f = 0.95; // inelasticity parameter
double v_0 = 1; // initial speed
double dotProduct(double * v1, double * v2) {
return v1[0]*v2[0] + v1[1]*v2[1];
}
typedef int bool;
#define true 1
#define false 0
// define Neutron object
typedef struct neutron_t {
double r[2];
double v[2]; // position and velocity
bool absorbed;
bool escaped; // final fate
} Neutron;
void initNeutron(Neutron * self) { // constructor initalizes neutron
self->r[0] = 2;
self->r[1] = 0;
double theta = 2 * pi * rand();
self->v[0] = (v_0 * cos(theta));
self->v[1] = (v_0 * sin(theta));
self->absorbed = false;
self->escaped = false;
}
void move_to_a(Neutron * self) { // move to inner wall of shield
double r2 = dotProduct(self->r, self->r);
double v2 = dotProduct(self->v, self->v);
double rv = dotProduct(self->r, self->v);
double t = sqrt((a*a - r2)/v2 + rv*rv/v2/v2) - rv/v2;
self->r[0] += t * self->v[0];
self->r[1] += t * self->v[1];
}
void move_free(Neutron * self) { // free motion in shield
double r2 = dotProduct(self->r, self->r);
double v2 = dotProduct(self->v, self->v);
double rv = dotProduct(self->r, self->v);
double t_max;
if (rv < 0 && abs(rv) >= v2 * (r2 - a*a))
t_max = abs(rv)/v2 - sqrt(rv*rv/v2/v2 - (r2 - a*a)/v2);
else
t_max = sqrt((b*b - r2)/v2 + rv*rv/v2/v2) - rv/v2;
double t = - tau * log(rand());
if (t > t_max)
t = t_max;
self->r[0] += t * self->v[0];
self->r[1] += t * self->v[1];
}
void interact(Neutron * self) { // interaction with nucleus in shield
if (rand() < p_a) { // absorbed by nucleus
self->absorbed = true;
} else { // scatter from nucleus
double v_mag = f * sqrt(dotProduct(self->v, self->v));
double theta = 2 * pi * rand();
self->v[0] = v_mag * cos(theta);
self->v[1] = v_mag * sin(theta);
}
}
bool step_succeeded(Neutron * self) { // attempt a Monte Carlo transport step
if (self->absorbed || self->escaped)
return false;
double r_mag = sqrt(dotProduct(self->r, self->r));
if (r_mag > b) { // outside the shield
self->escaped = true;
return false;
}
if (r_mag < a) // inside core region
move_to_a(self);
move_free(self); // inside the shield
interact(self);
return true;
}
int main() {
printf( "Neutron Reactor Shielding Monte Carlo\n");
printf( " -------------------------------------\n");
printf( " Enter number of particles to simulate:\n ");
int n;
scanf("%d", &n);
// set random number generator seed
srand(time(NULL));
//printf( " Using ".rg.get_algorithm());
//printf( " and seed ". rg.get_seed());
FILE *file;
file = fopen("neutron.data","w");
int step_sum = 0;
int absorbed = 0;
int escaped = 0;
int i = 0;
for (i; i < n; i++) {
Neutron * neutron;
initNeutron(neutron);
int steps = 0;
do {
fprintf(file,"%f\t%f\n",neutron->r[0],neutron->r[1]);
++steps;
} while (step_succeeded(neutron));
fprintf(file,"\n");
step_sum += steps;
if (neutron->absorbed)
++absorbed;
if (neutron->escaped)
++escaped;
}
fclose(file);
double averageSteps = (double)step_sum/(double)n;
printf( " Number escaped = %i\n", escaped);
printf( " Number absorbed = %i\n", absorbed);
printf( " Averge number of steps = %f\n", averageSteps);
printf( " %i trajectories in file neutron.data\n", n);
}