import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.visualization import quantity_support
from synphot.blackbody import BlackBody1D
from synphot.units import FLAM

bb = BlackBody1D(temperature=5778*u.K)
wav = np.arange(1000, 110000) * u.AA
flux = bb(wav).to(FLAM, u.spectral_density(wav))

with quantity_support():
    plt.figure()
    plt.semilogx(wav, flux)
    plt.axvline(bb.lambda_max.to_value(u.AA), ls='--')
    plt.show()