Python astropy.units.yr() Examples

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_valid_quantity_operations1(self):
        """Check adding/substracting/comparing a time-valued quantity works
        with a TimeDelta.  Addition/subtraction should give TimeDelta"""
        t0 = TimeDelta(106400., format='sec')
        q1 = 10.*u.second
        t1 = t0 + q1
        assert isinstance(t1, TimeDelta)
        assert t1.value == t0.value+q1.to_value(u.second)
        q2 = 1.*u.day
        t2 = t0 - q2
        assert isinstance(t2, TimeDelta)
        assert allclose_sec(t2.value, t0.value-q2.to_value(u.second))
        # now comparisons
        assert t0 > q1
        assert t0 < 1.*u.yr
        # and broadcasting
        q3 = np.arange(12.).reshape(4, 3) * u.hour
        t3 = t0 + q3
        assert isinstance(t3, TimeDelta)
        assert t3.shape == q3.shape
        assert allclose_sec(t3.value, t0.value + q3.to_value(u.second)) 

def test_latex_units():
    """
    Check to make sure that Latex and AASTex writers attempt to fall
    back on the **unit** attribute of **Column** if the supplied
    **latexdict** does not specify units.
    """
    t = table.Table([table.Column(name='date', data=['a', 'b']),
               table.Column(name='NUV exp.time', data=[1, 2])])
    latexdict = copy.deepcopy(ascii.latexdicts['AA'])
    latexdict['units'] = {'NUV exp.time': 's'}
    out = StringIO()
    expected = '''\
\\begin{table}{cc}
\\tablehead{\\colhead{date} & \\colhead{NUV exp.time}\\\\ \\colhead{ } & \\colhead{s}}
\\startdata
a & 1 \\\\
b & 2
\\enddata
\\end{table}
'''.replace('\n', os.linesep)

    ascii.write(t, out, format='aastex', latexdict=latexdict)
    assert out.getvalue() == expected
    # use unit attribute instead
    t['NUV exp.time'].unit = u.s
    t['date'].unit = u.yr
    out = StringIO()
    ascii.write(t, out, format='aastex', latexdict=ascii.latexdicts['AA'])
    assert out.getvalue() == expected.replace(
        'colhead{s}', r'colhead{$\mathrm{s}$}').replace(
        'colhead{ }', r'colhead{$\mathrm{yr}$}') 

def haloPosition(self,currentTime):
        """Finds orbit positions of spacecraft in a halo orbit in rotating frame
        
        This method returns the telescope L2 Halo orbit position vector in an ecliptic, 
        rotating frame as dictated by the Circular Restricted Three Body-Problem. 
        The origin of this frame is the centroid of the Sun and Earth-Moon system.
        
        Args:
            currentTime (astropy Time array):
                Current absolute mission time in MJD

        Returns:
            r_halo (astropy Quantity nx3 array):
                Observatory orbit positions vector in an ecliptic, rotating frame 
                in units of AU
        
        """
        
        # Find the time between Earth equinox and current time(s)
        dt = (currentTime - self.equinox).to('yr').value
        t_halo = dt % self.period_halo
        
        # Interpolate to find correct observatory position(s)
        r_halo = self.r_halo_interp_L2(t_halo).T*u.AU
        
        return r_halo 

def haloVelocity(self,currentTime):
        """Finds orbit velocity of spacecraft in a halo orbit in rotating frame
        
        This method returns the telescope L2 Halo orbit velocity vector in an ecliptic, 
        rotating frame as dictated by the Circular Restricted Three Body-Problem. 
        
        Args:
            currentTime (astropy Time array):
                Current absolute mission time in MJD

        Returns:
            v_halo (astropy Quantity nx3 array):
                Observatory orbit velocity vector in an ecliptic, rotating frame 
                in units of AU/year
        
        """
        
        # Find the time between Earth equinox and current time(s)
        
        dt = (currentTime - self.equinox).to('yr').value
        t_halo = dt % self.period_halo
        
        # Interpolate to find correct observatory velocity(-ies)
        v_halo = self.v_halo_interp(t_halo).T
        v_halo = v_halo*u.au/u.year
        
        return v_halo 

