Generating isotropic steps with lengths distributed as exp(-x/lambda)

Question

I'm writing a code about radiation transportation and have generated isotropic vectors as shown by my code. I'm now unsure about how to generate isotropic steps with lengths distributed as exp(-x/lambda). Any suggestions?

def isotropic_unit_vectors():
nparticles = 1000
lambda_a = 45
lambda_s = 0.3

sigma_a = 1/lambda_a
sigma_s = 1/lambda_s
sigma_T = sigma_a + sigma_s

lambda_T = 1/sigma_T
u = np.random.uniform(low = 0, high =1, size = 1000)
step = -lambda_T*np.log(u)
theta = 2*np.pi*np.random.rand(nparticles)
phi = np.arccos(2*np.random.rand(nparticles)-1)

dx = step*np.cos(theta)*np.cos(phi)
dy = step*np.cos(theta)*np.sin(phi)
dz = step*np.sin(theta)

rand = np.column_stack((dx,dy,dz))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(dx, dy, dz, c='r', marker='.')

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
return rand

edit: I want multiple particles to start at 0,0,0 then move randomly.

Answer 1

Ok, you have right bits and pieces, but they are not quite right together. Here what it might look like, Python 3.7, Anaconda WIndows 10 x64

import numpy as np

def isotropic_unit_vectors():
    """Generates troika of directional cosines isotropically on the sphere"""

    cs_theta = 2.0 * np.random.random() - 1.0
    sn_theta = np.sqrt((1.0 - cs_theta)*(1.0 + cs_theta))
    phi      = 2.0 * np.pi * np.random.random()

    wx = sn_theta * np.cos(phi)
    wy = sn_theta * np.sin(phi)
    wz = cs_theta

    return (wx, wy, wz)

def random_step(?):
    u = 1.0 - np.random.random() # NB! avoids log(0)
    return -? * np.log(u)

nparticles = 100

?_a = 45   # absorbtion
?_s = 0.3  # scattering

# cross-sections
?_a = 1/?_a
?_s = 1/?_s
?_t = ?_a + ?_s

?_t = 1/?_t

trajectories = {} # dictionary for all trajectories

for k in range(0, nparticles):

    trajectory = []

    # initialize particle position
    x = 0.0
    y = 0.0
    z = 0.0

    print(f"Partile number {k+1}")
    print(f"    Position {x}, {y}, {z}")

    trajectory.append((x, y, z))

    while True:
        (wx, wy, wz) = isotropic_unit_vectors() # sample direction

        step = random_step(?_t)

        # make step and find new position
        x += wx * step
        y += wy * step
        z += wz * step

        print(f"    Position {x}, {y}, {z}")

        trajectory.append((x, y, z))

        # simulate absorbtion
        if np.random.random() < ?_a/?_t:
            trajectories[k] = trajectory
            break # absorbtion, end of trajectory, next particle would be tracked
        else:
            pass # scattering, continue with this trajectory

Note, code is NOT vectorized, could be made better and faster with NumPy vectors.

UPDATE

Added code to keep all trajectories in the dictionary indexed by #. Trajectory itself is list of tuples, each tuple contains 3 coordinates of the walk

UPDATE II

If you add code after you got trajectories dictionary, it will plot one trajectory with vertices. Vertix color is defined by distance from (0,0,0).

from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt

fig = plt.figure()
ax = plt.axes(projection='3d')

# Data for a three-dimensional line
k = 7 # print particle #8    

trj = trajectories[k]

xline = [x for x, y, z in trj]
yline = [y for x, y, z in trj]
zline = [z for x, y, z in trj]

# plot trajectory
ax.plot3D(xline, yline, zline, 'gray')

# plot vertices
xpts = [x for x, y, z in trj]
ypts = [y for x, y, z in trj]
zpts = [z for x, y, z in trj]
rpts = [np.sqrt(x*x+y*y+z*z) for x, y, z in trj]

ax.scatter3D(xpts, ypts, zpts, c=rpts, cmap='hsv')

plt.show()

Answer 2

I assume you want to generate random number with CDF; F(x) = exp(-x/lambda), for valid support of x such that 0 <= F(x) <= 1.

Fact: CDF of any random number is uniformly distributed as uniform(0,1) (the CDF Inverse method). That is, generate a uniform(0,1) random number and set it equal to the CDF, written in terms of x. Then solve for x (See step 2).

Method You can follow a process as follows:

1. Generate a random number uniformly between (0, 1), let's call is u1.
2. Set u1 = exp(-x/lambda) => x = -lambda*log(u1)
3. Save x in a list or something.
4. Go back to step 1 until desired number of random numbers are generated.

Note: If exp(-x/lambda) is the PDF and not the CDF, you can find the CDF by integrating the PDF and do the same steps above, with a change in step 2, depending on what the CDF is. For step 1, you can generate a unif(0,1) in Py using random.random().

Hope this helps.

Answer 3

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def isotropic_step_generator(l, size):
    # l is the attenuation length (it's the l in the distribution exp(-x/l))
    
    theta = theta = np.arccos(1-2*np.random.uniform(0, 1, size))
    phi = np.random.uniform(0, 2*np.pi, size)
    radius = -l * np.log(np.random.uniform(size=size))
    
    x = radius * np.sin(theta) * np.cos(phi)
    y = radius * np.sin(theta) * np.sin(phi)
    z = radius * np.cos(theta)
    
    return x, y, z

empty_array = np.empty((1000,3))

i=0
while i<1000:
    empty_array[i] = isotropic_step_generator(45, 1) # attenuation length for water is 45 cm
    i+=1

x_array = empty_array[:, 0]
y_array = empty_array[:, 1]
z_array = empty_array[:, 2]

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x_array, y_array, z_array, s=7)
plt.show()

The above code generates isotropic steps with lengths distributed as exp(-x/lambda) and plots them for water. Note that you must use arccos for theta to prevent bunching of points towards the poles.