Monte Carlo Python Plotting Question

[clope023]

If this is better served in a section like Nuclear Engineering please let me know.

Homework Statement

Simulate the projection of a radiation source onto the plane of a collimated detector. Check if said projection is circular.

The attempt at a solution

What I want to do is create a list of vector objects that will be 'moved' from their initial position to the final position at the plane of the detector. The source point is at the origin. The collimator entrance is 2mm (in the y-direction) away and the detector plane is 4mm (in the y-direction) away.
I want the final position of the vectors to be randomly distributed in the plane and check if the x and z elements of those vectors satisfy the equation of the circle with the radius of that circle being the radius of the circular collimator and then plot those positions. I want to be using vpython's vector class for some of the vector operations though if simple nested lists are just as good I will use those too. This is the code I have so far.

import numpy as np
import random
import matplotlib.pyplot as plt
from visual import vector
r = 2.5 #radius of collimator in mm
l = 50.8 #length of collimator in mm
n = 1000 #number of test points
lst = []
for i in range(n):
lst.append(vector(0,0,0)) #creates 1000 vectors at the origin
lst1 = []
rr = r*random.random()
for j in lst:
lst1.append(j + vector(rr, l, rr))
def toList(v):
return [v.x, v.y, v.z] #returns the components of the vector and places them in a list
lst2 = []
for k in lst1:
lst2.append(toList(k)) #places those list of vector components into another list
lst3 = []
for u in lst2:
lst3.append(u.pop())
# pop function removes and returns the last element in the list

At this point I have a list filled with a 1000 3 element lists, and I want to organize the x (0th elements) and z (2nd elements) components and put them in a 2d plot and then check if they satisfy the equation of the circle like so:

plt.plot(x,z,'o')
plt.show()
lst4 = []
lst5 = []
for v in lst3:
if v[0]**2 + v[2]**2 <= r**2:
lst4.append(v)
if v[0]**2 + v[2]**2 > r**2:
lst5.append(v)

So right now this where my logic takes me, am I on the right track? Are my commands right to do what I want to do? I appreciate any and all help, thanks very much.

 

[ilh]

Yes, you're on the right track. You have the basic steps outlined for what you need to do, so now it's just a matter of implementing them in code. For example, you can use the numpy package to quickly calculate the squared distance between each vector and the origin, which will determine whether or not it falls within the circle. You can also use matplotlib to plot your results.