def eclip2rot(self,TL,sInd,currentTime):
        """Rotates star position vectors from ecliptic to rotating frame in CRTBP
        
        This method returns a star's position vector in the rotating frame of 
        the Circular Restricted Three Body Problem.  
        
        Args:
            TL (TargetList module):
                TargetList class object
            sInd (integer):
                Integer index of the star of interest
            currentTime (astropy Time):
                Current absolute mission time in MJD

        Returns:
            star_rot (astropy Quantity 1x3 array):
                Star position vector in rotating frame in units of AU
        """
        
        star_pos = TL.starprop(sInd,currentTime).to('au')
        theta    = (np.mod(currentTime.value,self.equinox.value[0])*u.d).to('yr') / u.yr * (2.*np.pi) * u.rad
        
        if currentTime.size == 1:
            star_rot = np.array([np.dot(self.rot(theta, 3), 
                star_pos[x,:].to('AU').value) for x in range(len(star_pos))])[0]*u.AU
        else:
            star_rot = np.array([np.dot(self.rot(theta[x], 3), 
                star_pos[x,:].to('AU').value) for x in range(len(star_pos))])*u.AU

        return star_rot 

def test_valid(self):

            self.aup.a = 3 * u.m
            self.aup.a = [1, 2, 3] * u.cm
            self.aup.a = np.ones((2, 2)) * u.pc

            self.aup.b = 3 * u.m / u.yr
            self.aup.b = [1, 2, 3] * u.cm / u.s
            self.aup.b = np.ones((2, 2)) * u.pc / u.Myr 

def test_compose_no_duplicates():
    new = u.kg / u.s ** 3 * u.au ** 2.5 / u.yr ** 0.5 / u.sr ** 2
    composed = new.compose(units=u.cgs.bases)
    assert len(composed) == 1 

def test_differential_init_errors(self, omit_coslat):
        self._setup(omit_coslat)
        with pytest.raises(u.UnitsError):
            self.USD_cls(0.*u.deg, 10.*u.deg/u.yr) 

def find_known_plans(self):
        """
        Find and return list of known RV stars and list of stars with earthlike planets
        """
        TL = self.TargetList
        SU = self.SimulatedUniverse

        c = 28.4 *u.m/u.s
        Mj = 317.8 * u.earthMass
        Mpj = SU.Mp/Mj                     # planet masses in jupiter mass units
        Ms = TL.MsTrue[SU.plan2star]
        Teff = TL.stellarTeff(SU.plan2star)
        mu = const.G*(SU.Mp + Ms)
        T = (2.*np.pi*np.sqrt(SU.a**3/mu)).to(u.yr)
        e = SU.e

        t_filt = np.where((Teff.value > 3000) & (Teff.value < 6800))[0]    # planets in correct temp range

        K = (c / np.sqrt(1 - e[t_filt])) * Mpj[t_filt] * np.sin(SU.I[t_filt]) * Ms[t_filt]**(-2/3) * T[t_filt]**(-1/3)

        K_filter = (T[t_filt].to(u.d)/10**4).value
        K_filter[np.where(K_filter < 0.03)[0]] = 0.03
        k_filt = t_filt[np.where(K.value > K_filter)[0]]               # planets in the correct K range

        a_filt = k_filt[np.where((SU.a[k_filt] > .95*u.AU) & (SU.a[k_filt] < 1.67*u.AU))[0]]   # planets in habitable zone
        r_filt = a_filt[np.where(SU.Rp.value[a_filt] < 1.75)[0]]                               # rocky planets
        self.earth_candidates = np.union1d(self.earth_candidates, r_filt).astype(int)

        known_stars = np.unique(SU.plan2star[k_filt])
        known_rocky = np.unique(SU.plan2star[r_filt])
        return known_stars.astype(int), known_rocky.astype(int) 

def test_earth_barycentric_velocity_multi_d():
    # Might as well test it with a multidimensional array too.
    t = Time('2016-03-20T12:30:00') + np.arange(8.).reshape(2, 2, 2) * u.yr / 2.
    ep, ev = get_body_barycentric_posvel('earth', t)
    # note: assert_quantity_allclose doesn't like the shape mismatch.
    # this is a problem with np.testing.assert_allclose.
    assert quantity_allclose(ep.get_xyz(xyz_axis=-1),
                             [[-1., 0., 0.], [+1., 0., 0.]]*u.AU,
                             atol=0.06*u.AU)
    expected = u.Quantity([0.*u.one,
                           np.cos(23.5*u.deg),
                           np.sin(23.5*u.deg)]) * ([[-30.], [30.]] * u.km / u.s)
    assert quantity_allclose(ev.get_xyz(xyz_axis=-1), expected,
                             atol=2.*u.km/u.s) 

