convert density in g/cm^3 to unit-cell volume

        :param rho: density in g/cm^3
        :return: unit-cell volume in A^3

        :note: internal function
496286
def _get_v(self, rho):
        
        v_mol = self.mass / rho  # in cm^3
        v = vol_mol2uc(v_mol * 1.e-6, self.z)
        print(v)
        return v
395334
def _vmomentdensity(self,R,z,n,m,o,nsigma=None,mc=False,nmc=10000,
                       _returnmc=False,_vrs=None,_vts=None,_vzs=None,
                       _rawgausssamples=False,
                       gl=False,ngl=_DEFAULTNGL,_returngl=False,_glqeval=None,
                       _return_actions=False,_jr=None,_lz=None,_jz=None,
                       _return_freqs=False,
                       _rg=None,_kappa=None,_nu=None,_Omega=None,
                       _sigmaR1=None,_sigmaz1=None,
                       **kwargs):
        """"""
        if isinstance(R,numpy.ndarray):
            return numpy.array([self._vmomentdensity(r,zz,n,m,o,nsigma=nsigma,
                                                    mc=mc,nmc=nmc,
                                                    gl=gl,ngl=ngl,**kwargs) for r,zz in zip(R,z)])
        if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                actionAngle.actionAngleAdiabaticGrid)):
            if n % 2 == 1. or o % 2 == 1.:
                return 0. #we know this must be the case
        if nsigma == None:
            nsigma= _NSIGMA
        if _sigmaR1 is None:
            sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr)
        else:
            sigmaR1= _sigmaR1
        if _sigmaz1 is None:
            sigmaz1= self._sz*numpy.exp((self._refr-R)/self._hsz)
        else:
            sigmaz1= _sigmaz1
        thisvc= potential.vcirc(self._pot,R,use_physical=False)
        #Use the asymmetric drift equation to estimate va
        gamma= numpy.sqrt(0.5)
        va= sigmaR1**2./2./thisvc\
            *(gamma**2.-1. #Assume close to flat rotation curve, sigphi2/sigR2 =~ 0.5
               +R*(1./self._hr+2./self._hsr))
        if math.fabs(va) > sigmaR1: va = 0.#To avoid craziness near the center
        if gl:
            if ngl % 2 == 1:
                raise ValueError("ngl must be even")
            if not _glqeval is None and ngl != _glqeval.shape[0]:
                _glqeval= None
            #Use Gauss-Legendre integration for all
            if ngl == _DEFAULTNGL:
                glx, glw= self._glxdef, self._glwdef
                glx12, glw12= self._glxdef12, self._glwdef12
            elif ngl == _DEFAULTNGL2:
                glx, glw= self._glxdef2, self._glwdef2
                glx12, glw12= self._glxdef, self._glwdef
            else:
                glx, glw= numpy.polynomial.legendre.leggauss(ngl)
                glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2)
            #Evaluate everywhere
            if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                    actionAngle.actionAngleAdiabaticGrid)):
                vRgl= nsigma*sigmaR1/2.*(glx+1.)
                vzgl= nsigma*sigmaz1/2.*(glx+1.)
                vRglw= glw
                vzglw= glw
            else:
                vRgl= nsigma*sigmaR1/2.*(glx12+1.)
                #vRgl= 1.5/2.*(glx12+1.)
                vRgl= list(vRgl)
                vRgl.extend(-nsigma*sigmaR1/2.*(glx12+1.))
                #vRgl.extend(-1.5/2.*(glx12+1.))
                vRgl= numpy.array(vRgl)
                vzgl= nsigma*sigmaz1/2.*(glx12+1.)
                #vzgl= 1.5/2.*(glx12+1.)
                vzgl= list(vzgl)
                vzgl.extend(-nsigma*sigmaz1/2.*(glx12+1.))
                #vzgl.extend(-1.5/2.*(glx12+1.))
                vzgl= numpy.array(vzgl)
                vRglw= glw12
                vRglw= list(vRglw)
                vRglw.extend(glw12)
                vRglw= numpy.array(vRglw)
                vzglw= glw12
                vzglw= list(vzglw)
                vzglw.extend(glw12)
                vzglw= numpy.array(vzglw)
            if 'vTmax' in kwargs: vTmax = kwargs['vTmax']
            else:                 vTmax = 1.5
            vTgl= vTmax/2.*(glx+1.)
            #Tile everything
            vTgl= numpy.tile(vTgl,(ngl,ngl,1)).T
            vRgl= numpy.tile(numpy.reshape(vRgl,(1,ngl)).T,(ngl,1,ngl))
            vzgl= numpy.tile(vzgl,(ngl,ngl,1))
            vTglw= numpy.tile(glw,(ngl,ngl,1)).T #also tile weights
            vRglw= numpy.tile(numpy.reshape(vRglw,(1,ngl)).T,(ngl,1,ngl))
            vzglw= numpy.tile(vzglw,(ngl,ngl,1))
            #evaluate
            if _glqeval is None and _jr is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R+numpy.zeros(ngl*ngl*ngl),
                                           vRgl.flatten(),
                                           vTgl.flatten(),
                                           z+numpy.zeros(ngl*ngl*ngl),
                                           vzgl.flatten(),
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            elif not _jr is None and _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                                                 log=True,
                                                                 _return_actions=True,
                                                                 _return_freqs=True,
                                                                 use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            elif not _jr is None and not _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                           rg=_rg,kappa=_kappa,nu=_nu,
                                           Omega=_Omega,
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(ngl,ngl,ngl))
            else:
                logqeval= _glqeval
            if _returngl:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        logqeval)
            elif _return_actions and _return_freqs:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        jr,lz,jz,
                        rg,kappa,nu,Omega)
            elif _return_actions:
                return (numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                  *vTglw*vRglw*vzglw)*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2,
                        jr,lz,jz)
            else:
                return numpy.sum(numpy.exp(logqeval)*vRgl**n*vTgl**m*vzgl**o
                                 *vTglw*vRglw*vzglw*sigmaR1*sigmaz1*0.125*vTmax*nsigma**2)
        elif mc:
            mvT= (thisvc-va)/gamma/sigmaR1
            if _vrs is None:
                vrs= numpy.random.normal(size=nmc)
            else:
                vrs= _vrs
            if _vts is None:
                vts= numpy.random.normal(size=nmc)+mvT
            else:
                if _rawgausssamples:
                    vts= _vts+mvT
                else:
                    vts= _vts
            if _vzs is None:
                vzs= numpy.random.normal(size=nmc)
            else:
                vzs= _vzs
            Is= _vmomentsurfaceMCIntegrand(vzs,vrs,vts,numpy.ones(nmc)*R,
                                           numpy.ones(nmc)*z,
                                           self,sigmaR1,gamma,sigmaz1,mvT,
                                           n,m,o)
            if _returnmc:
                if _rawgausssamples:
                    return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o),
                        vrs,vts-mvT,vzs)
                else:
                    return (numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o),
                            vrs,vts,vzs)
            else:
                return numpy.mean(Is)*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
        else: #pragma: no cover because this is too slow; a warning is shown
            warnings.warn("Calculations using direct numerical integration using tplquad is not recommended and extremely slow; it has also not been carefully tested",galpyWarning)
            return integrate.tplquad(_vmomentsurfaceIntegrand,
                                     1./gamma*(thisvc-va)/sigmaR1-nsigma,
                                     1./gamma*(thisvc-va)/sigmaR1+nsigma,
                                     lambda x: 0., lambda x: nsigma,
                                     lambda x,y: 0., lambda x,y: nsigma,
                                     (R,z,self,sigmaR1,gamma,sigmaz1,n,m,o),
                                     **kwargs)[0]*sigmaR1**(2.+n+m)*gamma**(1.+m)*sigmaz1**(1.+o)
