Python astropy.units.rad() Examples

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__) 

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) 

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') 

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)") 

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 

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 

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 

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 

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 

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__) 

def test_one_argument_ufunc_inplace(self, value):
        # without scaling
        s = value * u.rad
        check = s
        np.sin(s, out=s)
        assert check is s
        assert check.unit == u.dimensionless_unscaled
        # with scaling
        s2 = (value * u.rad).to(u.deg)
        check2 = s2
        np.sin(s2, out=s2)
        assert check2 is s2
        assert check2.unit == u.dimensionless_unscaled
        assert_allclose(s.value, s2.value) 

def test_init_array_nocopy(self):

        phi = Angle([8, 9] * u.hourangle)
        theta = Angle([5, 6] * u.deg)
        r = Distance([1, 2] * u.kpc)

        s1 = PhysicsSphericalRepresentation(phi=phi, theta=theta, r=r, copy=False)

        phi[:] = [1, 2] * u.rad
        theta[:] = [3, 4] * u.arcmin
        r[:] = [8, 9] * u.Mpc

        assert_allclose_quantity(phi, s1.phi)
        assert_allclose_quantity(theta, s1.theta)
        assert_allclose_quantity(r, s1.r) 

def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
        # Note that here we need to make sure that x and y are in the same
        # units otherwise this can lead to issues since rotation is not well
        # defined.
        if inputs_unit['x'] != inputs_unit['y']:
            raise UnitsError("Units of 'x' and 'y' inputs should match")
        return {'x_0': inputs_unit['x'],
                'y_0': inputs_unit['x'],
                'r_eff': inputs_unit['x'],
                'theta': u.rad,
                'amplitude': outputs_unit['z']} 

def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
        # Note that here we need to make sure that x and y are in the same
        # units otherwise this can lead to issues since rotation is not well
        # defined.
        if inputs_unit['x'] != inputs_unit['y']:
            raise UnitsError("Units of 'x' and 'y' inputs should match")
        return {'x_0': inputs_unit['x'],
                'y_0': inputs_unit['x'],
                'a': inputs_unit['x'],
                'b': inputs_unit['x'],
                'theta': u.rad,
                'amplitude': outputs_unit['z']} 

def evaluate(x, amplitude, frequency, phase):
        """One dimensional Sine model function"""
        # Note: If frequency and x are quantities, they should normally have
        # inverse units, so that argument ends up being dimensionless. However,
        # np.sin of a dimensionless quantity will crash, so we remove the
        # quantity-ness from argument in this case (another option would be to
        # multiply by * u.rad but this would be slower overall).
        argument = TWOPI * (frequency * x + phase)
        if isinstance(argument, Quantity):
            argument = argument.value
        return amplitude * np.sin(argument) 

def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
        # Note that here we need to make sure that x and y are in the same
        # units otherwise this can lead to issues since rotation is not well
        # defined.
        if inputs_unit['x'] != inputs_unit['y']:
            raise UnitsError("Units of 'x' and 'y' inputs should match")
        return {'x_mean': inputs_unit['x'],
                'y_mean': inputs_unit['x'],
                'x_stddev': inputs_unit['x'],
                'y_stddev': inputs_unit['x'],
                'theta': u.rad,
                'amplitude': outputs_unit['z']} 

def test_unwrap(self):
        q = [0., 3690., -270., 690.] * u.deg
        out = np.unwrap(q)
        expected = (np.unwrap(q.to_value(u.rad)) * u.rad).to(q.unit)
        assert out.unit == expected.unit
        assert np.allclose(out, expected, atol=1*u.urad, rtol=0)
        with pytest.raises(u.UnitsError):
            np.unwrap([1., 2.]*u.m)
        with pytest.raises(u.UnitsError):
            np.unwrap(q, discont=1.*u.m) 

def _to_radian(value):
    """ Convert ``value`` to radian. """
    if isinstance(value, u.Quantity):
        return value.to(u.rad)
    else:
        return np.deg2rad(value) 

def test_mixed_string_and_quantity():
    a1 = Angle(['1d', 1. * u.deg])
    assert_array_equal(a1.value, [1., 1.])
    assert a1.unit == u.deg

    a2 = Angle(['1d', 1 * u.rad * np.pi, '3d'])
    assert_array_equal(a2.value, [1., 180., 3.])
    assert a2.unit == u.deg 

def test_init_array_nocopy(self):

        lon = Longitude([8, 9] * u.hourangle)
        lat = Latitude([5, 6] * u.deg)
        distance = Distance([1, 2] * u.kpc)

        s1 = SphericalRepresentation(lon=lon, lat=lat, distance=distance, copy=False)

        lon[:] = [1, 2] * u.rad
        lat[:] = [3, 4] * u.arcmin
        distance[:] = [8, 9] * u.Mpc

        assert_allclose_quantity(lon, s1.lon)
        assert_allclose_quantity(lat, s1.lat)
        assert_allclose_quantity(distance, s1.distance) 