def test_valid_quantity_input(self):
        """Test Time formats that are allowed to take quantity input."""
        q = 2450000.125*u.day
        t1 = Time(q, format='jd', scale='utc')
        assert t1.value == q.value
        q2 = q.to(u.second)
        t2 = Time(q2, format='jd', scale='utc')
        assert t2.value == q.value == q2.to_value(u.day)
        q3 = q-2400000.5*u.day
        t3 = Time(q3, format='mjd', scale='utc')
        assert t3.value == q3.value
        # test we can deal with two quantity arguments, with different units
        qs = 24.*36.*u.second
        t4 = Time(q3, qs, format='mjd', scale='utc')
        assert t4.value == (q3+qs).to_value(u.day)

        qy = 1990.*u.yr
        ty1 = Time(qy, format='jyear', scale='utc')
        assert ty1.value == qy.value
        ty2 = Time(qy.to(u.day), format='jyear', scale='utc')
        assert ty2.value == qy.value
        qy2 = 10.*u.yr
        tcxc = Time(qy2, format='cxcsec')
        assert tcxc.value == qy2.to_value(u.second)
        tgps = Time(qy2, format='gps')
        assert tgps.value == qy2.to_value(u.second)
        tunix = Time(qy2, format='unix')
        assert tunix.value == qy2.to_value(u.second)
        qd = 2000.*365.*u.day
        tplt = Time(qd, format='plot_date', scale='utc')
        assert tplt.value == qd.value 

def test_no_quantity_input_allowed(self):
        """Time formats that are not allowed to take Quantity input."""
        qy = 1990.*u.yr
        for fmt in ('iso', 'yday', 'datetime', 'byear',
                    'byear_str', 'jyear_str'):
            with pytest.raises(ValueError):
                Time(qy, format=fmt, scale='utc') 

def test_preserve_dtype(self):
        """Test that if an explicit dtype is given, it is used, while if not,
        numbers are converted to float (including decimal.Decimal, which
        numpy converts to an object; closes #1419)
        """
        # If dtype is specified, use it, but if not, convert int, bool to float
        q1 = u.Quantity(12, unit=u.m / u.s, dtype=int)
        assert q1.dtype == int

        q2 = u.Quantity(q1)
        assert q2.dtype == float
        assert q2.value == float(q1.value)
        assert q2.unit == q1.unit

        # but we should preserve any float32 or even float16
        a3_32 = np.array([1., 2.], dtype=np.float32)
        q3_32 = u.Quantity(a3_32, u.yr)
        assert q3_32.dtype == a3_32.dtype
        a3_16 = np.array([1., 2.], dtype=np.float16)
        q3_16 = u.Quantity(a3_16, u.yr)
        assert q3_16.dtype == a3_16.dtype
        # items stored as objects by numpy should be converted to float
        # by default
        q4 = u.Quantity(decimal.Decimal('10.25'), u.m)
        assert q4.dtype == float

        q5 = u.Quantity(decimal.Decimal('10.25'), u.m, dtype=object)
        assert q5.dtype == object 

def test_units_conversion():
    assert_allclose(u.kpc.to(u.Mpc), 0.001)
    assert_allclose(u.Mpc.to(u.kpc), 1000)
    assert_allclose(u.yr.to(u.Myr), 1.e-6)
    assert_allclose(u.AU.to(u.pc), 4.84813681e-6)
    assert_allclose(u.cycle.to(u.rad), 6.283185307179586) 

def test_gravitational_redshift():
    someloc = EarthLocation(lon=-87.7*u.deg, lat=37*u.deg)
    sometime = Time('2017-8-21 18:26:40')
    zg0 = someloc.gravitational_redshift(sometime)

    # should be of order ~few mm/s change per week
    zg_week = someloc.gravitational_redshift(sometime + 7 * u.day)
    assert 1.*u.mm/u.s < abs(zg_week - zg0) < 1*u.cm/u.s

    # ~cm/s over a half-year
    zg_halfyear = someloc.gravitational_redshift(sometime + 0.5 * u.yr)
    assert 1*u.cm/u.s < abs(zg_halfyear - zg0) < 1*u.dm/u.s

    # but when back to the same time in a year, should be tenths of mm
    # even over decades
    zg_year = someloc.gravitational_redshift(sometime - 20 * u.year)
    assert .1*u.mm/u.s < abs(zg_year - zg0) < 1*u.mm/u.s

    # Check mass adjustments.
    # If Jupiter and the moon are ignored, effect should be off by ~ .5 mm/s
    masses = {'sun': constants.G*constants.M_sun,
              'jupiter': 0*constants.G*u.kg,
              'moon': 0*constants.G*u.kg}
    zg_moonjup = someloc.gravitational_redshift(sometime, masses=masses)
    assert .1*u.mm/u.s < abs(zg_moonjup - zg0) < 1*u.mm/u.s
    # Check that simply not including the bodies gives the same result.
    assert zg_moonjup == someloc.gravitational_redshift(sometime,
                                                        bodies=('sun',))
    # And that earth can be given, even not as last argument
    assert zg_moonjup == someloc.gravitational_redshift(
        sometime, bodies=('earth', 'sun',))

    # If the earth is also ignored, effect should be off by ~ 20 cm/s
    # This also tests the conversion of kg to gravitational units.
    masses['earth'] = 0*u.kg
    zg_moonjupearth = someloc.gravitational_redshift(sometime, masses=masses)
    assert 1*u.dm/u.s < abs(zg_moonjupearth - zg0) < 1*u.m/u.s

    # If all masses are zero, redshift should be 0 as well.
    masses['sun'] = 0*u.kg
    assert someloc.gravitational_redshift(sometime, masses=masses) == 0

    with pytest.raises(KeyError):
        someloc.gravitational_redshift(sometime, bodies=('saturn',))

    with pytest.raises(u.UnitsError):
        masses = {'sun': constants.G*constants.M_sun,
                  'jupiter': constants.G*constants.M_jup,
                  'moon': 1*u.km,  # wrong units!
                  'earth': constants.G*constants.M_earth}
        someloc.gravitational_redshift(sometime, masses=masses) 