1207919
def mass_3d(self, R, Rs, rho0, r_core):
        

        Rs = float(Rs)
        b = r_core * Rs ** -1
        c = R * Rs ** -1

        M0 = 4*np.pi*Rs**3 * rho0

        return M0 * (1+b**2) ** -1 * (0.5*np.log(1+c**2) + b**2*np.log(c*b**-1 + 1) - b*np.arctan(c))
828954
def mass_3d(self, R, Rs, rho0, r_trunc):
        

        x = R * Rs ** -1

        func = (r_trunc ** 2 * (-2 * x * (1 + r_trunc ** 2) + 4 * (1 + x) * r_trunc * np.arctan(x / r_trunc) -
                                2 * (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs) + 2 * (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs * (1 + x)) +
                                2 * (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs * r_trunc) -
                                (1 + x) * (-1 + r_trunc ** 2) * np.log(Rs ** 2 * (x ** 2 + r_trunc ** 2)))) / (2. * (1 + x) * (1 + r_trunc ** 2) ** 2)

        m_3d = 4*np.pi*Rs ** 3 * rho0 * func
        return m_3d
1269933
def _convert_volume_from(self, volume_from):
        
        if ':' in volume_from:
            container, permissions = volume_from.split(':')
        else:
            container = volume_from
            permissions = 'rw'

        if permissions not in ('ro', 'rw'):
            raise ValueError("only permissions supported for volumes_from are rw and ro.")

        return "{0}:{1}".format(container, permissions)
