Python astropy.units.mag() Examples

def test_littleh():
    H0_70 = 70*u.km/u.s/u.Mpc
    h70dist = 70 * u.Mpc/u.littleh

    assert_quantity_allclose(h70dist.to(u.Mpc, u.with_H0(H0_70)), 100*u.Mpc)

    # make sure using the default cosmology works
    cosmodist = cosmology.default_cosmology.get().H0.value * u.Mpc/u.littleh
    assert_quantity_allclose(cosmodist.to(u.Mpc, u.with_H0()), 100*u.Mpc)

    # Now try a luminosity scaling
    h1lum = .49 * u.Lsun * u.littleh**-2
    assert_quantity_allclose(h1lum.to(u.Lsun, u.with_H0(H0_70)), 1*u.Lsun)

    # And the trickiest one: magnitudes.  Using H0=10 here for the round numbers
    H0_10 = 10*u.km/u.s/u.Mpc
    # assume the "true" magnitude M = 12.
    # Then M - 5*log_10(h)  = M + 5 = 17
    withlittlehmag = 17 * (u.mag - u.MagUnit(u.littleh**2))
    assert_quantity_allclose(withlittlehmag.to(u.mag, u.with_H0(H0_10)), 12*u.mag) 

def test_ilshift_magnitude(self):
        # test in-place operation and conversion
        mag_fnu_cgs = u.mag(u.erg/u.s/u.cm**2/u.Hz)
        m = np.arange(10.0) * u.mag(u.Jy)
        jy = m.physical
        m2 = m << mag_fnu_cgs
        assert np.all(m2 == m.to(mag_fnu_cgs))
        m2 = m
        m <<= mag_fnu_cgs
        assert m is m2  # Check it was done in-place!
        assert np.all(m.value == m2.value)
        assert m.unit == mag_fnu_cgs
        # Check it works if equivalencies are in-place.
        with u.add_enabled_equivalencies(u.spectral_density(5500*u.AA)):
            st = jy.to(u.ST)
            m <<= u.STmag

        assert m is m2
        assert_quantity_allclose(m.physical, st)
        assert m.unit == u.STmag 

def test_input_units(data):
    t, y, dy, params = data

    t_unit = u.day
    y_unit = u.mag

    with pytest.raises(u.UnitConversionError):
        BoxLeastSquares(t * t_unit, y * y_unit, dy * u.one)
    with pytest.raises(u.UnitConversionError):
        BoxLeastSquares(t * t_unit, y * u.one, dy * y_unit)
    with pytest.raises(u.UnitConversionError):
        BoxLeastSquares(t * t_unit, y, dy * y_unit)
    model = BoxLeastSquares(t*t_unit, y * u.one, dy)
    assert model.dy.unit == model.y.unit
    model = BoxLeastSquares(t*t_unit, y * y_unit, dy)
    assert model.dy.unit == model.y.unit
    model = BoxLeastSquares(t*t_unit, y*y_unit)
    assert model.dy is None 

def test_lshift_magnitude(self):
        mag = 1. << u.ABmag
        assert isinstance(mag, u.Magnitude)
        assert mag.unit == u.ABmag
        assert mag.value == 1.
        # same test for an array, which should produce a view
        a2 = np.arange(10.)
        q2 = a2 << u.ABmag
        assert isinstance(q2, u.Magnitude)
        assert q2.unit == u.ABmag
        assert np.all(q2.value == a2)
        a2[9] = 0.
        assert np.all(q2.value == a2)
        # a different magnitude unit
        mag = 10. << u.STmag
        assert isinstance(mag, u.Magnitude)
        assert mag.unit == u.STmag
        assert mag.value == 10. 

def test_addition_subtraction(self, other):
        """Check physical units are changed appropriately"""
        lu1 = u.mag(u.Jy)
        other_pu = getattr(other, 'physical_unit', u.dimensionless_unscaled)

        lu_sf = lu1 + other
        assert lu_sf.is_equivalent(lu1.physical_unit * other_pu)

        lu_sr = other + lu1
        assert lu_sr.is_equivalent(lu1.physical_unit * other_pu)

        lu_df = lu1 - other
        assert lu_df.is_equivalent(lu1.physical_unit / other_pu)

        lu_dr = other - lu1
        assert lu_dr.is_equivalent(other_pu / lu1.physical_unit) 

