Python units.year方法代码示例

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def test_numerical_limits(distance):
    """
    Tests the numerical stability of the default settings for the finite
    difference transformation calculation.  This is *known* to fail for at
    >~1kpc, but this may be improved in future versions.
    """
    time = Time('J2017') + np.linspace(-.5, .5, 100)*u.year

    icoo = ICRS(ra=0*u.deg, dec=10*u.deg, distance=distance,
                pm_ra_cosdec=0*u.marcsec/u.yr, pm_dec=0*u.marcsec/u.yr,
                radial_velocity=0*u.km/u.s)
    gcoo = icoo.transform_to(GCRS(obstime=time))
    rv = gcoo.radial_velocity.to('km/s')

    # if its a lot bigger than this - ~the maximal velocity shift along
    # the direction above with a small allowance for noise - finite-difference
    # rounding errors have ruined the calculation
    assert np.ptp(rv) < 65*u.km/u.s 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def test_mission_is_over(self):
        r"""Test mission_is_over method.

        Approach: Allocate time until mission completes.  Check that the mission terminated at
        the right time.
        """
        life = 0.1 * u.year
        tk = self.fixture(missionLife=life.to(u.year).value, missionPortion=1.0)
        sim = self.allmods[0](scriptfile=self.script1)
        allModes = sim.OpticalSystem.observingModes
        Obs = sim.Observatory
        OS = sim.OpticalSystem
        det_mode = list(filter(lambda mode: mode['detectionMode'] == True, allModes))[0]

        # 1) mission not over
        tk.exoplanetObsTime = 0*u.d
        tk.currentTimeAbs = tk.missionStart
        tk.currentTimeNorm = 0*u.d
        self.assertFalse(tk.mission_is_over(OS, Obs, det_mode)) #the mission has just begun

        # 2) exoplanetObsTime exceeded
        tk.exoplanetObsTime = 1.1*tk.missionLife.to('day')*tk.missionPortion # set exoplanetObsTime to failure condition
        self.assertTrue(tk.mission_is_over(OS, Obs, det_mode))
        tk.exoplanetObsTime = 0.*tk.missionLife.to('day')*tk.missionPortion # reset exoplanetObsTime

        # 3) missionLife exceeded
        tk.currentTimeNorm = 1.1*tk.missionLife.to('day')
        tk.currentTimeAbs = tk.missionStart + 1.1*tk.missionLife.to('day')
        self.assertTrue(tk.mission_is_over(OS, Obs, det_mode))
        tk.currentTimeNorm = 0*u.d
        tk.currentTimeAbs = tk.missionStart

        # 4) OBendTimes Exceeded
        tk.OBendTimes = [10]*u.d
        tk.OBnumber = 0
        tk.currentTimeNorm = tk.OBendTimes[tk.OBnumber] + 1*u.d
        tk.currentTimeAbs = tk.missionStart + tk.currentTimeNorm
        self.assertTrue(tk.mission_is_over(OS, Obs, det_mode))
        tk.currentTimeAbs = 0*u.d
        tk.currentTimeAbs = tk.missionStart 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
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 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def _add_gaia_figure_elements(self, tpf, fig, magnitude_limit=18):
        """Make the Gaia Figure Elements"""
        # Get the positions of the Gaia sources
        c1 = SkyCoord(tpf.ra, tpf.dec, frame='icrs', unit='deg')
        # Use pixel scale for query size
        pix_scale = 21.0
        # We are querying with a diameter as the radius, overfilling by 2x.
        from astroquery.vizier import Vizier
        Vizier.ROW_LIMIT = -1
        result = Vizier.query_region(c1, catalog=["I/345/gaia2"],
                                     radius=Angle(np.max(tpf.shape[1:]) * pix_scale, "arcsec"))
        no_targets_found_message = ValueError('Either no sources were found in the query region '
                                              'or Vizier is unavailable')
        too_few_found_message = ValueError('No sources found brighter than {:0.1f}'.format(magnitude_limit))
        if result is None:
            raise no_targets_found_message
        elif len(result) == 0:
            raise too_few_found_message
        result = result["I/345/gaia2"].to_pandas()
        result = result[result.Gmag < magnitude_limit]
        if len(result) == 0:
            raise no_targets_found_message
        radecs = np.vstack([result['RA_ICRS'], result['DE_ICRS']]).T
        coords = tpf.wcs.all_world2pix(radecs, 1) ## TODO, is origin supposed to be zero or one?
        year = ((tpf.astropy_time[0].jd - 2457206.375) * u.day).to(u.year)
        pmra = ((np.nan_to_num(np.asarray(result.pmRA)) * u.milliarcsecond/u.year) * year).to(u.arcsec).value
        pmdec = ((np.nan_to_num(np.asarray(result.pmDE)) * u.milliarcsecond/u.year) * year).to(u.arcsec).value
        result.RA_ICRS += pmra
        result.DE_ICRS += pmdec

        # Gently size the points by their Gaia magnitude
        sizes = 10000.0 / 2**(result['Gmag']/2)

        plt.scatter(coords[:, 0]+tpf.column, coords[:, 1]+tpf.row, c='firebrick', alpha=0.5, edgecolors='r', s=sizes)
        plt.scatter(coords[:, 0]+tpf.column, coords[:, 1]+tpf.row, c='None', edgecolors='r', s=sizes)
        plt.xlim([tpf.column, tpf.column+tpf.shape[1]])
        plt.ylim([tpf.row, tpf.row+tpf.shape[2]])

        return fig 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def semi_amplitude(Msini, P, Mtotal, e, Msini_units='jupiter'):
    """Compute Doppler semi-amplitude

    Args:
        Msini (float): mass of planet [Mjup]
        P (float): Orbital period [days]
        Mtotal (float): Mass of star + mass of planet [Msun]
        e (float): eccentricity
        Msini_units (Optional[str]): Units of Msini {'earth','jupiter'}
            default: 'jupiter'

    Returns:
        Doppler semi-amplitude [m/s]

    """

    # convert inputs to array so they work with units
    P = np.array(P)
    Msini = np.array(Msini)
    Mtotal = np.array(Mtotal)
    e = np.array(e)

    P = (P * u.d).to(u.year).value
    if Msini_units.lower() == 'jupiter':
        pass
    elif Msini_units.lower() == 'earth':
        Msini = (Msini * u.M_earth).to(u.M_jup).value
    else:
        raise Exception("Msini_units must be 'earth', or 'jupiter'")

    K = K_0*(1 - e**2)**-0.5*Msini*P**(-1.0/3.0)*Mtotal**(-2.0 / 3.0)

    return K 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def test_gcrs_diffs():
    time = Time('2017-01-01')
    gf = GCRS(obstime=time)
    sung = get_sun(time)  # should have very little vhelio

    # qtr-year off sun location should be the direction of ~ maximal vhelio
    qtrsung = get_sun(time-.25*u.year)

    # now we use those essentially as directions where the velocities should
    # be either maximal or minimal - with or perpendiculat to Earh's orbit
    msungr = CartesianRepresentation(-sung.cartesian.xyz).represent_as(SphericalRepresentation)
    suni = ICRS(ra=msungr.lon, dec=msungr.lat, distance=100*u.au,
                pm_ra_cosdec=0*u.marcsec/u.yr, pm_dec=0*u.marcsec/u.yr,
                radial_velocity=0*u.km/u.s)
    qtrsuni = ICRS(ra=qtrsung.ra, dec=qtrsung.dec, distance=100*u.au,
                   pm_ra_cosdec=0*u.marcsec/u.yr, pm_dec=0*u.marcsec/u.yr,
                   radial_velocity=0*u.km/u.s)

    # Now we transform those parallel- and perpendicular-to Earth's orbit
    # directions to GCRS, which should shift the velocity to either include
    # the Earth's velocity vector, or not (for parallel and perpendicular,
    # respectively).
    sung = suni.transform_to(gf)
    qtrsung = qtrsuni.transform_to(gf)

    # should be high along the ecliptic-not-sun sun axis and
    # low along the sun axis
    assert np.abs(qtrsung.radial_velocity) > 30*u.km/u.s
    assert np.abs(qtrsung.radial_velocity) < 40*u.km/u.s
    assert np.abs(sung.radial_velocity) < 1*u.km/u.s

    suni2 = sung.transform_to(ICRS)
    assert np.all(np.abs(suni2.data.differentials['s'].d_xyz) < 3e-5*u.km/u.s)
    qtrisun2 = qtrsung.transform_to(ICRS)
    assert np.all(np.abs(qtrisun2.data.differentials['s'].d_xyz) < 3e-5*u.km/u.s) 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def _as_absolute_time_if_needed(self, name, times):
        """
        Convert the provided times to absolute times using the current _tstart
        value, if needed.
        """
        if self._tstart is not None:
            # Some time formats/scales can't represent dates/times too far
            # off from the present, so we need to mask values offset by
            # more than 100,000 yr (the periodogram algorithm can return
            # transit times of e.g 1e300 for some periods).
            reset = np.abs(times.to_value(u.year)) > 100000
            times[reset] = 0
            times = self._tstart + times
            times[reset] = np.nan
        return times 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def __init__(self, missionStart=60634, missionLife=0.1, 
        missionPortion=1, OBduration=np.inf, missionSchedule=None,
        cachedir=None, **specs):

        _outspec = {}
   
        #start the outspec
        self._outspec = {}

        # get cache directory
        self.cachedir = get_cache_dir(cachedir)
        self._outspec['cachedir'] = self.cachedir
        specs['cachedir'] = self.cachedir 

        # load the vprint function (same line in all prototype module constructors)
        self.vprint = vprint(specs.get('verbose', True))
        
        # illegal value checks
        assert missionLife >= 0, "Need missionLife >= 0, got %f"%missionLife
        # arithmetic on missionPortion fails if it is outside the legal range
        assert missionPortion > 0 and missionPortion <= 1, \
                "Require missionPortion in the interval [0,1], got %f"%missionPortion
        # OBduration must be positive nonzero
        assert OBduration*u.d > 0*u.d, "Required OBduration positive nonzero, got %f"%OBduration
        
        # set up state variables
        # tai scale specified because the default, utc, requires accounting for leap
        # seconds, causing warnings from astropy.time when time-deltas are added
        self.missionStart = Time(float(missionStart), format='mjd', scale='tai')#the absolute date of mission start must have scale tai
        self.missionPortion = float(missionPortion)#the portion of missionFinishNorm the instrument can observe for
        
        # set values derived from quantities above
        self.missionLife = float(missionLife)*u.year#the total amount of time since mission start that can elapse MUST BE IN YEAR HERE FOR OUTSPEC
        self.missionFinishAbs = self.missionStart + self.missionLife.to('day')#the absolute time the mission can possibly end
        
        # initialize values updated by functions
        self.currentTimeNorm = 0.*u.day#the current amount of time since mission start that has elapsed
        self.currentTimeAbs = self.missionStart#the absolute mission time
        
        # initialize observing block times arrays. #An Observing Block is a segment of time over which observations may take place
        self.missionSchedule = missionSchedule
        self.init_OB(str(missionSchedule), OBduration*u.d)
        
        # initialize time spend using instrument
        self.exoplanetObsTime = 0*u.day
        
        # populate outspec
        for att in self.__dict__:
            if att not in ['vprint','_outspec']:
                dat = self.__dict__[att]
                self._outspec[att] = dat.value if isinstance(dat,(u.Quantity,Time)) else dat 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def generate_fZ(self, Obs, TL, TK, mode, hashname):
        """Calculates fZ values for all stars over an entire orbit of the sun
            
        Args:
            Obs (module):
                Observatory module
            TL (module):
                Target List Module
            TK (TimeKeeping object):
                TimeKeeping object
            mode (dict):
                Selected observing mode
            hashname (string):
                hashname describing the files specific to the current json script
                
        Updates Attributes:
            fZ_startSaved[1000, TL.nStars] (astropy Quantity array):
                Surface brightness of zodiacal light in units of 1/arcsec2 for each star over 1 year at discrete points defined by resolution
        """
        #Generate cache Name#########################################################
        cachefname = hashname+'starkfZ'

        #Check if file exists########################################################
        if os.path.isfile(cachefname):#check if file exists
            self.vprint("Loading cached fZ from %s"%cachefname)
            try:
                with open(cachefname, "rb") as ff:
                    tmpfZ = pickle.load(ff)
            except UnicodeDecodeError:
                with open(cachefname, "rb") as ff:
                    tmpfZ = pickle.load(ff,encoding='latin1')
            return tmpfZ

        #IF the Completeness vs dMag for Each Star File Does Not Exist, Calculate It
        else:
            self.vprint("Calculating fZ")
            #OS = self.OpticalSystem#Testing to be sure I can remove this
            #WA = OS.WA0#Testing to be sure I can remove this
            sInds= np.arange(TL.nStars)
            startTime = np.zeros(sInds.shape[0])*u.d + TK.currentTimeAbs#Array of current times
            resolution = [j for j in range(1000)]
            fZ = np.zeros([sInds.shape[0], len(resolution)])
            dt = 365.25/len(resolution)*u.d
            for i in xrange(len(resolution)):#iterate through all times of year
                time = startTime + dt*resolution[i]
                fZ[:,i] = self.fZ(Obs, TL, sInds, time, mode)
            
            with open(cachefname, "wb") as fo:
                pickle.dump(fZ,fo)
                self.vprint("Saved cached 1st year fZ to %s"%cachefname)
            return fZ 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
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) 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import year [as 别名]
def test_apply_space_motion():
    # use this 12 year period because it's a multiple of 4 to avoid the quirks
    # of leap years while having 2 leap seconds in it
    t1 = Time('2000-01-01T00:00')
    t2 = Time('2012-01-01T00:00')

    # Check a very simple case first:
    frame = ICRS(ra=10.*u.deg, dec=0*u.deg,
                 distance=10.*u.pc,
                 pm_ra_cosdec=0.1*u.deg/u.yr,
                 pm_dec=0*u.mas/u.yr,
                 radial_velocity=0*u.km/u.s)

    # Cases that should work (just testing input for now):
    c1 = SkyCoord(frame, obstime=t1, pressure=101*u.kPa)
    applied1 = c1.apply_space_motion(new_obstime=t2)
    applied2 = c1.apply_space_motion(dt=12*u.year)

    assert isinstance(applied1.frame, c1.frame.__class__)
    assert isinstance(applied2.frame, c1.frame.__class__)
    assert_allclose(applied1.ra, applied2.ra)
    assert_allclose(applied1.pm_ra, applied2.pm_ra)
    assert_allclose(applied1.dec, applied2.dec)
    assert_allclose(applied1.distance, applied2.distance)

    # ensure any frame attributes that were there before get passed through
    assert applied1.pressure == c1.pressure

    # there were 2 leap seconds between 2000 and 2010, so the difference in
    # the two forms of time evolution should be ~2 sec
    adt = np.abs(applied2.obstime - applied1.obstime)
    assert 1.9*u.second < adt.to(u.second) < 2.1*u.second

    c2 = SkyCoord(frame)
    applied3 = c2.apply_space_motion(dt=6*u.year)
    assert isinstance(applied3.frame, c1.frame.__class__)
    assert applied3.obstime is None

    # this should *not* be .6 deg due to space-motion on a sphere, but it
    # should be fairly close
    assert 0.5*u.deg < applied3.ra-c1.ra < .7*u.deg

    # the two cases should only match somewhat due to it being space motion, but
    # they should be at least this close
    assert quantity_allclose(applied1.ra-c1.ra, (applied3.ra-c1.ra)*2, atol=1e-3*u.deg)
    # but *not* this close
    assert not quantity_allclose(applied1.ra-c1.ra, (applied3.ra-c1.ra)*2, atol=1e-4*u.deg)

    with pytest.raises(ValueError):
        c2.apply_space_motion(new_obstime=t2)