Python units.rad方法代码示例

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def test_log_occulterResults(self):
        """
        Test that log_occulter_results returns proper dictionary with keys
        """

        atts_list = ['slew_time', 'slew_angle', 'slew_dV', 'slew_mass_used', 'scMass']
        for mod in self.allmods:
            if 'log_occulterResults' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    if 'SotoStarshade' in mod.__name__:
                        obj = mod(f_nStars=4, **copy.deepcopy(self.spec))
                    else:
                        obj = mod(**copy.deepcopy(self.spec))
                DRM = {}
                slewTimes = np.ones((5,))*u.day
                sInds = np.arange(5)
                sd = np.ones((5,))*u.rad
                dV = np.ones((5,))*u.m/u.s

                DRM = obj.log_occulterResults(DRM, slewTimes, sInds, sd, dV)

                for att in atts_list:
                    self.assertTrue(att in DRM, 'Missing key in log_occulterResults for %s' % mod.__name__) 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def test_calc_Phi(self):
        """
        Tests that phase function returns appropriate values.
        """

        for mod in self.allmods:
            if 'calc_Phi' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    obj = mod()
                betas = np.linspace(0.0, np.pi, 100)*u.rad
                Phi = obj.calc_Phi(betas)

                self.assertTrue(len(Phi) == len(betas),'length of phase function values returned does not match input phase angles for %s'%mod.__name__)
                self.assertTrue(np.all(np.isfinite(Phi)),'calc_Phi returned infinite value for %s'%mod.__name__)
                self.assertTrue(np.all(Phi <= 1.0),'calc_Phi returned value > 1 for %s'%mod.__name__)
                self.assertTrue(np.all(Phi >= 0.0),'calc_Phi returned negative value for %s'%mod.__name__) 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def test_outside_IWA_filter(self):
        n0 = self.targetlist.nStars
        #Test default IWA = 0
        self.targetlist.outside_IWA_filter()
        n1 = self.targetlist.nStars

        self.assertEqual( n0, n1 )

        #Test particular case
        self.opticalsystem.IWA = 10 * u.arcsec
        self.targetlist.outside_IWA_filter()
        #assert self.targetlist.nStars == 417 #not a useful test
        n1 = self.targetlist.nStars #reference 
        #introduce two stars with planet below 10 arcsec. should be removed
        self.targetlist.dist[10] = 21 * u.pc #rrange is 1e-3 to 200au, so this planet is below the IWA of 10 arcsec 
        self.targetlist.dist[12] = 22 * u.pc
        self.targetlist.outside_IWA_filter()

        self.assertEqual( self.targetlist.nStars , n1 - 2 )

        #Test limiting case of IWA = PI/2
        self.opticalsystem.IWA = 3.14/2 * u.rad
        with self.assertRaises(IndexError):
            self.targetlist.outside_IWA_filter()
        #self.assertEqual(targetlist.nStars, 0) #Note that nStars is now zero so I can no longer filter out stars. This is why the limiting case of dMagLim = 0 should be done last 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def __init__(self, cachedir=None, **specs):
        
        #start the outspec
        self._outspec = {}

        # 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))
        
        #Define Phase Function Inverse
        betas = np.linspace(start=0.,stop=np.pi,num=1000,endpoint=True)*u.rad
        Phis = self.calc_Phi(betas)
        self.betaFunction = PchipInterpolator(-Phis,betas) #the -Phis ensure the function monotonically increases

        return 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def calc_Phi(self, beta):
        """Calculate the phase function. Prototype method uses the Lambert phase 
        function from Sobolev 1975.
        
        Args:
            beta (astropy Quantity array):
                Planet phase angles at which the phase function is to be calculated,
                in units of rad
                
        Returns:
            Phi (ndarray):
                Planet phase function
        
        """
        
        beta = beta.to('rad').value
        Phi = (np.sin(beta) + (np.pi - beta)*np.cos(beta))/np.pi
        
        return Phi 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def setup_class(self):
        self.impact_reco = ImPACTReconstructor(root_dir=".")
        self.horizon_frame = AltAz()

        self.h1 = HillasParametersContainer(
            x=1 * u.deg,
            y=1 * u.deg,
            r=1 * u.deg,
            phi=Angle(0 * u.rad),
            intensity=100,
            length=0.4 * u.deg,
            width=0.4 * u.deg,
            psi=Angle(0 * u.rad),
            skewness=0,
            kurtosis=0,
        )

    # @pytest.mark.skip('need a dataset for this to work') 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def create_initial_guess(center_x, center_y, radius, telescope_description):
    geometry = telescope_description.camera.geometry
    optics = telescope_description.optics

    focal_length = optics.equivalent_focal_length.to_value(u.m)
    pixel_area = geometry.pix_area[0].to_value(u.m ** 2)
    pixel_radius = np.sqrt(pixel_area / np.pi) / focal_length

    mirror_radius = np.sqrt(optics.mirror_area.to_value(u.m ** 2) / np.pi)

    initial_guess = {}
    initial_guess["impact_parameter"] = mirror_radius / 2
    initial_guess["phi"] = 0
    initial_guess["radius"] = radius.to_value(u.rad)
    initial_guess["center_x"] = center_x.to_value(u.rad)
    initial_guess["center_y"] = center_y.to_value(u.rad)
    initial_guess["ring_width"] = 3 * pixel_radius
    initial_guess["optical_efficiency_muon"] = 0.1

    return initial_guess 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def test_wrap_at():
    a = Angle([-20, 150, 350, 360] * u.deg)
    assert np.all(a.wrap_at(360 * u.deg).degree == np.array([340., 150., 350., 0.]))
    assert np.all(a.wrap_at(Angle(360, unit=u.deg)).degree == np.array([340., 150., 350., 0.]))
    assert np.all(a.wrap_at('360d').degree == np.array([340., 150., 350., 0.]))
    assert np.all(a.wrap_at('180d').degree == np.array([-20., 150., -10., 0.]))
    assert np.all(a.wrap_at(np.pi * u.rad).degree == np.array([-20., 150., -10., 0.]))

    # Test wrapping a scalar Angle
    a = Angle('190d')
    assert a.wrap_at('180d') == Angle('-170d')

    a = Angle(np.arange(-1000.0, 1000.0, 0.125), unit=u.deg)
    for wrap_angle in (270, 0.2, 0.0, 360.0, 500, -2000.125):
        aw = a.wrap_at(wrap_angle * u.deg)
        assert np.all(aw.degree >= wrap_angle - 360.0)
        assert np.all(aw.degree < wrap_angle)

        aw = a.to(u.rad).wrap_at(wrap_angle * u.deg)
        assert np.all(aw.degree >= wrap_angle - 360.0)
        assert np.all(aw.degree < wrap_angle) 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def test_input_units(self):

        self.model._input_units = {'x': u.deg}

        assert_quantity_allclose(self.model(3 * u.deg, 4), 12 * u.deg)
        assert_quantity_allclose(self.model(4 * u.rad, 2), 8 * u.rad)
        assert_quantity_allclose(self.model(4 * u.rad, 2 * u.s), 8 * u.rad * u.s)

        with pytest.raises(UnitsError) as exc:
            self.model(4 * u.s, 3)
        assert exc.value.args[0] == ("MyTestModel: Units of input 'x', s (time), could not be "
                                     "converted to required input units of deg (angle)")
        with pytest.raises(UnitsError) as exc:
            self.model(3, 3)
        assert exc.value.args[0] == ("MyTestModel: Units of input 'x', (dimensionless), could "
                                     "not be converted to required input units of deg (angle)") 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def test_gen_angles(self):
        """
        Test generation of orientation angles.
        
        We expect long. and periapse to be uniformly distributed and
        inclination to be sinusoidally distributed.

        Test method: chi squares
        """

        pp = PlanetPopulation(**self.spec)

        x = 10000
        I, O, w = pp.gen_angles(x)

        #O & w are expected to be uniform
        for param,param_range in zip([O,w],[pp.Orange,pp.wrange]):
            h = np.histogram(param.to('rad').value,100,density=True)
            chi2 = scipy.stats.chisquare(h[0],[1.0/np.diff(param_range.to('rad').value)[0]]*len(h[0]))
            self.assertGreater(chi2[1], 0.95)

        #I is expected to be sinusoidal
        hI = np.histogram(I.to(u.rad).value,100,density=True)
        Ix = np.diff(hI[1])/2.+hI[1][:-1]
        Ip = np.sin(Ix)/2

        Ichi2 = scipy.stats.chisquare(hI[0],Ip)
        assert(Ichi2[1] > 0.95) 

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

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def calculate_slewTimes(self,TL,old_sInd,sInds,sd,obsTimes,currentTime):
        """Finds slew times and separation angles between target stars
        
        This method determines the slew times of an occulter spacecraft needed
        to transfer from one star's line of sight to all others in a given 
        target list.
        
        Args:
            TL (TargetList module):
                TargetList class object
            old_sInd (integer):
                Integer index of the most recently observed star
            sInds (integer ndarray):
                Integer indeces of the star of interest
            sd (astropy Quantity):
                Angular separation between stars in rad
            currentTime (astropy Time):
                Current absolute mission time in MJD
                
        Returns:
            astropy Quantity:
                Time to transfer to new star line of sight in units of days
        """
    
        self.ao = self.thrust/self.scMass
        slewTime_fac = (2.*self.occulterSep/np.abs(self.ao)/(self.defburnPortion/2. - 
            self.defburnPortion**2./4.)).decompose().to('d2')

        if old_sInd is None:
            slewTimes = np.zeros(TL.nStars)*u.d
        else:
            # calculate slew time
            slewTimes = np.sqrt(slewTime_fac*np.sin(abs(sd)/2.)) #an issue exists if sd is negative
            
            #The following are debugging 
            assert np.where(np.isnan(slewTimes))[0].shape[0] == 0, 'At least one slewTime is nan'
        
        return slewTimes 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def log_occulterResults(self,DRM,slewTimes,sInd,sd,dV):
        """Updates the given DRM to include occulter values and results
        
        Args:
            DRM (dict):
                Design Reference Mission, contains the results of one complete
                observation (detection and characterization)
            slewTimes (astropy Quantity):
                Time to transfer to new star line of sight in units of days
            sInd (integer):
                Integer index of the star of interest
            sd (astropy Quantity):
                Angular separation between stars in rad
            dV (astropy Quantity):
                Delta-V used to transfer to new star line of sight in units of m/s
                
        Returns:
            dict:
                Design Reference Mission dictionary, contains the results of one complete
                observation (detection and characterization)
        
        """
        
        DRM['slew_time'] = slewTimes.to('day')
        DRM['slew_angle'] = sd.to('deg')
        
        slew_mass_used = slewTimes*self.defburnPortion*self.flowRate
        DRM['slew_dV'] = (slewTimes*self.ao*self.defburnPortion).to('m/s')
        DRM['slew_mass_used'] = slew_mass_used.to('kg')
        self.scMass = self.scMass - slew_mass_used
        DRM['scMass'] = self.scMass.to('kg')
        
        return DRM 

# 需要导入模块: from astropy import units [as 别名]
# 或者: from astropy.units import rad [as 别名]
def _lonlat_to_healpy(lon, lat, lonlat=False):
    # We use in-place operations below to avoid making temporary arrays - this
    # is safe because the lon/lat arrays returned from healpix_to_lonlat are
    # new and not used elsewhere.
    if lonlat:
        return lon.to(u.deg).value, lat.to(u.deg).value
    else:
        lat, lon = lat.to(u.rad).value, lon.to(u.rad).value
        if np.isscalar(lon):
            return PI_2 - lat, lon
        else:
            lat = np.subtract(PI_2, lat, out=lat)
            return lat, lon