'Evaluate the model.
Parameters
x : float, `~numpy.ndarray`, or `~astropy.units.Quantity`
Frequency at which to compute the blackbody. If no units are given,
this defaults to Hz.
temperature : float, `~numpy.ndarray`, or `~astropy.units.Quantity`
Temperature of the blackbody. If no units are given, this defaults
to Kelvin.
bolometric_flux : float, `~numpy.ndarray`, or `~astropy.units.Quantity`
Desired integral for the blackbody.
Returns
y : number or ndarray
Blackbody spectrum. The units are determined from the units of
``bolometric_flux``.'

def evaluate(self, x, temperature, bolometric_flux):

    if isinstance(temperature, u.Quantity):
        temperature = temperature.to(u.K, equivalencies=u.temperature())
    else:
        temperature = u.Quantity(temperature, u.K)
    fnu = (((((np.pi * u.sr) * blackbody_nu(x, temperature)) / const.sigma_sb) / (temperature ** 4)).to((1 / u.Hz)) * bolometric_flux)
    if hasattr(bolometric_flux, u'unit'):
        return fnu
    else:
        return fnu.value
