Intensity function integral in python

Question

There is a function that determines the intensity of the Fraunhofer diffraction pattern of a circular aperture ... ( additional information )

The integral of the function at a distance x = [-3.8317, 3.8317] should be about 83.8% (assuming I0 is 100), and when you increase the distance to [-13.33, 13.33], it should be about 95%. But when I use integral in python, the answer is wrong. I do not know what is going on in my code :(

from scipy.integrate import quad
from scipy import special as sp
I0=100.0
dist=3.8317
I= quad(lambda x:( I0*((2*sp.j1(x)/x)**2))  , -dist, dist)[0]
print I

The result of the integral cannot be more than 100 (I0), because it is the diffraction of I0 ... I don’t know .. maybe scaling ... maybe a method! :(

Answer 1

, -, . , :

scipy.integrate.quad , . 0 (, !), .

- . . , , .

, , . , , , (/) . , .

2 pi r ( x r ). :

f = lambda(r): r*(sp.j1(r)/r)**2

f = lambda(r): sp.j1(r)**2/r

:

f = lambda(r): r * (sp.j0(r) + sp.jn(2,r))

, - . (. !).

( , .) ( ):

fullpower = quad(f, 1e-9, np.inf)[0]

:

pwr = quad(f, 1e-9, 3.8317)[0] / fullpower

0.839 ( 84%). (13.33):

pwr = quad(f, 1e-9, 13.33)

0,954.

, , 1e-9 0. , . 1-9 1-12, . , , , 1-30, . ( , .)

:

import matplotlib.pyplot as plt
import numpy as np

x = linspace(0.01, 20, 1000)
intg = np.array([ quad(f, 1e-9, xx)[0] for xx in x])

plt.plot(x, intg/fullpower)
plt.grid('on')
plt.show()

:

, , Airy .

---

: I0 ( , , W/m2), ( /2, W). 100 . .

, :

P (x) = P0 (1 - J0 (x) ^ 2 - J1 (x) ^ 2),

P0 - .

Answer 2

, , Sympy:

import sympy as sy

sy.init_printing()  # LaTeX like pretty printing in IPython

x,d = sy.symbols("x,d", real=True)

I0=100
dist=3.8317
f = I0*((2*sy.besselj(1,x)/x)**2)  # the integrand
F = f.integrate((x, -d, d))  # symbolic integration
print(F.evalf(subs={d:dist}))  # numeric evalution

F :

1600*d*besselj(0, Abs(d))**2/3 + 1600*d*besselj(1, Abs(d))**2/3 - 800*besselj(1, Abs(d))**2/(3*d)

besselj(0,r), sp.j0(r).

Answer 3

jacobian x = 0. points ":

f = lambda x:( I0*((2*sp.j1(x)/x)**2))
I = quad(f, -dist, dist, points = [0])

( ?)

331.4990321315221