Why does my defined function not accept an array as an input?

Question

I am working with a double integral, and the inner integral has variable bounds. I wrote a function, using SciPy's quad integration, that allows me to evaluate this integral. However, what I want is to just evaluate the inner integral, so that all I am left with is a single, non-evaluated integral with respect to some variable. I want to then plot this "half-way" evaluated double integral vs. a range of that variable so that I can see some trend. However, when I input an array of that variable (it is just 0-10000 but increments of 1, it provides the following error message:

"ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"

I have defined functions before that allow me to input an array (which outputs that function at however many points were in the array), so I am unsure as to why this message is popping up now. I think it has something to do with me using SciPy's "quad" integration when defining the function. How do I go around this so that I can input an array?

import numpy as np
from scipy import integrate
from scipy.integrate import quad
import matplotlib.pyplot as plt


#This is the array of values of the variable I ultimately want to have the function plotted 
against
timearray = np.arange(0,10000,1)

#This below defines the general function, that is with respect to two variables (p and t)
def thomtest(p,t):
    d = 3.086e22
    c = 2.998e10
    return (3*((np.cos(p + (2*t/(d*p))))**2))/(8*(t+((p**2)*d/(2*c))))


#The function below evaluates just the inner-integral
def phib(t):
    d = 3.086e22
    c = 2.998e10
    return quad(thomtest,0.00001*(c*t)/d,np.pi, args=(t))[0]


#This evaluates the outer-integral, giving me the complete numerical answer
doubleintegral = quad(phib,0,((3.086e22)/(2.998e10)))


#This below is what gives me the error message: "ValueError: The truth 
#value of an array with more than one element is ambiguous.
#Use a.any() or a.all()". Apparently I cannot input an array
print(phib(timearray))

Attached is a picture of the full error message

Answer 1

In [2]: doubleintegral                                                                                               
Out[2]: (0.9892936902920587, 1.3643899787751934e-08)
In [3]: phib(0)                                                                                                      
Out[3]: 6.377997354641736e-10
In [4]: phib(np.arange(0,5))                                                                                         
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-89b7b858c2f5> in <module>
----> 1 phib(np.arange(0,5))

<ipython-input-1-4c64e69e4f62> in phib(t)
     10     d = 3.086e22
     11     c = 2.998e10
---> 12     return quad(thomtest,0.00001*(c*t)/d,np.pi, args=(t))[0]
     13 
     14 doubleintegral = quad(phib,0,((3.086e22)/(2.998e10)))

/usr/local/lib/python3.6/dist-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    336 
    337     # check the limits of integration: \int_a^b, expect a < b
--> 338     flip, a, b = b < a, min(a, b), max(a, b)
    339 
    340     if weight is None:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
In [5]:       

The error comes when it test the integration bounds, making sure that the first is smaller than the second.

But with an array t, that test is multivalued, and cannot be used:

In [6]: d = 3.086e22  
   ...: c = 2.998e10                                                                                                 
In [7]: 0.00001*(c*np.arange(5))/d                                                                                   
Out[7]: 
array([0.00000000e+00, 9.71484122e-18, 1.94296824e-17, 2.91445237e-17,
       3.88593649e-17])

With quad you have to run your code, phib and doubleintegral once for each element of t. In numpy we try to avoid iterations, but quad only works with scalar bounds. So good old iteration is required.