Генерация распределения интенсивности рассеянного света от точечной частицы с центром на сфере с использованием питона

Question

Я пытаюсь смоделировать частицу в трехмерном пространстве, рассеивая свет с распределением интенсивности, пропорциональным 1+cos(theta)^2, где theta измеряется от плоскости xy. Частица является точечной, и поэтому задача состоит в том, чтобы вычислить массив точек (которые будут преобразованы в векторы положения, действующие как лучи света), которые распределены с плотностью, следующей за распределением 1 + cos (тета) ^ 2.

Моей первой мыслью было использование алгоритма, который равномерно распределяет точки на сфере, а затем настраивает его так, чтобы он соответствовал распределению 1 + cos (theta) ^ 2. Однако наиболее точный из этих методов использует случайные величины для определения точек. Поскольку я хочу иметь возможность включать произвольное количество точек, я бы предпочел, чтобы последний метод не использовал случайную выборку (поскольку небольшая случайная выборка приведет к плохому распределению).

Мне кажется, что это очень стандартная проблема, хотя я не смог найти ничего, что имело бы к ней непосредственное отношение.

Вот код для чего-то, что я собрал вместе (я не опытный кодер, так что извините за плохую форму):

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


initial_vector = np.array([0,0,-1])

def phi_rotation(phi):
    #  rotation matrix about y axis  anticlockwise with y towards observer  
    return np.array([[np.cos(phi), 0, np.sin(phi)], [0, 1, 0], [-np.sin(phi), 0, np.cos(phi)]])

def theta_rotation(theta):
    #  rotation matrix about x axis clockwise with x pointing towards (not conventional theta from z axis)
    return np.array([[1, 0, 0], [0, np.cos(theta), -np.sin(theta)], [0, np.sin(theta), np.cos(theta)]])


def generate_equi_intensity_ring(theta, param):
    vectors = []
    first_vector = theta_rotation(theta).dot(np.array([0,0,1]))
    dphi = param/(abs(np.cos(theta))*np.sqrt(1+np.cos(theta)**2))

    for phi in np.linspace(0, 2*np.pi-dphi, 2*np.pi/dphi):
        vectors.append(phi_rotation(phi).dot(first_vector))

    return vectors    

def generate_full_array(param):
    theta = -np.pi/2
    thetas = [theta]
    while theta < np.pi/2:
        dtheta = param/np.sqrt(1+np.cos(theta)**2)
        theta += dtheta
        thetas.append(theta)

    equi_intensitiy_arrays = []

    for theta in thetas:
        vectors = generate_equi_intensity_ring(theta, param)  
        if len(vectors) != 0:
            equi_intensitiy_arrays.append(vectors)    

    return np.concatenate(equi_intensitiy_arrays)



vectors = generate_full_array(0.08) 
x, y, z = zip(*vectors)
print(len(x)) # prints number of rays created    

##### PLOTTING COMPARISON OF DENSITY WITH THETA TO 1+COS(THETA)^2 ####
from scipy import stats

xyz = np.vstack([x,y,z])
density = stats.gaussian_kde(xyz)(xyz)     
idx = density.argsort()
x, y ,z = np.array(x), np.array(y) ,np.array(z) 
x, y, z, density = x[idx], y[idx], z[idx], density[idx]

thetas = np.arccos(y)
plt.plot(np.concatenate([thetas, thetas+np.pi]), 2/max(density) * np.concatenate([density for i in range(2)]), ".")
plt.plot(np.arange(0, 2*np.pi, 0.1), [1+np.sin(t)**2 for t in np.arange(0, 2*np.pi, 0.1)], linewidth=2)

Я масштабировал угловое разделение векторов вокруг оси y, чтобы масштабировать как 1/cos(тета), а также 1/(интенсивность) = 1/(1+cos(тета)^2), это сохраняет плотность точек равномерной и останавливает группирование на полюсах. Это также следует из моего определения тета выше и размера телесного угла, созданного dphi и dtheta.

Однако, когда я сравниваю результат с 1 + cos (тета)**2, это не совсем верно: вычисленная плотность точек (синий) и 1 + cos (тета)^2 (зеленый). ось у = плотность, ось х = тета (0,2pi)

(Примечание: я поменял местами мое определение тэты, измеренное от плоскости xy до измерения от оси z (стандарт), для удобства в некоторых местах, так что cos'становятся грехами и т. Д.)

Ясно, что метод не совсем верный, так как вычисленная плотность не разделена (я думаю, что полюса не одинаково сгенерированы относительно друг друга) и неправильно масштабируется.

Любые предложения будут наиболее ценными!

Спасибо