def test_differential_arithmetic(self, omit_coslat):
        self._setup(omit_coslat)
        s = self.s

        o_lon = self.USD_cls(1.*u.arcsec, 0.*u.arcsec)
        o_lon_by_2 = o_lon / 2.
        assert type(o_lon_by_2) is self.USD_cls
        assert_representation_allclose(o_lon_by_2.to_cartesian(s) * 2.,
                                       o_lon.to_cartesian(s), atol=1e-10*u.one)
        s_lon = s + o_lon
        s_lon2 = s + 2 * o_lon_by_2
        assert type(s_lon) is SphericalRepresentation
        assert_representation_allclose(s_lon, s_lon2, atol=1e-10*u.one)
        o_lon_rec = o_lon_by_2 + o_lon_by_2
        assert type(o_lon_rec) is self.USD_cls
        assert representation_equal(o_lon, o_lon_rec)
        assert_representation_allclose(s + o_lon, s + o_lon_rec,
                                       atol=1e-10*u.one)
        o_lon_0 = o_lon - o_lon
        assert type(o_lon_0) is self.USD_cls
        for c in o_lon_0.components:
            assert np.all(getattr(o_lon_0, c) == 0.)

        o_lon2 = self.USD_cls(1.*u.mas/u.yr, 0.*u.mas/u.yr)
        kks = u.km/u.kpc/u.s
        assert_quantity_allclose(o_lon2.norm(s)[0], 4.74047*kks, atol=1e-4*kks)
        assert_representation_allclose(o_lon2.to_cartesian(s) * 1000.*u.yr,
                                       o_lon.to_cartesian(s), atol=1e-10*u.one)
        s_off = s + o_lon
        s_off2 = s + o_lon2 * 1000.*u.yr
        assert_representation_allclose(s_off, s_off2, atol=1e-10*u.one)

        factor = 1e5 * u.radian/u.arcsec
        if not omit_coslat:
            factor = factor / np.cos(s.lat)
        s_off_big = s + o_lon * factor

        assert_representation_allclose(
            s_off_big, SphericalRepresentation(s.lon + 90.*u.deg,
                                               0.*u.deg, 1e5),
            atol=5.*u.one)

        o_lon3c = CartesianRepresentation(0., 4.74047, 0., unit=kks)
        # This looses information!!
        o_lon3 = self.USD_cls.from_cartesian(o_lon3c, base=s)
        expected0 = self.USD_cls(1.*u.mas/u.yr, 0.*u.mas/u.yr)
        assert_differential_allclose(o_lon3[0], expected0)
        # Part of motion kept.
        part_kept = s.cross(CartesianRepresentation(0, 1, 0, unit=u.one)).norm()
        assert_quantity_allclose(o_lon3.norm(s), 4.74047*part_kept*kks,
                                 atol=1e-10*kks)
        # (lat[0]=0, so works for both normal and CosLat differential)
        s_off_big2 = s + o_lon3 * 1e5 * u.yr * u.radian/u.mas
        expected0 = SphericalRepresentation(90.*u.deg, 0.*u.deg,
                                            1e5*u.one)
        assert_representation_allclose(s_off_big2[0], expected0, atol=5.*u.one) 

