Python G.decompose方法代码示例

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def mass_enclosed(self, x):
        """
        Estimate the mass enclosed within the given position by assumine the potential
        is spherical. This is basic, and assumes spherical symmetry.

        Parameters
        ----------
        x : array_like, numeric
            Position to estimate the enclossed mass.
        """

        # Fractional step-size in radius
        h = 0.01

        # Radius
        r = np.sqrt(np.sum(x**2, axis=-1))

        epsilon = h*x/r[...,np.newaxis]

        dPhi_dr_plus = self.value(x + epsilon)
        dPhi_dr_minus = self.value(x - epsilon)
        diff = dPhi_dr_plus - dPhi_dr_minus

        if self.units is None:
            raise ValueError("No units specified when creating potential object.")
        Gee = G.decompose(self.units).value

        return np.abs(r*r * diff / Gee / (2.*h))

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def __init__(self, m, a, units):
        self.parameters = OrderedDict()
        self.parameters['m'] = m
        self.parameters['a'] = a

        super(KuzminPotential, self).__init__(units=units)
        self.G = G.decompose(units).value

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def _read_units(self):
        """
        Read and parse the SCFPAR file containing simulation parameters
        and initial conditions. Right now, only parse out the simulation units.
        """
        pars = dict()

        parfile = os.path.join(self.path, "SCFPAR")
        with open(parfile) as f:
            lines = f.readlines()

            # find what G is set to
            for i,line in enumerate(lines):
                if line.split()[1].strip() == "G":
                    break
            pars['G'] = float(lines[i].split()[0])
            pars['length'] = float(lines[i+10].split()[0])
            pars['mass'] = float(lines[i+11].split()[0])

            self.x0 = np.array(map(float, lines[19].split()))*u.kpc
            self.v0 = np.array(map(float, lines[20].split()))*u.km/u.s

        _G = G.decompose(bases=[u.kpc,u.M_sun,u.Myr]).value
        X = (_G / pars['length']**3 * pars['mass'])**-0.5

        length_unit = u.Unit("{0} kpc".format(pars['length']))
        mass_unit = u.Unit("{0} M_sun".format(pars['mass']))
        time_unit = u.Unit("{:08f} Myr".format(X))

        # units = dict(length=length_unit,
        #              mass=mass_unit,
        #              time=time_unit,
        #              speed=length_unit/time_unit,
        #              dimensionless=u.dimensionless_unscaled)
        return UnitSystem(length_unit, mass_unit, time_unit, length_unit/time_unit, u.radian)

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
def test_against_gary():
    cos_coeff = np.array([[[1.509, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [-2.606, 0.0, 0.665, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [6.406, 0.0, -0.66, 0.0, 0.044, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [-5.859, 0.0, 0.984, 0.0, -0.03, 0.0, 0.001]], [[-0.086, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [-0.221, 0.0, 0.129, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [1.295, 0.0, -0.14, 0.0, -0.012, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]], [[-0.033, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [-0.001, 0.0, 0.006, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]], [[-0.02, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]])
    sin_coeff = np.zeros_like(cos_coeff)
    nmax = sin_coeff.shape[0]-1
    lmax = sin_coeff.shape[1]-1

    max_x = 10.
    grid = np.linspace(-max_x,max_x,8)
    grid = grid[grid != 0.]
    xyz = np.ascontiguousarray(np.vstack(map(np.ravel, np.meshgrid(grid,grid,grid))).T)
    # xyz = np.array([[1.,1.,1.]])

    grad = np.zeros_like(xyz)
    acceleration(xyz, grad, sin_coeff, cos_coeff, nmax, lmax)
    grad *= -1.

    val = np.zeros(len(grad))
    value(xyz, val, sin_coeff, cos_coeff, nmax, lmax)

    # gary
    import gary.potential as gp
    from gary.units import galactic
    G = _G.decompose(galactic).value
    potential = gp.SCFPotential(m=1/G, r_s=1.,
                                sin_coeff=sin_coeff, cos_coeff=cos_coeff,
                                units=galactic)
    gary_val = potential.value(xyz)
    gary_grad = potential.gradient(xyz)

    np.testing.assert_allclose(gary_val, val, rtol=1E-5)
    np.testing.assert_allclose(gary_grad, grad, rtol=1E-5)

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def __init__(self, parameters, origin=None, R=None,
                 ndim=3, units=None):

        self._units = self._validate_units(units)
        self.parameters = self._prepare_parameters(parameters, self.units)

        try:
            self.G = G.decompose(self.units).value
        except u.UnitConversionError:
            self.G = 1. # TODO: this is a HACK and could lead to user confusion

        self.ndim = ndim

        if origin is None:
            origin = np.zeros(self.ndim)
        self.origin = self._remove_units(origin)

        if R is not None and ndim not in [2, 3]:
            raise NotImplementedError('Gala potentials currently only support '
                                      'rotations when ndim=2 or ndim=3.')

        if R is not None:
            if isinstance(R, Rotation):
                R = R.as_dcm()
            R = np.array(R)

            if R.shape != (ndim, ndim):
                raise ValueError('Rotation matrix passed to potential {0} has '
                                 'an invalid shape: expected {1}, got {2}'
                                 .format(self.__class__.__name__,
                                         (ndim, ndim), R.shape))
        self.R = R

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def mass_enclosed(self, q, t=0.):
        """
        Estimate the mass enclosed within the given position by assuming the potential
        is spherical.

        Parameters
        ----------
        q : `~gala.dynamics.PhaseSpacePosition`, `~astropy.units.Quantity`, array_like
            Position(s) to estimate the enclossed mass.

        Returns
        -------
        menc : `~astropy.units.Quantity`
            Mass enclosed at the given position(s). If the input position
            has shape ``q.shape``, the output energy will have shape
            ``q.shape[1:]``.
        """
        q = self._remove_units_prepare_shape(q)
        orig_shape, q = self._get_c_valid_arr(q)
        t = self._validate_prepare_time(t, q)

        # small step-size in direction of q
        h = 1E-3 # MAGIC NUMBER

        # Radius
        r = np.sqrt(np.sum(q**2, axis=1))

        epsilon = h*q/r[:, np.newaxis]

        dPhi_dr_plus = self._energy(q + epsilon, t=t)
        dPhi_dr_minus = self._energy(q - epsilon, t=t)
        diff = (dPhi_dr_plus - dPhi_dr_minus)

        if isinstance(self.units, DimensionlessUnitSystem):
            Gee = 1.
        else:
            Gee = G.decompose(self.units).value

        Menc = np.abs(r*r * diff / Gee / (2.*h))
        Menc = Menc.reshape(orig_shape[1:])

        sgn = 1.
        if 'm' in self.parameters and self.parameters['m'] < 0:
            sgn = -1.

        return sgn * Menc * self.units['mass']

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
def _units_from_file(scfpar):
    """ Generate a unit system from an SCFPAR file. """

    with open(scfpar) as f:
        lines = f.readlines()
        length = float(lines[16].split()[0])
        mass = float(lines[17].split()[0])

    GG = G.decompose(bases=[u.kpc,u.M_sun,u.Myr]).value
    X = (GG / length**3 * mass)**-0.5

    length_unit = u.Unit("{0} kpc".format(length))
    mass_unit = u.Unit("{0} M_sun".format(mass))
    time_unit = u.Unit("{:08f} Myr".format(X))

    return dict(length=length_unit,
                mass=mass_unit,
                time=time_unit)

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def __init__(self, parameters, units=None):
        # make sure the units specified are a UnitSystem instance
        if units is not None and not isinstance(units, UnitSystem):
            units = UnitSystem(*units)

        elif units is None:
            units = DimensionlessUnitSystem()

        self.units = units

        # in case the user specified an ordered dict
        self.parameters = OrderedDict()
        for k,v in parameters.items():
            if hasattr(v, 'unit'):
                self.parameters[k] = v.decompose(self.units)
            else:
                self.parameters[k] = v*u.one

        try:
            self.G = G.decompose(self.units).value
        except u.UnitConversionError:
            self.G = 1. # TODO: this is a HACK and could lead to user confusion

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def mass_enclosed(self, q, t=0.):
        """
        Estimate the mass enclosed within the given position by assuming the potential
        is spherical.

        Parameters
        ----------
        x : array_like, numeric
            Position to estimate the enclossed mass.

        Returns
        -------
        menc : `~astropy.units.Quantity`
            The potential energy or value of the potential. If the input
            position has shape ``q.shape``, the output energy will have
            shape ``q.shape[1:]``.
        """
        q = self._prefilter_pos(q)

        # Fractional step-size in radius
        h = 0.01

        # Radius
        r = np.sqrt(np.sum(q**2, axis=0))

        epsilon = h*q/r[np.newaxis]

        dPhi_dr_plus = self.value(q + epsilon, t=t)
        dPhi_dr_minus = self.value(q - epsilon, t=t)
        diff = dPhi_dr_plus - dPhi_dr_minus

        if isinstance(self.units, DimensionlessUnitSystem):
            raise ValueError("No units specified when creating potential object.")
        Gee = G.decompose(self.units).value

        return np.abs(r*r * diff / Gee / (2.*h)) * self.units['mass']

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    def mass_enclosed(self, q, t=0.):
        """
        Estimate the mass enclosed within the given position by assuming the potential
        is spherical.

        Parameters
        ----------
        x : array_like, numeric
            Position to estimate the enclossed mass.

        Returns
        -------
        menc : `~numpy.ndarray`
            The mass. Will have the same shape as the input position,
            array, ``q``, but without the coordinate axis, ``axis=0``
        """

        q = np.ascontiguousarray(atleast_2d(q, insert_axis=1))

        # Fractional step-size in radius
        h = 0.01

        # Radius
        r = np.sqrt(np.sum(q**2, axis=0))

        epsilon = h*q/r[np.newaxis]

        dPhi_dr_plus = self.value(q + epsilon, t=t)
        dPhi_dr_minus = self.value(q - epsilon, t=t)
        diff = dPhi_dr_plus - dPhi_dr_minus

        if self.units is None:
            raise ValueError("No units specified when creating potential object.")
        Gee = G.decompose(self.units).value

        return np.abs(r*r * diff / Gee / (2.*h))

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
def total_rv():

    filenamer = os.path.join(plot_path, "rel_r.png")
    filenamev = os.path.join(plot_path, "rel_v.png")

    figr,axesr = plt.subplots(4,1,figsize=(10,14),
                              sharex=True)
    figv,axesv = plt.subplots(4,1,figsize=(10,14),
                              sharex=True)

    nparticles = 2000
    for k,_m in enumerate(range(6,9+1)):
        mass = "2.5e{}".format(_m)
        m = float(mass)
        print(mass)

        sgr = SgrSimulation(sgr_path.format(_m),snapfile)
        p = sgr.particles(n=nparticles, expr=expr)
        s = sgr.satellite()

        X = np.vstack((s._X[...,:3], p._X[...,:3].copy()))
        V = np.vstack((s._X[...,3:], p._X[...,3:].copy()))
        integrator = LeapfrogIntegrator(sgr.potential._acceleration_at,
                                        np.array(X), np.array(V),
                                        args=(X.shape[0], np.zeros_like(X)))
        ts, rs, vs = integrator.run(t1=sgr.t1, t2=sgr.t2, dt=-1.)

        s_orbit = np.vstack((rs[:,0][:,np.newaxis].T, vs[:,0][:,np.newaxis].T)).T
        p_orbits = np.vstack((rs[:,1:].T, vs[:,1:].T)).T
        t_idx = np.array([np.argmin(np.fabs(ts - t)) for t in p.tub])

        m_t = (-s.mdot*ts + s.m0)[:,np.newaxis]
        s_R = np.sqrt(np.sum(s_orbit[...,:3]**2, axis=-1))
        s_V = np.sqrt(np.sum(s_orbit[...,3:]**2, axis=-1))
        r_tide = sgr.potential._tidal_radius(m_t, s_orbit[...,:3])
        v_disp = s_V * r_tide / s_R

        # cartesian basis to project into
        x_hat = s_orbit[...,:3] / np.sqrt(np.sum(s_orbit[...,:3]**2, axis=-1))[...,np.newaxis]
        _y_hat = s_orbit[...,3:] / np.sqrt(np.sum(s_orbit[...,3:]**2, axis=-1))[...,np.newaxis]
        z_hat = np.cross(x_hat, _y_hat)
        y_hat = -np.cross(x_hat, z_hat)

        # translate to satellite position
        rel_orbits = p_orbits - s_orbit
        rel_pos = rel_orbits[...,:3]
        rel_vel = rel_orbits[...,3:]

        # project onto each
        X = np.sum(rel_pos * x_hat, axis=-1)
        Y = np.sum(rel_pos * y_hat, axis=-1)
        Z = np.sum(rel_pos * z_hat, axis=-1)
        RR = np.sqrt(X**2 + Y**2 + Z**2)

        VX = np.sum(rel_vel * x_hat, axis=-1)
        VY = np.sum(rel_vel * y_hat, axis=-1)
        VZ = np.sum(rel_vel * z_hat, axis=-1)
        VV = (np.sqrt(VX**2 + VY**2 + VZ**2)*u.kpc/u.Myr).to(u.km/u.s).value
        v_disp = (v_disp*u.kpc/u.Myr).to(u.km/u.s).value

        _tcross = r_tide / np.sqrt(G.decompose(usys).value*m/r_tide)
        for ii,jj in enumerate(t_idx):
            #tcross = r_tide[jj,0] / _v[jj,ii]
            tcross = _tcross[jj]
            bnd = int(tcross / 2)

            ix1,ix2 = jj-bnd, jj+bnd
            if ix1 < 0: ix1 = 0
            if ix2 > max(sgr.t1,sgr.t2): ix2 = -1
            axesr[k].plot(ts[ix1:ix2],
                          RR[ix1:ix2,ii],
                          linestyle='-', alpha=0.1, marker=None, color='#555555', zorder=-1)

            axesv[k].plot(ts[ix1:ix2],
                          VV[ix1:ix2,ii],
                          linestyle='-', alpha=0.1, marker=None, color='#555555', zorder=-1)

        axesr[k].plot(ts, r_tide*2., marker=None)

        axesr[k].set_xlim(ts.min(), ts.max())
        axesv[k].set_xlim(ts.min(), ts.max())

        axesr[k].set_ylim(0,max(r_tide)*7)
        axesv[k].set_ylim(0,max(v_disp)*7)

        # axes[1,k].set_xlabel(r"$x_1$")

        # if k == 0:
        #     axes[0,k].set_ylabel(r"$x_2$")
        #     axes[1,k].set_ylabel(r"$x_3$")

        axesr[k].text(3000, max(r_tide)*5, r"$2.5\times10^{}M_\odot$".format(_m))
        axesv[k].text(3000, max(v_disp)*5, r"$2.5\times10^{}M_\odot$".format(_m))

    axesr[-1].set_xlabel("time [Myr]")
    axesv[-1].set_xlabel("time [Myr]")

    figr.suptitle("Relative distance", fontsize=26)
    figr.tight_layout()
    figr.subplots_adjust(top=0.92, hspace=0.025, wspace=0.1)
#.........这里部分代码省略.........

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
def test_compose_into_arbitrary_units():
    # Issue #1438
    from astropy.constants import G
    G.decompose([u.kg, u.km, u.Unit("15 s")])

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
from matplotlib import animation
import matplotlib.pyplot as plt
import numpy as np
from astropy.constants import G
import astropy.units as u
from scipy.signal import argrelmin
from streamteam.integrate import DOPRI853Integrator
from streamteam.dynamics import lyapunov

# Create logger
logger = logging.getLogger(__name__)

plot_path = "plots"

usys = [u.kpc, u.Myr, u.radian, u.M_sun]
_G = G.decompose(usys).value
# standard prolate:
prolate_params = (_G, 1.63E11, 3., 6., 0.2,
                  5.8E9, 0.25,
                  0.204542433, 0.5, 8.5)

# standard oblate:
oblate_params = (_G, 1.63E11, 3., 6., 0.2,
                 5.8E9, 0.25,
                 0.204542433, 1.5, 8.5)

def zotos_potential(R, z, G, Md, alpha, b, h, Mn, cn, v0, beta, ch):
    Rsq = R*R

    Vd = -G*Md / np.sqrt(b**2 + Rsq + (alpha + np.sqrt(h**2+z**2))**2)
    Vn = -G*Mn / np.sqrt(Rsq+z**2+cn**2)

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
    elif args.quiet:
        logger.setLevel(logging.ERROR)
    else:
        logger.setLevel(logging.INFO)

    # Contains user-specified parameters for SCF
    scfpars = dict()

    if args.rscale is None:
        ru = 0.43089*(args.mass/2.5e9)**(1/3.)
    else:
        ru = args.rscale
    scfpars['rscale'] = ru
    scfpars['mass'] = args.mass

    _G = G.decompose(bases=[u.kpc,u.M_sun,u.Myr]).value
    X = (_G / ru**3 * args.mass)**-0.5

    length_unit = u.Unit("{0} kpc".format(ru))
    mass_unit = u.Unit("{0} M_sun".format(args.mass))
    time_unit = u.Unit("{:08f} Myr".format(X))

    scfpars['dt'] = args.dt / (1*time_unit).to(u.Myr).value
    scfpars['nsteps'] = args.nsteps
    scfpars['ncen'] = args.ncen
    scfpars['nsnap'] = args.nsnap
    scfpars['ntide'] = args.ntide
    if args.scfbi_path is None:
        scfpars['SCFBIpath'] = os.path.abspath(os.path.join(project_path, 'src', 'SCFBI'))

    main(name=args.name, pos=args.x, vel=args.v, scfpars=scfpars,

# 需要导入模块: from astropy.constants import G [as 别名]
# 或者: from astropy.constants.G import decompose [as 别名]
 def __init__(self, m, a, units):
     self.parameters = dict(m=m, a=a)
     self.units = units
     self.G = G.decompose(units).value