Python astropy.units.g() Examples

def test_complex_fractional_rounding_errors():
    # See #3788

    kappa = 0.34 * u.cm**2 / u.g
    r_0 = 886221439924.7849 * u.cm
    q = 1.75
    rho_0 = 5e-10 * u.solMass / u.solRad**3
    y = 0.5
    beta = 0.19047619047619049
    a = 0.47619047619047628
    m_h = 1e6*u.solMass

    t1 = 2 * c.c / (kappa * np.sqrt(np.pi))
    t2 = (r_0**-q) / (rho_0 * y * beta * (a * c.G * m_h)**0.5)

    result = ((t1 * t2)**-0.8)

    assert result.unit.physical_type == 'length'
    result.to(u.solRad) 

def _EdS_age(self, z):
        """ Age of the universe in Gyr at redshift ``z``.

        For Omega_radiation = 0 (T_CMB = 0; massless neutrinos)
        the age can be directly calculated as an elliptic integral.
        See, e.g.,
            Thomas and Kantowski, arXiv:0003463

        Parameters
        ----------
        z : array_like
          Input redshifts.

        Returns
        -------
        t : `~astropy.units.Quantity`
          The age of the universe in Gyr at each input redshift.
        """
        if isiterable(z):
            z = np.asarray(z)

        return (2./3) * self._hubble_time * (1+z)**(-3./2) 

def lookback_distance(self, z):
        """
        The lookback distance is the light travel time distance to a given
        redshift. It is simply c * lookback_time.  It may be used to calculate
        the proper distance between two redshifts, e.g. for the mean free path
        to ionizing radiation.

        Parameters
        ----------
        z : array_like
          Input redshifts.  Must be 1D or scalar

        Returns
        -------
        d : `~astropy.units.Quantity`
          Lookback distance in Mpc
        """
        return (self.lookback_time(z) * const.c).to(u.Mpc) 

def test_critical_density():
    from astropy.constants import codata2014

    # WMAP7 but with Omega_relativistic = 0
    # These tests will fail if astropy.const starts returning non-mks
    #  units by default; see the comment at the top of core.py.
    # critical_density0 is inversely proportional to G.
    tcos = core.FlatLambdaCDM(70.4, 0.272, Tcmb0=0.0)
    fac = (const.G / codata2014.G).to(u.dimensionless_unscaled).value
    assert allclose(tcos.critical_density0 * fac,
                    9.309668456020899e-30 * (u.g / u.cm**3))
    assert allclose(tcos.critical_density0,
                    tcos.critical_density(0))
    assert allclose(
        tcos.critical_density([1, 5]) * fac,
        [2.70352772e-29, 5.53739080e-28] * (u.g / u.cm**3))
    assert allclose(
        tcos.critical_density([1., 5.]) * fac,
        [2.70352772e-29, 5.53739080e-28] * (u.g / u.cm**3)) 

def test_units():
    class WithUnits(Container):
        inverse_length = Field(5 / u.m, "foo")
        time = Field(1 * u.s, "bar", unit=u.s)
        grammage = Field(2 * u.g / u.cm ** 2, "baz", unit=u.g / u.cm ** 2)

    c = WithUnits()

    with tempfile.NamedTemporaryFile() as f:
        with HDF5TableWriter(f.name, "data") as writer:
            writer.write("units", c)

        with tables.open_file(f.name, "r") as f:

            assert f.root.data.units.attrs["inverse_length_UNIT"] == "m-1"
            assert f.root.data.units.attrs["time_UNIT"] == "s"
            assert f.root.data.units.attrs["grammage_UNIT"] == "cm-2 g" 

def test_units():
    plt.figure()

    with quantity_support():
        buff = io.BytesIO()

        plt.plot([1, 2, 3] * u.m, [3, 4, 5] * u.kg, label='label')
        plt.plot([105, 210, 315] * u.cm, [3050, 3025, 3010] * u.g)
        plt.legend()
        # Also test fill_between, which requires actual conversion to ndarray
        # with numpy >=1.10 (#4654).
        plt.fill_between([1, 3] * u.m, [3, 5] * u.kg, [3050, 3010] * u.g)
        plt.savefig(buff, format='svg')

        assert plt.gca().xaxis.get_units() == u.m
        assert plt.gca().yaxis.get_units() == u.kg 

