Double Integral with Gauss- Hermite for one component

Question

I am trying to perform the following integral

$$\int_{0}^{2\pi}\int_{0}^{+\infty} \frac{r'\left(e^{-r'^2/2\sigma^2}\right)\left(r-r'\cos(\theta-\theta')\right)}{r^2+r'^2-2rr'\cos(\theta-\theta')}dr'dθ'$$

Using Gauss-Hermite for $r$ and Simpson 1/3 rule for $\theta$ with no success. I can't find my mistake but the output should look like Fig. 2. This was my code (sorry for my bad formatting, this is my first time uploading here).

$\sigma$ should be assumed as 1.

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as ss

def rt(d, r, theta ,sig):
    return r*(d-r*np.cos(theta))*np.exp(-r**2/(2*sig**2))/(d**2+r**2-2*d*r*np.cos(theta))
def intheta1(d, r, b, sig, N):
    h = b/N
    I = rt(d,r,0,sig) + rt(d,r,b,sig)
    for i in range(1, N, 2):
        I += 4*rt(d, r, i*h, sig)
    for j in range(2, N, 2):
        I += 2*rt(d, r, j*h, sig)
    return I*h/3

def intr1(d, b, sig, N, M):
    x, w = ss.roots_hermitenorm(N)
    s = 0
    for k in range(N):
        s += intheta1(d, x[k], b, sig, M)*w[k]
    return s/2

ps = np.linspace(0, 5, 1000)
qs = intr1(xs, 2*np.pi, 1, 1000, 90)

plt.plot(ps, qs)

Answer 1

I have been studying this type of numerical integration and I believe I understood my mistake. First of all I am using Gauss-Hermite which work with limits ${-\infty}$ to ${\infty}$ so using the fact that this function is even makes it so that to integrate from $0$ to ${\infty}$ I have to use np.abs() of my integration variable. Also, using Gauss-Hermite makes it so that I have to remove the exponential function. In this case I am using roots_hermitenorm() so I had to find a way to remove $\exp(-r^2/2)$ from the expression.
I got to these answer which is currently working flawlessly.

python
def integral_theta(r, rline, theta, sigma):
    rline = np.abs(rline)
    return np.exp((-(rline)**2*(1-sigma**2))/2/sigma**2) * rline * (r - rline*np.cos(theta))/(r**2 + rline**2 - 2*r*rline*np.cos(theta))

def i_theta(r, rline, sigma):
    a, b = 0, 2*np.pi
    N = 100
    h = (b-a)/N
    s_odd = 0
    for k in range(1,N,2):
        s_odd += integral_theta(r, rline, a+k*h, sigma)
    s_even = 0
    for j in range(2, N-1,2):
        s_even += integral_theta(r, rline, a+j*h, sigma)

    return h/3*(integral_theta(r, rline, a, sigma) + integral_theta(r, rline, b, sigma) + 4*s_odd + 2*s_even)

def i_r(r, sigma):
    M = 1000
    x, w = ss.roots_hermitenorm(M)
    s = 0
    for h in range(M):
        s += i_theta(r, x[h], sigma)*w[h]
    return s/2/np.pi/sigma**2