Q:

Integrate the black body spectrum to get the sun's bolometric luminosity

I tried to integrate the black body spectrum (function BBS) to get the sun's bolometric luminosity (Lbol in main), which should about 3.85*10**26 watts. But I get only 1/3 of that.
import numpy as np
from scipy.integrate import quad

global h, c, k # ISU
h = 6.62607e-34
c = 2.998e8
k = 1.38065e-23

global mu_min, mu_max
mu_min, mu_max = 3e10, 3e18
# hertz, corresponds to 1 ångström to 1e8 ångström
# while the sun's spectrum peak at 5000 ångström

global Rsun
Rsun = 6.955e8 # meter

def BBS(mu, tempe):
    i = 2.*h/(c**2.) * (mu**3.) / (np.exp(h/k*mu/tempe)-1.)
    return i

def Teff2Lbol(Teff):
    I = quad(BBS, mu_min, mu_max, args=(Teff,))[0]
    return I

def main():
    T = 5800 # Kelvin
    Lbol = Teff2Lbol(T) * (4*np.pi*Rsun**2.)

A:

Your code is correct but your physics is definitely not. The spectral radiance is measured in units of W m-2 Hz-1 sr-1. The sr-1 is because the radiance is per unit of solid angle and you have to integrate over the whole hemisphere that covers the emitting point. When computing that integral, one has to keep in mind that black bodies are Lambertian, i.e. they emit according to the cosine law: I(theta) = I0*cos(theta), where theta is the angle between the surface normal and the direction of radiation.
To get the total radiance per unit of surface area, you have to multiply the integral over the frequencies by the integral over the upper hemisphere of cos(theta) dOmega (in spherical coordinates). It is easy to compute that integral analytically and its value is exactly pi. Therefore, you have to redefine Teff2Lbol as:
def Teff2Lbol(Teff):
    I = quad(BBS, mu_min, mu_max, args=(Teff,))[0]
    return np.pi*I

Also, note that you are integrating over 8 orders of magnitude on the frequency scale when 98% of the radiated energy actually lies somewhere between 1013 Hz and 1015 Hz. Fortunately, QUADPACK is a very good integrator, capable of dealing with such cases.

