Function integration returns error "only length-1 arrays can be converted to Python scalars"

I'm trying to integrate a function into a given range that shows the flow of particles in a zero-angle (theta) direction as a function of the energy E of the particles. I've tried several ways and got different errors but this is what always persists at the end. My knowledge of Python is limited and I am not understanding how to get the integral.

My role at the moment is like this:

from numpy import radians, cos, arange
from matplotlib.pyplot import plot, yscale, xscale, xlabel, ylabel
from scipy.integrate import quad

def integral(self):
     theta=0
     E = arange(1, 5000, 1)
     costh = cos(radians(theta))
     a = 18 / (E * costh + 145)
     b = (E + 2.7 / costh)**-2.7
     c = (E + 5) / (E + 5 / costh)
     #flux = a*b*c*1**4
     return ((a*b*c*1**4))
     #print(flux)

     """
     plot(flux)
     yscale('log')
     xscale('log')
     xlabel('E')
     ylabel('Fluxo')
     """

#integral(0)

A, err= quad(integral, 500, 1000)

If I remove the self from the function argument I get the error:

"integral () takes 0 positional arguments but 1 was given"

I do not know how to solve.

Does anyone know the solution to this case?