def test_quantity_subclass():
    """Check that subclasses are recognized.

    This sadly is not done by matplotlib.units itself, though
    there is a PR to change it:
    https://github.com/matplotlib/matplotlib/pull/13536
    """
    plt.figure()

    with quantity_support():
        plt.scatter(Angle([1, 2, 3], u.deg), [3, 4, 5] * u.kg)
        plt.scatter([105, 210, 315] * u.arcsec, [3050, 3025, 3010] * u.g)
        plt.plot(Angle([105, 210, 315], u.arcsec), [3050, 3025, 3010] * u.g)

        assert plt.gca().xaxis.get_units() == u.deg
        assert plt.gca().yaxis.get_units() == u.kg 

def test_compose_fractional_powers():
    # Warning: with a complicated unit, this test becomes very slow;
    # e.g., x = (u.kg / u.s ** 3 * u.au ** 2.5 / u.yr ** 0.5 / u.sr ** 2)
    # takes 3 s
    x = u.m ** 0.5 / u.yr ** 1.5

    factored = x.compose()

    for unit in factored:
        assert x.decompose() == unit.decompose()

    factored = x.compose(units=u.cgs)

    for unit in factored:
        assert x.decompose() == unit.decompose()

    factored = x.compose(units=u.si)

    for unit in factored:
        assert x.decompose() == unit.decompose() 

def test_comparison_valid_units(self, ufunc):
        q_i1 = np.array([-3.3, 2.1, 10.2]) * u.kg / u.s
        q_i2 = np.array([10., -5., 1.e6]) * u.g / u.Ms
        q_o = ufunc(q_i1, q_i2)
        assert not isinstance(q_o, u.Quantity)
        assert q_o.dtype == bool
        assert np.all(q_o == ufunc(q_i1.value, q_i2.to_value(q_i1.unit)))
        q_o2 = ufunc(q_i1 / q_i2, 2.)
        assert not isinstance(q_o2, u.Quantity)
        assert q_o2.dtype == bool
        assert np.all(q_o2 == ufunc((q_i1 / q_i2)
                                    .to_value(u.dimensionless_unscaled), 2.))
        # comparison with 0., inf, nan is OK even for dimensional quantities
        # (though ignore numpy runtime warnings for comparisons with nan).
        with catch_warnings(RuntimeWarning):
            for arbitrary_unit_value in (0., np.inf, np.nan):
                ufunc(q_i1, arbitrary_unit_value)
                ufunc(q_i1, arbitrary_unit_value*np.ones(len(q_i1)))
            # and just for completeness
            ufunc(q_i1, np.array([0., np.inf, np.nan])) 

def test_subclass_conversion(self, flu_unit, tlu_unit, physical_unit):
        """Check various LogUnit subclasses are equivalent and convertible
        to each other if they correspond to equivalent physical units."""
        values = np.linspace(0., 10., 6)
        flu = flu_unit(physical_unit)

        tlu = tlu_unit(physical_unit)
        assert flu.is_equivalent(tlu)
        assert_allclose(flu.to(tlu), flu.function_unit.to(tlu.function_unit))
        assert_allclose(flu.to(tlu, values),
                        values * flu.function_unit.to(tlu.function_unit))

        tlu2 = tlu_unit(u.Unit(100.*physical_unit))
        assert flu.is_equivalent(tlu2)
        # Check that we round-trip.
        assert_allclose(flu.to(tlu2, tlu2.to(flu, values)), values, atol=1.e-15)

        tlu3 = tlu_unit(physical_unit.to_system(u.si)[0])
        assert flu.is_equivalent(tlu3)
        assert_allclose(flu.to(tlu3, tlu3.to(flu, values)), values, atol=1.e-15)

        tlu4 = tlu_unit(u.g)
        assert not flu.is_equivalent(tlu4)
        with pytest.raises(u.UnitsError):
            flu.to(tlu4, values) 

