Python Ошибка: Ошибка типа: только массивы размера 1 могут быть преобразованы в Python скаляры

Question

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

from scipy.integrate import nquad
import numpy as np
from scipy.special import genlaguerre, gamma
from scipy.constants import alpha
import matplotlib.pyplot as plt
import cmath

#Constants
epsilon = 13.1 # dielectric constant of the material
gamma_C = 0.5  # donor impurity linewidth 
nr = 3.2       # refractive index of semiconductor
flux = 0       # Phi in eqn 8 magnetic flux
R = 5.0        # radius of the quantum ring in nm
r = np.linspace(0, 6 * R)
rho = r / R
m_0 = 0.0067*0.511 # electron effective mass
h = 4.13e-15       # Planck constant in eV
hbar =  6.58e-16   # reduced Planck constant in eV

#Photon energy
hnu = np.linspace(0, 100) #in eV

#Function that calculates the integrand
def func(rho, theta):
    betai = gama**2/2
    betaf = np.sqrt(1+gama**4/2)
    return R/np.pi*((gama * rho)**(betai + betaf) *
            np.exp(-1/2*(gama * rho)**2) *
          (gama*rho)**2/2  * genlaguerre(1, betaf) *np.cos(theta) * cmath.exp(-1j*theta))

def cross_section(hnu, gama):
    #function that calculates the photoionisation cross section
    betai = np.sqrt( gama**4/4)
    betaf = np.sqrt(1+gama**4/2)
    Ei = gama**2*(1+betai)-gama**4/2
    Ef = gama**2*(3+betaf)-gama**4/2
    return (nr/epsilon * 4*np.pi/3 * alpha * hnu *
            (abs( np.sqrt(1/2**betai*gamma(betai + 1))*
            np.sqrt(gamma(2)/2**betaf*gamma(betaf + 2)) *
            nquad(func, [[0, np.infty], [0, np.pi]])[0])**2 * 
             hbar * gamma_C/(Ef - Ei - hnu)**2 + ( hbar * gamma_C)**2))

#Plot
plt.figure();plt.clf()

for gama in [1.0, 1.5, 2.0]:
    plt.plot(hnu, cross_section(hnu, gama))
plt.legend(['$\gamma = 1.0$', '$\gamma = 1.5$', '$\gamma = 2.0$'] )
plt.ylabel('Photoionization cross\n section $\sigma (10^{-14}cm^{2}$)')
plt.xlabel('Photon energy $h\\nu (meV)$ ')    

Но я получил следующую ошибку:

    runfile('/home/daniel/Arquivos_Python/crosssection.py', wdir='/home/daniel/Arquivos_Python')
Reloaded modules: __mp_main__
Traceback (most recent call last):

  File "/home/daniel/Arquivos_Python/crosssection.py", line 53, in <module>
    plt.plot(hnu, cross_section(hnu, gama))

  File "/home/daniel/Arquivos_Python/crosssection.py", line 47, in cross_section
    hbar * gamma_C/(Ef - Ei - hnu)**2 + ( hbar * gamma_C)**2))

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 805, in nquad
    return _NQuad(func, ranges, opts, full_output).integrate(*args)

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 860, in integrate
    **opt)

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 341, in quad
    points)

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 448, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 860, in integrate
    **opt)

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 341, in quad
    points)

  File "/home/daniel/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 450, in _quad
    return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

TypeError: only size-1 arrays can be converted to Python scalars

<Figure size 432x288 with 0 Axes>

Я искал и есть некоторые вопросы с той же ошибкой на сайте, но для очень простые программы с простыми списками, и ответы не помогают мне понять, как решить эту проблему.