Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from scipy import integrate
- import math
- import numpy as np
- a = 0.250
- s02 = 214.0
- a_s = 0.0163
- def integrand(r, R, s02, a_s, a):
- return 2.0 * r * (r/a)**(-0.1) * (1.0 + (r**2/a**2))**(-2.45)\
- *(math.sqrt(r**2 - R**2))**(-1.0) * (a_s/(1 + (R-0.0283)**2/a_s**2 ))
- def bounds_R(s02, a_s, a):
- return [0, np.inf]
- def bounds_r(R, s02, a_s, a):
- return [R, np.inf]
- result = integrate.nquad(integrand, [bounds_r(R, s02, a_s, a), bounds_R(s02, a_s, a)])
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement