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> 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] = 0; this->r[1] = 0; double theta = 2 * pi * rand()/RAND_MAX; 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 && fabs(rv) >= v2 * (r2 - a*a)) t_max = fabs(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())/(RAND_MAX); 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()/RAND_MAX < 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())/ RAND_MAX; 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); (double) rand()/RAND_MAX; 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=(Neutron *)malloc(sizeof(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