def test_angle_to_is_angle():
    with pytest.warns(IllegalSecondWarning):
        a = Angle('00:00:60', u.deg)
    assert isinstance(a, Angle)
    assert isinstance(a.to(u.rad), Angle) 

def test_to_string_radian_with_precision():
    """
    Regression test for a bug that caused ``to_string`` to crash for angles in
    radians when specifying the precision.
    """

    # Check that specifying the precision works
    a = Angle(3., unit=u.rad)
    assert a.to_string(precision=3, sep='fromunit') == '3.000rad' 

def test_to_string_formats():
    a = Angle(1.113355, unit=u.deg)
    assert a.to_string(format='latex') == r'$1^\circ06{}^\prime48.078{}^{\prime\prime}$'
    assert a.to_string(format='unicode') == '1°06′48.078″'

    a = Angle(1.113355, unit=u.hour)
    assert a.to_string(format='latex') == r'$1^\mathrm{h}06^\mathrm{m}48.078^\mathrm{s}$'
    assert a.to_string(format='unicode') == '1ʰ06ᵐ48.078ˢ'

    a = Angle(1.113355, unit=u.radian)
    assert a.to_string(format='latex') == r'$1.11336\mathrm{rad}$'
    assert a.to_string(format='unicode') == '1.11336rad' 

def test_icrs_basic(tmpdir):
    wrap_angle = Angle(1.5, unit=units.rad)
    ra = Longitude(25, unit=units.deg, wrap_angle=wrap_angle)
    dec = Latitude(45, unit=units.deg)

    tree = {'coord': ICRS(ra=ra, dec=dec)}

    assert_roundtrip_tree(tree, tmpdir) 

def test_init_array_nocopy(self):

        lon = Longitude([8, 9] * u.hourangle)
        lat = Latitude([5, 6] * u.deg)

        s1 = UnitSphericalRepresentation(lon=lon, lat=lat, copy=False)

        lon[:] = [1, 2] * u.rad
        lat[:] = [3, 4] * u.arcmin

        assert_allclose_quantity(lon, s1.lon)
        assert_allclose_quantity(lat, s1.lat) 

def randomly_sample_sphere(ntosample, randomseed=12345):
    """
    Generates a set of spherical coordinates uniformly distributed over the
    sphere in a way that gives the same answer for the same seed.  Also
    generates a random distance vector on [0, 1] (no units)

    This simply returns (lon, lat, r) instead of a representation to avoid
    failures due to the representation module.
    """
    with NumpyRNGContext(randomseed):
        lat = np.arcsin(np.random.rand(ntosample)*2-1)
        lon = np.random.rand(ntosample)*np.pi*2
        r = np.random.rand(ntosample)

    return lon*u.rad, lat*u.rad, r 

def test_cylindrical(self):

        s = CylindricalRepresentation(rho=[1, 2, 3] * u.pc,
                                      phi=[0., 90., -45.] * u.deg,
                                      z=[3, 4, 5] * u.kpc)
        e = s.unit_vectors()
        self.check_unit_vectors(e)
        sf = s.scale_factors()
        self.check_scale_factors(sf, s)

        s_rho = s + 1. * u.pc * e['rho']
        assert_quantity_allclose(s_rho.rho, s.rho + 1.*u.pc)
        assert_quantity_allclose(s_rho.phi, s.phi)
        assert_quantity_allclose(s_rho.z, s.z)
        s_rho2 = s + 1. * u.pc * sf['rho'] * e['rho']
        assert_representation_allclose(s_rho2, s_rho)

        s_phi = s + s.rho * 1e-5 * e['phi']
        assert_quantity_allclose(s_phi.rho, s.rho)
        assert_quantity_allclose(s_phi.phi, s.phi + 1e-5*u.rad)
        assert_quantity_allclose(s_phi.z, s.z)
        s_phi2 = s + 1e-5 * u.radian * sf['phi'] * e['phi']
        assert_representation_allclose(s_phi2, s_phi)

        s_z = s + 1. * u.pc * e['z']
        assert_quantity_allclose(s_z.rho, s.rho)
        assert_quantity_allclose(s_z.phi, s.phi, atol=1e-10*u.rad)
        assert_quantity_allclose(s_z.z, s.z + 1.*u.pc)
        s_z2 = s + 1. * u.pc * sf['z'] * e['z']
        assert_representation_allclose(s_z2, s_z) 