def gammaw_approx(self, f, P, rho, T):
        rp = P / 1013
        rt = 288 / (T)
        eta1 = 0.955 * rp * rt**0.68 + 0.006 * rho
        eta2 = 0.735 * rp * rt**0.50 + 0.0353 * rt**4 * rho

        def g(f, fi): return 1 + ((f - fi) / (f + fi))**2
        gammaw = (
            (3.98 * eta1 * np.exp(2.23 * (1 - rt))) /
            ((f - 22.235) ** 2 + 9.42 * eta1 ** 2) * g(f, 22.0) +
            (11.96 * eta1 * np.exp(0.70 * (1 - rt))) /
            ((f - 183.310) ** 2 + 11.14 * eta1 ** 2) +
            (0.081 * eta1 * np.exp(6.44 * (1 - rt))) /
            ((f - 321.226) ** 2 + 6.29 * eta1 ** 2) +
            (3.660 * eta1 * np.exp(1.60 * (1 - rt))) /
            ((f - 325.153) ** 2 + 9.22 * eta1 ** 2) +
            (25.37 * eta1 * np.exp(1.09 * (1 - rt))) / ((f - 380.000) ** 2) +
            (17.40 * eta1 * np.exp(1.46 * (1 - rt))) / ((f - 448.000) ** 2) +
            (844.6 * eta1 * np.exp(0.17 * (1 - rt))) / ((f - 557.000) ** 2) *
            g(f, 557.0) + (290.0 * eta1 * np.exp(0.41 * (1 - rt))) /
            ((f - 752.000) ** 2) * g(f, 752.0) +
            (8.3328e4 * eta2 * np.exp(0.99 * (1 - rt))) /
            ((f - 1780.00) ** 2) *
            g(f, 1780.0)) * f ** 2 * rt ** 2.5 * rho * 1e-4
        return gammaw 

def test_invariant_twoarg_array(self, ufunc):

        q_i1 = np.array([-3.3, 2.1, 10.2]) * u.kg / u.s
        q_i2 = np.array([10., -5., 1.e6]) * u.g / u.us
        q_o = ufunc(q_i1, q_i2)
        assert isinstance(q_o, u.Quantity)
        assert q_o.unit == q_i1.unit
        assert_allclose(q_o.value, ufunc(q_i1.value, q_i2.to_value(q_i1.unit))) 

def test_molar_mass_amu():
    x = 1 * (u.g/u.mol)
    y = 1 * u.u
    assert_allclose(x.to_value(u.u, u.molar_mass_amu()), y.value)
    assert_allclose(y.to_value(u.g/u.mol, u.molar_mass_amu()), x.value)
    with pytest.raises(u.UnitsError):
        x.to(u.u) 

def test_comparison_ufuncs_inplace(self, ufunc):
        q_i1 = np.array([-3.3, 2.1, 10.2]) * u.kg / u.s
        q_i2 = np.array([10., -5., 1.e6]) * u.g / u.Ms
        check = np.empty(q_i1.shape, bool)
        ufunc(q_i1.value, q_i2.to_value(q_i1.unit), out=check)

        result = np.empty(q_i1.shape, bool)
        q_o = ufunc(q_i1, q_i2, out=result)
        assert q_o is result
        assert type(q_o) is np.ndarray
        assert q_o.dtype == bool
        assert np.all(q_o == check) 

def test_decompose_bases():
    """
    From issue #576
    """

    from astropy.units import cgs
    from astropy.constants import e

    d = e.esu.unit.decompose(bases=cgs.bases)
    assert d._bases == [u.cm, u.g, u.s]
    assert d._powers == [Fraction(3, 2), 0.5, -1]
    assert d._scale == 1.0 

def test_unit_decomposition(self):
        lu = u.mag(u.Jy)
        assert lu.decompose() == u.mag(u.Jy.decompose())
        assert lu.decompose().physical_unit.bases == [u.kg, u.s]
        assert lu.si == u.mag(u.Jy.si)
        assert lu.si.physical_unit.bases == [u.kg, u.s]
        assert lu.cgs == u.mag(u.Jy.cgs)
        assert lu.cgs.physical_unit.bases == [u.g, u.s] 