117601
def mass_2d(self,R,Rs,rho0,r_trunc):

        

        x = R / Rs
        tau = r_trunc / Rs
        gx = self._g(x,tau)
        m_2d = 4 * rho0 * Rs * R ** 2 * gx / x ** 2 * np.pi
        return m_2d
1085691
def _get_rho(self, w, arr):
        
        lsb = 0
        while not (w & arr[lsb]):
            lsb += 1
        return lsb + 1
1207013
def _rho02alpha(self, rho0, Rs):

        

        theta_Rs = rho0 * (4 * Rs ** 2 * (1 + np.log(1. / 2.)))
        return theta_Rs
1104974
def _compute_term3(self, C, rrup):
        
        f3 = np.zeros_like(rrup)

        idx_between_70_130 = (rrup > 70) & (rrup <= 130)
        idx_greater_130 = rrup > 130

        f3[idx_between_70_130] = (
            C['c9'] * (np.log(rrup[idx_between_70_130]) - np.log(70))
        )

        f3[idx_greater_130] = (
            C['c9'] * (np.log(rrup[idx_greater_130]) - np.log(70)) +
            C['c10'] * (np.log(rrup[idx_greater_130]) - np.log(130))
        )

        return f3
689442
def mass_2d(self, R, Rs, rho0, r_core):

        

        x = R * Rs ** -1
        p = Rs * r_core ** -1
        gx = self._G(x, p)

        m_2d = 2 * np.pi * rho0 * Rs ** 3 * gx

        return m_2d
1301097
def lattice_array_to_unit_cell(lattice_array):
    """"""
    cell_lengths = np.sqrt(np.sum(lattice_array**2, axis=0))
    gamma_r = np.arccos(lattice_array[0][1] / cell_lengths[1])
    beta_r = np.arccos(lattice_array[0][2] / cell_lengths[2])
    alpha_r = np.arccos(
        lattice_array[1][2] * np.sin(gamma_r) / cell_lengths[2]
        + np.cos(beta_r) * np.cos(gamma_r)
        )
    cell_angles = [
        np.rad2deg(alpha_r), np.rad2deg(beta_r), np.rad2deg(gamma_r)
        ]
    return np.append(cell_lengths, cell_angles)
553588
def average(self, rho):
        """
        def marginal(f, rho):
            remaining = len(f.shape)
            if remaining == 0:
                return rho
            rho = sum([f[i]*rho[i] for i in range(rho.shape[0])])
            f = np.sum(f, 0)
            return marginal(f, rho)

        return marginal(self.distribution, rho)
852224
def rho2theta(self, rho0, gamma):
        
        fac = np.sqrt(np.pi) * special.gamma(1. / 2 * (-1 + gamma)) / special.gamma(gamma / 2.) * 2 / (3 - gamma) * rho0

        #fac = theta_E**(gamma - 1)
        theta_E = fac**(1. / (gamma - 1))
        return theta_E
426831
def density_water(temp):
    
    ut.check_range([temp, ">0", "Temperature in Kelvin"])
    rhointerpolated = interpolate.CubicSpline(WATER_DENSITY_TABLE[0],
                                                    WATER_DENSITY_TABLE[1])
    return np.asscalar(rhointerpolated(temp))
207657
def from_hdf5(cls, f):
        
        if isinstance(f, str):
            import h5py
            f = h5py.File(f)

        pos = quantity_from_hdf5(f['pos'])
        vel = quantity_from_hdf5(f['vel'])

        frame = None
        if 'frame' in f:
            g = f['frame']

            frame_mod = g.attrs['module']
            frame_cls = g.attrs['class']
            frame_units = [u.Unit(x.decode('utf-8')) for x in g['units']]

            if u.dimensionless_unscaled in frame_units:
                units = DimensionlessUnitSystem()
            else:
                units = UnitSystem(*frame_units)

            pars = dict()
            for k in g['parameters']:
                pars[k] = quantity_from_hdf5(g['parameters/'+k])

            exec("from {0} import {1}".format(frame_mod, frame_cls))
            frame_cls = eval(frame_cls)

            frame = frame_cls(units=units, **pars)

        return cls(pos=pos, vel=vel, frame=frame)
160338
def rho_crit(self):
        
        h = self.cosmo.H(0).value / 100.
        return 3 * h ** 2 / (8 * np.pi * const.G) * 10 ** 10 * const.Mpc / const.M_sun
1225017
def densities_close(rho0: Density, rho1: Density,
                    tolerance: float = TOLERANCE) -> bool:
    
    return vectors_close(rho0.vec, rho1.vec, tolerance)