def test_physical_spherical(self):

        s = PhysicsSphericalRepresentation(phi=[0., 6., 21.] * u.hourangle,
                                           theta=[90., 120., 5.] * u.deg,
                                           r=[1, 2, 3] * u.kpc)

        e = s.unit_vectors()
        self.check_unit_vectors(e)
        sf = s.scale_factors()
        self.check_scale_factors(sf, s)

        s_phi = s + s.r * 1e-5 * np.sin(s.theta) * e['phi']
        assert_quantity_allclose(s_phi.phi, s.phi + 1e-5*u.rad,
                                 atol=1e-10*u.rad)
        assert_quantity_allclose(s_phi.theta, s.theta, atol=1e-10*u.rad)
        assert_quantity_allclose(s_phi.r, s.r)
        s_phi2 = s + 1e-5 * u.radian * sf['phi'] * e['phi']
        assert_representation_allclose(s_phi2, s_phi)

        s_theta = s + s.r * 1e-5 * e['theta']
        assert_quantity_allclose(s_theta.phi, s.phi)
        assert_quantity_allclose(s_theta.theta, s.theta + 1e-5*u.rad,
                                 atol=1e-10*u.rad)
        assert_quantity_allclose(s_theta.r, s.r)
        s_theta2 = s + 1.e-5 * u.radian * sf['theta'] * e['theta']
        assert_representation_allclose(s_theta2, s_theta)

        s_r = s + 1. * u.pc * e['r']
        assert_quantity_allclose(s_r.phi, s.phi, atol=1e-10*u.rad)
        assert_quantity_allclose(s_r.theta, s.theta, atol=1e-10*u.rad)
        assert_quantity_allclose(s_r.r, s.r + 1.*u.pc)
        s_r2 = s + 1. * u.pc * sf['r'] * e['r']
        assert_representation_allclose(s_r2, s_r) 

def test_spherical(self):
        s = SphericalRepresentation(lon=[0., 6., 21.] * u.hourangle,
                                    lat=[0., -30., 85.] * u.deg,
                                    distance=[1, 2, 3] * u.kpc)
        e = s.unit_vectors()
        self.check_unit_vectors(e)
        sf = s.scale_factors()
        self.check_scale_factors(sf, s)

        s_lon = s + s.distance * 1e-5 * np.cos(s.lat) * e['lon']
        assert_quantity_allclose(s_lon.lon, s.lon + 1e-5*u.rad,
                                 atol=1e-10*u.rad)
        assert_quantity_allclose(s_lon.lat, s.lat, atol=1e-10*u.rad)
        assert_quantity_allclose(s_lon.distance, s.distance)
        s_lon2 = s + 1e-5 * u.radian * sf['lon'] * e['lon']
        assert_representation_allclose(s_lon2, s_lon)

        s_lat = s + s.distance * 1e-5 * e['lat']
        assert_quantity_allclose(s_lat.lon, s.lon)
        assert_quantity_allclose(s_lat.lat, s.lat + 1e-5*u.rad,
                                 atol=1e-10*u.rad)
        assert_quantity_allclose(s_lon.distance, s.distance)
        s_lat2 = s + 1.e-5 * u.radian * sf['lat'] * e['lat']
        assert_representation_allclose(s_lat2, s_lat)

        s_distance = s + 1. * u.pc * e['distance']
        assert_quantity_allclose(s_distance.lon, s.lon, atol=1e-10*u.rad)
        assert_quantity_allclose(s_distance.lat, s.lat, atol=1e-10*u.rad)
        assert_quantity_allclose(s_distance.distance, s.distance + 1.*u.pc)
        s_distance2 = s + 1. * u.pc * sf['distance'] * e['distance']
        assert_representation_allclose(s_distance2, s_distance) 

def test_sphere_cart():
    """
    Tests the spherical <-> cartesian transform functions
    """
    from astropy.utils import NumpyRNGContext
    from astropy.coordinates import spherical_to_cartesian, cartesian_to_spherical

    x, y, z = spherical_to_cartesian(1, 0, 0)
    assert_allclose(x, 1)
    assert_allclose(y, 0)
    assert_allclose(z, 0)

    x, y, z = spherical_to_cartesian(0, 1, 1)
    assert_allclose(x, 0)
    assert_allclose(y, 0)
    assert_allclose(z, 0)

    x, y, z = spherical_to_cartesian(5, 0, np.arcsin(4. / 5.))
    assert_allclose(x, 3)
    assert_allclose(y, 4)
    assert_allclose(z, 0)

    r, lat, lon = cartesian_to_spherical(0, 1, 0)
    assert_allclose(r, 1)
    assert_allclose(lat, 0 * u.deg)
    assert_allclose(lon, np.pi / 2 * u.rad)

    # test round-tripping
    with NumpyRNGContext(13579):
        x, y, z = np.random.randn(3, 5)

    r, lat, lon = cartesian_to_spherical(x, y, z)
    x2, y2, z2 = spherical_to_cartesian(r, lat, lon)

    assert_allclose(x, x2)
    assert_allclose(y, y2)
    assert_allclose(z, z2)