def test_quantity_decomposition():
    lq = 10.*u.mag(u.Jy)
    assert lq.decompose() == lq
    assert lq.decompose().unit.physical_unit.bases == [u.kg, u.s]
    assert lq.si == lq
    assert lq.si.unit.physical_unit.bases == [u.kg, u.s]
    assert lq.cgs == lq
    assert lq.cgs.unit.physical_unit.bases == [u.g, u.s] 

def _float_or_none(x, digits=3):
    """ Helper function to format a variable that can be a float or None"""
    if x is None:
        return str(x)
    fmtstr = "{0:.{digits}g}".format(x, digits=digits)
    return fmtstr.format(x) 

def test_dimensionless_angles():
    # test that the angles_dimensionless option allows one to change
    # by any order in radian in the unit (#1161)
    rad1 = u.dimensionless_angles()
    assert u.radian.to(1, equivalencies=rad1) == 1.
    assert u.deg.to(1, equivalencies=rad1) == u.deg.to(u.rad)
    assert u.steradian.to(1, equivalencies=rad1) == 1.
    assert u.dimensionless_unscaled.to(u.steradian, equivalencies=rad1) == 1.
    # now quantities
    assert (1.*u.radian).to_value(1, equivalencies=rad1) == 1.
    assert (1.*u.deg).to_value(1, equivalencies=rad1) == u.deg.to(u.rad)
    assert (1.*u.steradian).to_value(1, equivalencies=rad1) == 1.
    # more complicated example
    I = 1.e45 * u.g * u.cm**2  # noqa
    Omega = u.cycle / (1.*u.s)
    Erot = 0.5 * I * Omega**2
    # check that equivalency makes this work
    Erot_in_erg1 = Erot.to(u.erg, equivalencies=rad1)
    # and check that value is correct
    assert_allclose(Erot_in_erg1.value, (Erot/u.radian**2).to_value(u.erg))

    # test build-in equivalency in subclass
    class MyRad1(u.Quantity):
        _equivalencies = rad1

    phase = MyRad1(1., u.cycle)
    assert phase.to_value(1) == u.cycle.to(u.radian) 

def _flat_age(self, z):
        """ Age of the universe in Gyr at redshift ``z``.

        For Omega_radiation = 0 (T_CMB = 0; massless neutrinos)
        the age can be directly calculated as an elliptic integral.
        See, e.g.,
            Thomas and Kantowski, arXiv:0003463

        Parameters
        ----------
        z : array_like
          Input redshifts.

        Returns
        -------
        t : `~astropy.units.Quantity`
          The age of the universe in Gyr at each input redshift.
        """
        if isiterable(z):
            z = np.asarray(z)

        # Use np.sqrt, np.arcsinh instead of math.sqrt, math.asinh
        # to handle properly the complex numbers for 1 - Om0 < 0
        prefactor = (2./3) * self._hubble_time / \
            np.lib.scimath.sqrt(1 - self._Om0)
        arg = np.arcsinh(np.lib.scimath.sqrt((1 / self._Om0 - 1 + 0j) /
                                             (1 + z)**3))
        return (prefactor * arg).real 

def test_equivalent_units():
    from astropy.units import imperial
    with u.add_enabled_units(imperial):
        units = u.g.find_equivalent_units()
        units_set = set(units)
        match = set(
            [u.M_e, u.M_p, u.g, u.kg, u.solMass, u.t, u.u, u.M_earth,
             u.M_jup, imperial.oz, imperial.lb, imperial.st, imperial.ton,
             imperial.slug])
        assert units_set == match

    r = repr(units)
    assert r.count('\n') == len(units) + 2 

