codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
#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); }
Private
[
?
]
Run code
Submit