human
Write me a the Python program you requested given the following:
Complete the below
def brightness_temperature(beam_area, disp):

    ''''Defines the conversion between Jy/beam and "brightness temperature",
    :math:`T_B`, in Kelvins.  The brightness temperature is a unit very
    commonly used in radio astronomy.  See, e.g., "Tools of Radio Astronomy"
    (Wilson 2009) eqn 8.16 and eqn 8.19 (these pages are available on `google
    books
    <http://books.google.com/books?id=9KHw6R8rQEMC&pg=PA179&source=gbs_toc_r&cad=4#v=onepage&q&f=false>`__).
    :math:`T_B \equiv S_\nu / \left(2 k \nu^2 / c^2 \right)`
    However, the beam area is essential for this computation: the brightness
    temperature is inversely proportional to the beam area
    Parameters
    beam_area : Beam Area equivalent
    Beam area in angular units, i.e. steradian equivalent
    disp : `~astropy.units.Quantity` with spectral units
    The observed `spectral` equivalent `~astropy.units.Unit` (e.g.,
    frequency or wavelength)
    Examples
    Arecibo C-band beam::
    >>> import numpy as np
    >>> from astropy import units as u
    >>> beam_sigma = 50*u.arcsec
    >>> beam_area = 2*np.pi*(beam_sigma)**2
    >>> freq = 5*u.GHz
    >>> equiv = u.brightness_temperature(beam_area, freq)
    >>> u.Jy.to(u.K, equivalencies=equiv)  # doctest: +FLOAT_CMP
    3.526294429423223
    >>> (1*u.Jy).to(u.K, equivalencies=equiv)  # doctest: +FLOAT_CMP
    <Quantity 3.526294429423223 K>
    VLA synthetic beam::
    >>> bmaj = 15*u.arcsec
    >>> bmin = 15*u.arcsec
    >>> fwhm_to_sigma = 1./(8*np.log(2))**0.5
    >>> beam_area = 2.*np.pi*(bmaj*bmin*fwhm_to_sigma**2)
    >>> freq = 5*u.GHz
    >>> equiv = u.brightness_temperature(beam_area, freq)
    >>> u.Jy.to(u.K, equivalencies=equiv)  # doctest: +FLOAT_CMP
    217.2658703625732'
    '''
assistant
#Here is the Python program you requested.
    beam = beam_area.to_value(si.sr)
    nu = disp.to(si.GHz, spectral())
    def convert_Jy_to_K(x_jybm):
        factor = ((((2 * _si.k_B) * si.K) * (nu ** 2)) / (_si.c ** 2)).to_value(astrophys.Jy)
        return ((x_jybm / beam) / factor)
    def convert_K_to_Jy(x_K):
        factor = (astrophys.Jy / (((2 * _si.k_B) * (nu ** 2)) / (_si.c ** 2))).to_value(si.K)
        return ((x_K * beam) / factor)
    return [(astrophys.Jy, si.K, convert_Jy_to_K, convert_K_to_Jy)]