def density(mass, radius, MR_units='earth'):
    """Compute density from mass and radius

    Args:
        mass (float): mass [MR_units]
        radius (float): radius [MR_units]
        MR_units (string): (optional) units of mass and radius. Must be 'earth', or 'jupiter' (default 'earth').

    Returns:
        float: density in g/cc
    """

    mass = np.array(mass)
    radius = np.array(radius)

    if MR_units.lower() == 'earth':
        uradius = u.R_earth
        umass = u.M_earth
    elif MR_units.lower() == 'jupiter':
        uradius = u.R_jup
        umass = u.M_jup
    else:
        raise Exception("MR_units must be 'earth', or 'jupiter'")

    vol = 4. / 3. * np.pi * (radius * uradius) ** 3
    rho = ((mass * umass / vol).to(u.g / u.cm ** 3)).value
    return rho 

def map_wet_term_radio_refractivity(lat, lon, p=50):
    """
    Method to determine the wet term of the radio refractivity


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points


    Returns
    -------
    N_wet: Quantity
        Wet term of the radio refractivity (-)



    References
    ----------
    [1] The radio refractive index: its formula and refractivity data
    https://www.itu.int/rec/R-REC-P.453/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    val = __model.map_wet_term_radio_refractivity(lat, lon, p)
    return prepare_output_array(val, type_output) * u.g / u.m**3 

def zenit_water_vapour_attenuation(
            self, lat, lon, p, f, V_t=None, h=None):
        f_ref = 20.6        # [GHz]
        p_ref = 780         # [hPa]
        if V_t is None:
            V_t = total_water_vapour_content(lat, lon, p, h).value
        rho_ref = V_t / 4     # [g/m3]
        t_ref = 14 * np.log(0.22 * V_t / 4) + 3    # [Celsius]

        gammaw_approx_vect = np.vectorize(self.gammaw_approx)
        return (0.0173 * V_t *
                gammaw_approx_vect(f, p_ref, rho_ref, t_ref + 273) /
                gammaw_approx_vect(f_ref, p_ref, rho_ref, t_ref + 273)) 

def gamma0_approx(f, P, rho, T):
    """
    Method to estimate the specific attenuation due to dry atmosphere using the
    approximate method descibed in Annex 2.

    Parameters
    ----------
    f : number or Quantity
        Frequency (GHz)
    P : number or Quantity
        Atmospheric pressure (hPa)
    rho : number or Quantity
        Water vapor density (g/m3)
    T : number or Quantity
        Absolute temperature (K)


    Returns
    -------
    gamma_w : Quantity
        Dry atmosphere specific attenuation (dB/km)

    References
    --------
    [1] Attenuation by atmospheric gases:
    https://www.itu.int/rec/R-REC-P.676/en
    """
    global __model
    type_output = type(f)
    f = prepare_quantity(f, u.GHz, 'Frequency')
    P = prepare_quantity(P, u.hPa, 'Atmospheric pressure')
    rho = prepare_quantity(rho, u.g / u.m**3, 'Water vapour density')
    T = prepare_quantity(T, u.K, 'Temperature')
    val = __model.gamma0_approx(f, P, rho, T)
    return prepare_output_array(val, type_output) * u.dB / u.km 

def gammaw_exact(f, P, rho, T):
    """
    Method to estimate the specific attenuation due to water vapour using
    the line-by-line method described in Annex 1 of the recommendation.


    Parameters
    ----------
    f : number or Quantity
        Frequency (GHz)
    P : number or Quantity
        Atmospheric pressure (hPa)
    rho : number or Quantity
        Water vapor density (g/m3)
    T : number or Quantity
        Absolute temperature (K)


    Returns
    -------
    gamma_w : Quantity
        Water vapour specific attenuation (dB/km)

    References
    --------
    [1] Attenuation by atmospheric gases:
    https://www.itu.int/rec/R-REC-P.676/en
    """
    global __model
    type_output = type(f)
    f = prepare_quantity(f, u.GHz, 'Frequency')
    P = prepare_quantity(P, u.hPa, 'Atmospheric pressure ')
    rho = prepare_quantity(rho, u.g / u.m**3, 'Water vapour density')
    T = prepare_quantity(T, u.K, 'Temperature')
    val = __model.gammaw_exact(f, P, rho, T)
    return prepare_output_array(val, type_output) * u.dB / u.km 

def gamma0_exact(f, P, rho, T):
    """
    Method to estimate the specific attenuation due to dry atmosphere using
    the line-by-line method described in Annex 1 of the recommendation.

    Parameters
    ----------
    f : number or Quantity
        Frequency (GHz)
    P : number or Quantity
        Atmospheric pressure (hPa)
    rho : number or Quantity
        Water vapor density (g/m3)
    T : number or Quantity
        Absolute temperature (K)


    Returns
    -------
    gamma_w : Quantity
        Dry atmosphere specific attenuation (dB/km)

    References
    --------
    [1] Attenuation by atmospheric gases:
    https://www.itu.int/rec/R-REC-P.676/en
    """
    global __model
    type_output = type(f)
    f = prepare_quantity(f, u.GHz, 'Frequency')
    P = prepare_quantity(P, u.hPa, 'Atmospheric pressure')
    rho = prepare_quantity(rho, u.g / u.m**3, 'Water vapour density')
    T = prepare_quantity(T, u.K, 'Temperature')
    val = __model.gamma0_exact(f, P, rho, T)
    return prepare_output_array(val, type_output) * u.dB / u.km 

def gamma_exact(f, P, rho, T):
    """
    Method to estimate the specific attenuation using the line-by-line method
    described in Annex 1 of the recommendation.


    Parameters
    ----------
    f : number or Quantity
        Frequency (GHz)
    P : number or Quantity
        Atmospheric pressure (hPa)
    rho : number or Quantity
        Water vapor density (g/m3)
    T : number or Quantity
        Absolute temperature (K)

    Returns
    -------
    gamma : Quantity
        Specific attenuation (dB/km)

    References
    --------
    [1] Attenuation by atmospheric gases:
    https://www.itu.int/rec/R-REC-P.676/en
    """
    global __model
    f = prepare_quantity(f, u.GHz, 'Frequency')
    type_output = type(f)

    P = prepare_quantity(P, u.hPa, 'Atmospheric pressure ')
    rho = prepare_quantity(rho, u.g / u.m**3, 'Water vapour density')
    T = prepare_quantity(T, u.K, 'Temperature')
    val = __model.gamma_exact(f, P, rho, T)
    return prepare_output_array(val, type_output) * u.dB / u.km 

def surface_water_vapour_density(lat, lon, p, alt=None):
    """
    Method to compute the surface water vapour density along a path  at any
    desired location on the surface of the Earth.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points
    p : number
        Percentage of time exceeded for p% of the average year
    alt : number, sequence, or numpy.ndarray
        Altitude of the receivers. If None, use the topographical altitude as
        described in recommendation ITU-R P.1511


    Returns
    -------
    rho: Quantity
       Surface water vapour density (g/m3)


    References
    ----------
    [1] Water vapour: surface density and total columnar content
    https://www.itu.int/rec/R-REC-P.836/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    alt = prepare_input_array(alt)
    alt = prepare_quantity(alt, u.km, 'Altitude of the receivers')
    val = __model.surface_water_vapour_density(lat, lon, p, alt)
    return prepare_output_array(val, type_output) * u.g / u.m**3 

def standard_water_vapour_density(h, h_0=2, rho_0=7.5):
    """
    Method to compute the water vapour density of an standard atmosphere at
    a given height.

    The reference standard atmosphere is based on the United States Standard
    Atmosphere, 1976, in which the atmosphere is divided into seven successive
    layers showing linear variation with temperature.


    Parameters
    ----------
    h : number or Quantity
        Height (km)
    h_0 : number or Quantity
        Scale height (km)
    rho_0 : number or Quantity
        Surface water vapour density (g/m^3)


    Returns
    -------
    rho: Quantity
        Water vapour density (g/m^3)


    References
    ----------
    [1] Reference Standard Atmospheres
    https://www.itu.int/rec/R-REC-P.835/en
    """
    global __model

    h = prepare_quantity(h, u.km, 'Height')
    h_0 = prepare_quantity(h_0, u.km, 'Scale height')
    rho_0 = prepare_quantity(
        rho_0,
        u.g / u.m**3,
        'Surface water vapour density')
    return __model.standard_water_vapour_density(h, h_0, rho_0) * u.g / u.m**3