def test_model(fit_mean, with_units, freq):
    rand = np.random.RandomState(0)
    t = 10 * rand.rand(40)
    params = 10 * rand.rand(3)

    y = np.zeros_like(t)
    if fit_mean:
        y += params[0]
    y += params[1] * np.sin(2 * np.pi * freq * (t - params[2]))

    if with_units:
        t = t * u.day
        y = y * u.mag
        freq = freq / u.day

    ls = LombScargle(t, y, center_data=False, fit_mean=fit_mean)
    y_fit = ls.model(t, freq)
    assert_quantity_allclose(y_fit, y) 

def test_model_units_mismatch(data):
    t, y, dy = data
    frequency = 1.0
    t_fit = t[:5]

    t = t * u.second
    t_fit = t_fit * u.second
    y = y * u.mag
    frequency = 1.0 / t.unit

    # this should fail because frequency and 1/t units do not match
    with pytest.raises(ValueError) as err:
        LombScargle(t, y).model(t_fit, frequency=1.0)
    assert str(err.value).startswith('Units of frequency not equivalent')

    # this should fail because t and t_fit units do not match
    with pytest.raises(ValueError) as err:
        LombScargle(t, y).model([1, 2], frequency)
    assert str(err.value).startswith('Units of t not equivalent')

    # this should fail because dy and y units do not match
    with pytest.raises(ValueError) as err:
        LombScargle(t, y, dy).model(t_fit, frequency)
    assert str(err.value).startswith('Units of dy not equivalent') 

def test_container_unit_conversion(self, lu_unit):
        """Check that conversion to logarithmic units (u.mag, u.dB, u.dex)
        is only possible when the physical unit is dimensionless."""
        values = np.linspace(0., 10., 6)
        lu1 = lu_unit(u.dimensionless_unscaled)
        assert lu1.is_equivalent(lu1.function_unit)
        assert_allclose(lu1.to(lu1.function_unit, values), values)

        lu2 = lu_unit(u.Jy)
        assert not lu2.is_equivalent(lu2.function_unit)
        with pytest.raises(u.UnitsError):
            lu2.to(lu2.function_unit, values) 

def test_lshift_errors(self):
        m = np.arange(10.0) * u.mag(u.Jy)
        with pytest.raises(u.UnitsError):
            m << u.STmag

        with pytest.raises(u.UnitsError):
            m << u.Jy

        with pytest.raises(u.UnitsError):
            m <<= u.STmag

        with pytest.raises(u.UnitsError):
            m <<= u.Jy 

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_unit_multiple_possible_equivalencies(self):
        lu = u.mag(u.Jy)
        assert lu.is_equivalent(pu_sample) 

def test_raise_to_power(self, power):
        """Check that raising LogUnits to some power is only possible when the
        physical unit is dimensionless, and that conversion is turned off when
        the resulting logarithmic unit (such as mag**2) is incompatible."""
        lu1 = u.mag(u.Jy)

        if power == 0:
            assert lu1 ** power == u.dimensionless_unscaled
        elif power == 1:
            assert lu1 ** power == lu1
        else:
            with pytest.raises(u.UnitsError):
                lu1 ** power

        # With dimensionless, though, it works, but returns a normal unit.
        lu2 = u.mag(u.dimensionless_unscaled)

        t = lu2**power
        if power == 0:
            assert t == u.dimensionless_unscaled
        elif power == 1:
            assert t == lu2
        else:
            assert not isinstance(t, type(lu2))
            assert t == lu2.function_unit**power
            # also check we roundtrip
            t2 = t**(1./power)
            assert t2 == lu2.function_unit
            with u.set_enabled_equivalencies(u.logarithmic()):
                assert_allclose(t2.to(u.dimensionless_unscaled, np.arange(3.)),
                                lu2.to(lu2.physical_unit, np.arange(3.))) 

def test_addition_subtraction_to_normal_units_fails(self, other):
        lu1 = u.mag(u.Jy)
        with pytest.raises(u.UnitsError):
            lu1 + other

        with pytest.raises(u.UnitsError):
            lu1 - other

        with pytest.raises(u.UnitsError):
            other - lu1 

def test_complicated_addition_subtraction(self):
        """for fun, a more complicated example of addition and subtraction"""
        dm0 = u.Unit('DM', 1./(4.*np.pi*(10.*u.pc)**2))
        lu_dm = u.mag(dm0)
        lu_absST = u.STmag - lu_dm
        assert lu_absST.is_equivalent(u.erg/u.s/u.AA) 