def test_differential_arithmetic(self, omit_coslat):
        self._setup(omit_coslat)
        s = self.s

        o_lon = self.SD_cls(1.*u.arcsec, 0.*u.arcsec, 0.*u.kpc)
        o_lon_by_2 = o_lon / 2.
        assert_representation_allclose(o_lon_by_2.to_cartesian(s) * 2.,
                                       o_lon.to_cartesian(s), atol=1e-10*u.kpc)
        assert_representation_allclose(s + o_lon, s + 2 * o_lon_by_2,
                                       atol=1e-10*u.kpc)
        o_lon_rec = o_lon_by_2 + o_lon_by_2
        assert_representation_allclose(s + o_lon, s + o_lon_rec,
                                       atol=1e-10*u.kpc)
        o_lon_0 = o_lon - o_lon
        for c in o_lon_0.components:
            assert np.all(getattr(o_lon_0, c) == 0.)
        o_lon2 = self.SD_cls(1*u.mas/u.yr, 0*u.mas/u.yr, 0*u.km/u.s)
        assert_quantity_allclose(o_lon2.norm(s)[0], 4.74*u.km/u.s,
                                 atol=0.01*u.km/u.s)
        assert_representation_allclose(o_lon2.to_cartesian(s) * 1000.*u.yr,
                                       o_lon.to_cartesian(s), atol=1e-10*u.kpc)
        s_off = s + o_lon
        s_off2 = s + o_lon2 * 1000.*u.yr
        assert_representation_allclose(s_off, s_off2, atol=1e-10*u.kpc)

        factor = 1e5 * u.radian/u.arcsec
        if not omit_coslat:
            factor = factor / np.cos(s.lat)
        s_off_big = s + o_lon * factor

        assert_representation_allclose(
            s_off_big, SphericalRepresentation(s.lon + 90.*u.deg, 0.*u.deg,
                                               1e5*s.distance),
            atol=5.*u.kpc)

        o_lon3c = CartesianRepresentation(0., 4.74047, 0., unit=u.km/u.s)
        o_lon3 = self.SD_cls.from_cartesian(o_lon3c, base=s)
        expected0 = self.SD_cls(1.*u.mas/u.yr, 0.*u.mas/u.yr, 0.*u.km/u.s)
        assert_differential_allclose(o_lon3[0], expected0)
        s_off_big2 = s + o_lon3 * 1e5 * u.yr * u.radian/u.mas
        assert_representation_allclose(
            s_off_big2, SphericalRepresentation(90.*u.deg, 0.*u.deg,
                                                1e5*u.kpc), atol=5.*u.kpc)

        with pytest.raises(TypeError):
            o_lon - s
        with pytest.raises(TypeError):
            s.to_cartesian() + o_lon 

def orbit(self, currentTime, eclip=False):
        """Finds observatory orbit positions vector in heliocentric equatorial (default)
        or ecliptic frame for current time (MJD).
        
        This method returns the telescope L2 Halo orbit position vector.
        
        Args:
            currentTime (astropy Time array):
                Current absolute mission time in MJD
            eclip (boolean):
                Boolean used to switch to heliocentric ecliptic frame. Defaults to 
                False, corresponding to heliocentric equatorial frame.
        
        Returns:
            r_obs (astropy Quantity nx3 array):
                Observatory orbit positions vector in heliocentric equatorial (default)
                or ecliptic frame in units of AU
        
        Note: Use eclip=True to get ecliptic coordinates.
        
        """
        
        # find time from Earth equinox and interpolated position
        dt = (currentTime - self.equinox).to('yr').value
        t_halo = dt % self.period_halo
        r_halo = self.r_halo_interp(t_halo).T
        # find Earth positions in heliocentric ecliptic frame
        r_Earth = self.solarSystem_body_position(currentTime, 'Earth',
                eclip=True).to('AU').value
        # adding Earth-Sun distances (projected in ecliptic plane)
        r_Earth_norm = np.linalg.norm(r_Earth[:,0:2], axis=1)
        r_halo[:,0] = r_halo[:,0] + r_Earth_norm
        # Earth ecliptic longitudes
        lon = np.sign(r_Earth[:,1])*np.arccos(r_Earth[:,0]/r_Earth_norm)
        # observatory positions vector in heliocentric ecliptic frame
        r_obs = np.array([np.dot(self.rot(-lon[x], 3), 
                r_halo[x,:]) for x in range(currentTime.size)])*u.AU
        
        assert np.all(np.isfinite(r_obs)), \
                "Observatory positions vector r_obs has infinite value."
        
        if not eclip:
            # observatory positions vector in heliocentric equatorial frame
            r_obs = self.eclip2equat(r_obs, currentTime)
        
        return r_obs