Плотностный график с питоном, создающим дифракционную картину с интегралами Бесселя, но он не прекращает работать

Я пытаюсь сделать круговую дифракционную картину, который имеет центральное место, окруженное серией колец. Для этого используется Интеграл Бесселя, который определен в коде.

Мои проблемы в том, что это занимает слишком много времени, как я ждал 10 минут для запуска кода, но не получил ничего для отображения. Я понимаю, что это потому, что мой Интеграл Бесселя имеет 1000 итераций в точке может ли кто-нибудь помочь с этим ?

На правильном ли я пути?

Я пытаюсь себя научите себя python и вычислительной физике через книгу Марка Ньюманса Computational Physics упражнение 5.4 вычислительной физики.Вот ссылка на эту главу. Это на странице 9. http://www-personal.umich.edu/~mejn / cp / chapters / int. pdf

Вот образ, который я пытаюсь создать.

концентрическое кольцо.

Мой Код: 

import numpy as np
import pylab as plt
import math as mathy

#N = number of slicices 
#h = b-a/N 

def J(m,x): #Bessel Integral
    def f(theta):
        return (1/mathy.pi)*mathy.cos(m*theta - x*mathy.sin(theta)) #I replaced np. with mathy. for this line

    N = 1000
    A = 0
    a=0
    b=mathy.pi
    h = ((b-a)/N)

    for k in range(1,N,2):

        A +=  4*f(a + h*k)

    for k in range(2,N,2):

        A +=  2*f(a + h*k)

    Idx =  (h/3)*(A + f(a)+f(b))

    return Idx

def I(lmda,r): #Intensity
    k = (mathy.pi/lmda)    
    return ((J(1,k*r))/(k*r))**2

wavelength = .5        # microm meters
I0 = 1
points = 500           
sepration = 0.2  

Intensity = np.empty([points,points],np.float)

for i in range(points):
    y = sepration*i
    for j in range(points):
        x = sepration*j
        r = np.sqrt((x)**2+(y)**2)

        if r < 0.000000000001:
            Intensity[i,j]= 0.5 #this is the lim as r  -> 0, I -> 0.5
        else: 
            Intensity[i,j] = I0*I(wavelength,r)

plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()
1 2
python matplotlib density-plot imshow

1 ответ:

Ваш код, кажется, работает нормально-если я уменьшу N до 100 (с 1000) и размер изображения (points) до 50 (С 500). Примерно через 4 часа времени выполнения я получаю следующее изображение:

Введите описание изображения здесь

Вот что мы получаем при профилировании вашего кода с помощью cProfile: 

$ python -m cProfile -s time bessel.py | head -n 10
         361598 function calls (359660 primitive calls) in 3.470 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   252399    2.250    0.000    2.250    0.000 bessel.py:24(f)
     2499    0.821    0.000    3.079    0.001 bessel.py:23(J)
        1    0.027    0.027    3.472    3.472 bessel.py:15(<module>)
     2499    0.015    0.000    3.093    0.001 bessel.py:45(I)
        1    0.013    0.013    0.013    0.013 backend_macosx.py:1(<module>)

Таким образом, похоже, что большая часть вашего времени выполнения тратится в f. Вы можете оптимизировать эту функцию или попробовать запустить код с помощью PyPy. Он является отличной оптимизации такого рода вещи. Вам необходимо установить их версию numpy (см. http://pypy.readthedocs.org/en/latest/getting-started.html#). но PyPy завершает ваш исходный код в 40-х годах на моей машине (с удалением графического материала).

EDIT

У меня нет plotlib, установленного в PyPy в моей системе, поэтому я заменил ваши вызовы plt в конце на 

#plt.imshow(Intensity,vmax=0.01,cmap='hot')
#plt.gray()
#plt.show()
np.savetxt('bessel.txt', Intensity)

И создал отдельную программу, которую я выполняю с обычным Python, содержащую: 

import numpy as np
import pylab as plt

Intensity = np.loadtxt('bessel.txt')

plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()

Это создает изображение ниже со следующими изменениями в коде: 

sepration = 0.008 # decrease separation

Intensity = np.empty([points,points],np.float)

for i in range(points):
    y = sepration*(i - points/2) # centre on image
    for j in range(points):
        x = sepration*(j - points/2)

Введите описание изображения здесь

2