def test_neg_pos(self):
        lu1 = u.mag(u.Jy)
        neg_lu = -lu1
        assert neg_lu != lu1
        assert neg_lu.physical_unit == u.Jy**-1
        assert -neg_lu == lu1
        pos_lu = +lu1
        assert pos_lu is not lu1
        assert pos_lu == lu1 

def test_function_values(self, value, unit):
        lq = u.Magnitude(value, unit)
        assert lq == value
        assert lq.unit.function_unit == u.mag
        assert lq.unit.physical_unit == getattr(unit, 'physical_unit',
                                                value.unit.physical_unit) 

def test_indirect_creation(self, unit):
        q1 = 2.5 * unit
        assert isinstance(q1, u.Magnitude)
        assert q1.value == 2.5
        assert q1.unit == unit
        pv = 100. * unit.physical_unit
        q2 = unit * pv
        assert q2.unit == unit
        assert q2.unit.physical_unit == pv.unit
        assert q2.to_value(unit.physical_unit) == 100.
        assert (q2._function_view / u.mag).to_value(1) == -5.
        q3 = unit / 0.4
        assert q3 == q1 

def test_quantity_view(self):
        # Cannot view as Quantity, since the unit cannot be represented.
        with pytest.raises(TypeError):
            self.lq.view(u.Quantity)
        # But a dimensionless one is fine.
        q2 = self.lq2.view(u.Quantity)
        assert q2.unit is u.mag
        assert np.all(q2.value == self.lq2.value)
        lq3 = q2.view(u.Magnitude)
        assert type(lq3.unit) is u.MagUnit
        assert lq3.unit.physical_unit == u.dimensionless_unscaled
        assert np.all(lq3 == self.lq2) 

def test_item_get_and_set(self):
        lq1 = u.Magnitude(np.arange(1., 11.)*u.Jy)
        assert lq1[9] == u.Magnitude(10.*u.Jy)
        lq1[2] = 100.*u.Jy
        assert lq1[2] == u.Magnitude(100.*u.Jy)
        with pytest.raises(u.UnitsError):
            lq1[2] = 100.*u.m
        with pytest.raises(u.UnitsError):
            lq1[2] = 100.*u.mag
        with pytest.raises(u.UnitsError):
            lq1[2] = u.Magnitude(100.*u.m)
        assert lq1[2] == u.Magnitude(100.*u.Jy) 

def test_slice_get_and_set(self):
        lq1 = u.Magnitude(np.arange(1., 10.)*u.Jy)
        lq1[2:4] = 100.*u.Jy
        assert np.all(lq1[2:4] == u.Magnitude(100.*u.Jy))
        with pytest.raises(u.UnitsError):
            lq1[2:4] = 100.*u.m
        with pytest.raises(u.UnitsError):
            lq1[2:4] = 100.*u.mag
        with pytest.raises(u.UnitsError):
            lq1[2:4] = u.Magnitude(100.*u.m)
        assert np.all(lq1[2] == u.Magnitude(100.*u.Jy)) 

def test_raise_to_power(self, power):
        """Check that raising LogQuantities to some power is only possible when
        the physical unit is dimensionless, and that conversion is turned off
        when the resulting logarithmic unit (say, mag**2) is incompatible."""
        lq = u.Magnitude(np.arange(1., 4.)*u.Jy)

        if power == 0:
            assert np.all(lq ** power == 1.)
        elif power == 1:
            assert np.all(lq ** power == lq)
        else:
            with pytest.raises(u.UnitsError):
                lq ** power

        # with dimensionless, it works, but falls back to normal quantity
        # (except for power=1)
        lq2 = u.Magnitude(np.arange(10.))

        t = lq2**power
        if power == 0:
            assert t.unit is u.dimensionless_unscaled
            assert np.all(t.value == 1.)
        elif power == 1:
            assert np.all(t == lq2)
        else:
            assert not isinstance(t, type(lq2))
            assert t.unit == lq2.unit.function_unit ** power
            with u.set_enabled_equivalencies(u.logarithmic()):
                with pytest.raises(u.UnitsError):
                    t.to(u.dimensionless_unscaled) 

def test_complicated_addition_subtraction(self):
        """For fun, a more complicated example of addition and subtraction."""
        dm0 = u.Unit('DM', 1./(4.*np.pi*(10.*u.pc)**2))
        DMmag = u.mag(dm0)
        m_st = 10. * u.STmag
        dm = 5. * DMmag
        M_st = m_st - dm
        assert M_st.unit.is_equivalent(u.erg/u.s/u.AA)
        assert np.abs(M_st.physical /
                      (m_st.physical*4.*np.pi*(100.*u.pc)**2) - 1.) < 1.e-15 

