Линии дифракционной решетки в Python ⇐ Python

[Anonymous]

Я пытаюсь воспроизвести рисунок дифракционной решетки (рис., описанный ниже) с помощью Python. Однако у меня с этим проблемы. Вот краткое изложение теории и проблемы, с которой я столкнулся:
Интенсивность дифракционной картины на экране, расстояние $x$ от центральной оси системы, определяется по формуле:
Формула интенсивности:


где u — расстояние от центральной оси (пунктирная линия на рисунке ниже)
Схема в задаче:


Мне нужно воспроизвести фигуру, имеющую центральную светлую бахрому и чередующиеся светлые и темные полосы (с уменьшением интенсивности) слева и справа от центральной полосы. На данный момент у меня есть следующий фрагмент кода:

Код: Выделить всё

import numpy as np
from scipy.integrate import trapz

def Intensity_img(u,x):
wavelength = 500*1e-3 # micro meters
f = 1e+6 # micro meters
alpha = np.pi/20
w = 1000

return np.cos(np.pi*u/w) * np.sin(alpha*u)*np.sin(2*np.pi*x*u/(wavelength*f))

x = np.arange(-50000,50100,100)
y = np.arange(-50000,50100,100)
# X, Y = np.meshgrid(x,y)
print(len(x))
a = -500
b = 500
N = 1000
u = np.arange(a,b+1,1)
print(len(u))
I_img = np.full((len(x),len(x)),0)
I_img = np.float128(I_img)

for i in range(len(x)-1):
y_img = Intensity_img(u,x[i])

I_img[i] = trapz(y_img,u)

I = I_img**2

Действительная часть функции интенсивности интегрируется до 0 (поскольку действительная часть представляет собой нечетную функцию на симметричном интервале). Мнимая часть должна быть интегрирована, и я использовал метод catchz из библиотеки scipy.
Когда я создаю график плотности этой функции, я получаю нечто совсем другое:

Две горизонтальные яркие полосы и центральная (горизонтальная) темная полоса


Будем очень признательны за любые полезные советы и подсказки. Я думаю, что делаю что-то не так в цикле for.

Подробнее здесь: https://stackoverflow.com/questions/737 ... -in-python