-
Notifications
You must be signed in to change notification settings - Fork 2
/
scattering.py
133 lines (107 loc) · 4.03 KB
/
scattering.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 5 15:58:08 2022
@author: Francesco Vascelli
"""
import numpy
from numpy import random
def random_walk(n_end, depth, steps):
"""This method performs a random walk in 2D for a neutron.
Parameters
n_end : destination of the neutron [back, lost, through]
depth : width of the reactor wall in units of mean free path
steps : number of interactions before the neutron stops
"""
#The neutron starts inside the reactor wall towards the positive x-axis.
xy_vector = [0,0,1,0]
for i in range(steps):
p = random.rand()
#The probability of going forward is double with respect to changing direction.
if p<float(1/6):
turn_right(xy_vector)
if p>=float(1/6) and p<float(1/3):
turn_left(xy_vector)
if p>=float(1/3):
go_forward(xy_vector)
#Stop the random walk if the neutron goes all the way through the reactor wall.
if xy_vector[2]==(depth+1):
n_end[2]+=1
break
#Stop the random walk if the neutron returns inside the reactor.
if xy_vector[2]==-1:
n_end[0]+=1
break
#Otherwise, the neutron is lost inside the reactor wall.
else:
n_end[1]+=1
def turn_left(xy_vector):
"""This method changes the direction of the neutron to the left
and takes a step of the random walk.
Parameters
xy_vector : position of the neutron in the xy plane.
The first two values are the position in x and y at the previous step,
the last two are the current position in x and y
"""
#The neutron has +x direction
if xy_vector[2]-xy_vector[0]==1:
xy_vector[0] += 1
xy_vector[3] += 1
#The neutron has -x direction
elif xy_vector[2]-xy_vector[0]==-1:
xy_vector[2] += 1
xy_vector[3] -= 1
#The neutron has +y direction
elif xy_vector[3]-xy_vector[1]==1:
xy_vector[1] += 1
xy_vector[2] -= 1
#The neutron has -y direction
elif xy_vector[3]-xy_vector[1]==-1:
xy_vector[2] += 1
xy_vector[3] += 1
def turn_right(xy_vector):
"""This method changes the direction of the neutron to the right
and takes a step of the random walk.
Parameters
xy_vector : position of the neutron in the xy plane.
The first two numbers are the position in x and y at the previous step,
the last two are the current position in x and y
"""
#The neutron has +x direction
if xy_vector[2]-xy_vector[0]==1:
xy_vector[0] += 1
xy_vector[3] -= 1
#The neutron has -x direction
elif xy_vector[2]-xy_vector[0]==-1:
xy_vector[2] += 1
xy_vector[3] += 1
#The neutron has +y direction
elif xy_vector[3]-xy_vector[1]==1:
xy_vector[1] += 1
xy_vector[2] += 1
#The neutron has -y direction
elif xy_vector[3]-xy_vector[1]==-1:
xy_vector[2] -= 1
xy_vector[3] += 1
def go_forward(xy_vector):
"""This method makes a step of the random walk in the same direction.
Parameters
xy_vector : position of the neutron in the xy plane.
The first two numbers are the position in x and y at the previous step,
the last two are the current position in x and y
"""
#The neutron has +x direction
if xy_vector[2]-xy_vector[0]==1:
xy_vector[0] += 1
xy_vector[2] += 1
#The neutron has -x direction
elif xy_vector[2]-xy_vector[0]==-1:
xy_vector[0] -= 1
xy_vector[2] -= 1
#The neutron has +y direction
elif xy_vector[3]-xy_vector[1]==1:
xy_vector[1] += 1
xy_vector[3] += 1
#The neutron has -y direction
elif xy_vector[3]-xy_vector[1]==-1:
xy_vector[1] -= 1
xy_vector[3] -= 1