def test_comparison(self):
        lq1 = u.Magnitude(np.arange(1., 4.)*u.Jy)
        lq2 = u.Magnitude(2.*u.Jy)
        assert np.all((lq1 > lq2) == np.array([True, False, False]))
        assert np.all((lq1 == lq2) == np.array([False, True, False]))
        lq3 = u.Dex(2.*u.Jy)
        assert np.all((lq1 > lq3) == np.array([True, False, False]))
        assert np.all((lq1 == lq3) == np.array([False, True, False]))
        lq4 = u.Magnitude(2.*u.m)
        assert not (lq1 == lq4)
        assert lq1 != lq4
        with pytest.raises(u.UnitsError):
            lq1 < lq4
        q5 = 1.5 * u.Jy
        assert np.all((lq1 > q5) == np.array([True, False, False]))
        assert np.all((q5 < lq1) == np.array([True, False, False]))
        with pytest.raises(u.UnitsError):
            lq1 >= 2.*u.m
        with pytest.raises(u.UnitsError):
            lq1 <= lq1.value * u.mag
        # For physically dimensionless, we can compare with the function unit.
        lq6 = u.Magnitude(np.arange(1., 4.))
        fv6 = lq6.value * u.mag
        assert np.all(lq6 == fv6)
        # but not some arbitrary unit, of course.
        with pytest.raises(u.UnitsError):
            lq6 < 2.*u.m 

def setup(self):
        self.mJy = np.arange(1., 5.).reshape(2, 2) * u.mag(u.Jy)
        self.m1 = np.arange(1., 5.5, 0.5).reshape(3, 3) * u.mag()
        self.mags = (self.mJy, self.m1) 

def test_always_ok(self, method):
        for mag in self.mags:
            res = getattr(mag, method)()
            assert np.all(res.value ==
                          getattr(mag._function_view, method)().value)
            if method in ('std', 'ptp', 'diff', 'ediff1d'):
                assert res.unit == u.mag()
            elif method == 'var':
                assert res.unit == u.mag**2
            else:
                assert res.unit == mag.unit 

def test_clip(self):
        for mag in self.mags:
            assert np.all(mag.clip(2. * mag.unit, 4. * mag.unit).value ==
                          mag.value.clip(2., 4.)) 

def test_model_parameters(data, nterms, fit_mean, center_data,
                          errors, with_units):
    if nterms == 0 and not fit_mean:
        return

    t, y, dy = data
    frequency = 1.5
    if with_units:
        t = t * u.day
        y = y * u.mag
        dy = dy * u.mag
        frequency = frequency / t.unit

    if errors == 'none':
        dy = None
    elif errors == 'partial':
        dy = dy[0]
    elif errors == 'full':
        pass
    else:
        raise ValueError(f"Unrecognized error type: '{errors}'")

    ls = LombScargle(t, y, dy,
                     nterms=nterms,
                     fit_mean=fit_mean,
                     center_data=center_data)
    tfit = np.linspace(0, 20, 10)
    if with_units:
        tfit = tfit * u.day

    model = ls.model(tfit, frequency)
    params = ls.model_parameters(frequency)
    design = ls.design_matrix(frequency, t=tfit)
    offset = ls.offset()

    assert len(params) == int(fit_mean) + 2 * nterms

    assert_quantity_allclose(offset + design.dot(params), model) 

def distmod(self):
        """The distance modulus as a `~astropy.units.Quantity`"""
        val = 5. * np.log10(self.to_value(u.pc)) - 5.
        return u.Quantity(val, u.mag, copy=False) 

def _distmod_to_pc(cls, dm):
        dm = u.Quantity(dm, u.mag)
        return cls(10 ** ((dm.value + 5) / 5.), u.pc, copy=False) 

def distmod(self, z):
        """ Distance modulus at redshift ``z``.

        The distance modulus is defined as the (apparent magnitude -
        absolute magnitude) for an object at redshift ``z``.

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

        Returns
        -------
        distmod : `~astropy.units.Quantity`
          Distance modulus at each input redshift, in magnitudes

        See Also
        --------
        z_at_value : Find the redshift corresponding to a distance modulus.
        """

        # Remember that the luminosity distance is in Mpc
        # Abs is necessary because in certain obscure closed cosmologies
        #  the distance modulus can be negative -- which is okay because
        #  it enters as the square.
        val = 5. * np.log10(abs(self.luminosity_distance(z).value)) + 25.0
        return u.Quantity(val, u.mag)