#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
const double pi = 3.1415926; /* pi */
double a = 1; /*inner radius of wall */
double b = 1.5; /* outer radius of wall */
double tau = 0.2; /* mean free time for v=1 */
double p_a = 0.05; /* probability of absorption */
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 * this) { /* initalize neutron */
this->r[0] = 2;
this->r[1] = 0;
double theta = 2 * pi * rand();
this->v[0] = (v_0 * cos(theta));
this->v[1] = (v_0 * sin(theta));
this->absorbed = false;
this->escaped = false;
}
void move_to_a(Neutron * this) { /* move to inner wall of shield*/
double r2 = dotProduct(this->r, this->r);
double v2 = dotProduct(this->v, this->v);
double rv = dotProduct(this->r, this->v);
double t = sqrt((a*a - r2)/v2 + rv*rv/v2/v2) - rv/v2;
this->r[0] += t * this->v[0];
this->r[1] += t * this->v[1];
}
void move_free(Neutron * this) { /* in shield*/
double r2 = dotProduct(this->r, this->r);
double v2 = dotProduct(this->v, this->v);
double rv = dotProduct(this->r, this->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;
this->r[0] += t * this->v[0];
this->r[1] += t * this->v[1];
}
void interact(Neutron * this) { /* interaction with nucleus in shield */
if (rand() < p_a) { /* absorbed by nucleus */
this->absorbed = true;
} else { /* scatter from nucleus */
double v_mag = f * sqrt(dotProduct(this->v, this->v));
double theta = 2 * pi * rand();
this->v[0] = v_mag * cos(theta);
this->v[1] = v_mag * sin(theta);
}
}
bool step_succeeded(Neutron * this) {
if (this->absorbed || this->escaped)
return false;
double r_mag = sqrt(dotProduct(this->r, this->r));
if (r_mag > b) { /* outside the shield */
this->escaped = true;
return false;
}
if (r_mag < a) /* inside core */
move_to_a(this);
move_free(this); /* inside the shield */
interact(this);
return true;
}
int main() {
printf( "Neutron Reactor Shielding simultion\n");
printf( " -------------------------------------\n");
printf( " Number of particles to simulate:\n ");
int n;
scanf("%d", &n);
srand(time(NULL));
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 = malloc( sizeof( Neutron ) );
if( neutron == NULL ) return;
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);
printf( " Number escaped = %i\n", escaped);
printf( " Number absorbed = %i\n", absorbed);
printf( " %i trajectories in file neutron.data\n", n);
}