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; /* 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 * self) { /* initalize 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) { /* 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) { 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 */ move_to_a(self); move_free(self); /* inside the shield */ interact(self); 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; 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); }
Private
[
?
]
Run code
Submit