escape velocity.
    
    Parameters
    ----------
    M : float
        Mass in solar masses.
    R : float
        Radius in solar radiu.

    Returns
    -------
    v_escape
        in km/s.
18567
def escape_velocity(M,R):
    

    ve = np.sqrt(2.*grav_const*M*msun_g/(R*rsun_cm))
    ve = ve*1.e-5
    return ve
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
696295
def mean_velocity(self, tol=0.05, Lx=None):
        '''
        
        '''
        dx = None
        dt = None
        p = None
        ind = None
        t_end = None

        if self.area == 0:
            return dict(dx=dx, dt=dt, p=p, ind=ind, t_end=t_end)

        x = self.x_position(Lx=Lx)

        # find the first and last valid indices
        ind_start = mlab.find(x.mask==False)[0]
        ind_last = mlab.find(x.mask==False)[-1]

        # if the original x value is within tol % of the final x value, include
        # all samples
        if x[ind_start] > (1 - tol) * x[ind_last] or x[ind_last] < 0:
            ind_stop = ind_last
        else: # otherwise, stop when x reaches (1 - tol) % of it's final value
            ind_stop = mlab.find(x > (1 - tol) * x[ind_last])[0]

        ind = [ind_start, ind_stop]

        # if we have at least 2 valid samples
        if len(ind) >=2:
            dx = np.diff(x[ind])[0]
            dt = np.diff(self.time[ind])[0] # ms

            # suppress polyfit warnings
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                # fit a line to the data
                p = np.polyfit(self.time[ind[0]:ind[1]], x[ind[0]:ind[1]], 1)

            # find time when the the line intercepts x[ind_last]
            ind_stop = mlab.find(self.time > \
                                 (x[ind_last] - p[1]) / p[0])
            if len(ind_stop):
                t_end = self.time[ind_stop[0]]
            else:
                t_end = self.time[-1]
        return dict(dx=dx, dt=dt, p=p, ind=ind, t_end=t_end)
866214
def mass_in_theta_E(self, theta_E):
        
        mass = self.arcsec2phys_lens(theta_E) ** 2 * np.pi * self.epsilon_crit
        return mass
393658
def c_mvar(self):
        
        if self.b == -4:
            return 1.0
        elif self.b == -3:
            return 0.0
        elif self.b == -2:
            return -1.0
        elif self.b == -1:
            return -2.0
        elif self.b == 0:
            return -3.0
92819
def mass_2d_lens(self, r, sigma0, Rs):
        
        X = r/Rs
        alpha_r = 2*sigma0 * Rs * X * (1-self._F(X)) / (X**2-1)
        mass_2d = alpha_r * r * np.pi
        return mass_2d
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
495516
def headloss_exp(FlowRate, Diam, KMinor):
    
    #Checking input validity
    ut.check_range([FlowRate, ">0", "Flow rate"], [Diam, ">0", "Diameter"],
                   [KMinor, ">=0", "K minor"])
    return KMinor * 8 / (gravity.magnitude * np.pi**2) * FlowRate**2 / Diam**4
132468
def trans_v(self, structure):
        
        nsites = structure.num_sites
        volume = structure.volume
        natoms = structure.composition.num_atoms
        weight = float(structure.composition.weight)
        mass_density = 1.6605e3 * nsites * weight / (natoms * volume)
        if self.g_vrh < 0:
            raise ValueError("k_vrh or g_vrh is negative, "
                             "sound velocity is undefined")
        return (1e9 * self.g_vrh / mass_density) ** 0.5
284425
def _mdens(self,m):
        """"""
        return (self.a/m)**self.alpha/(1.+m/self.a)**(self.betaminusalpha)
1119665
def density(self, *args):
        
        M = self.mass(*args) * MSUN
        V = 4./3 * np.pi * (self.radius(*args) * RSUN)**3
        return  M/V
735116
def period(A,M1,M2):
    

    A *= rsun_cm
    print(A)
    velocity = np.sqrt(grav_const*msun_g*(M1+M2)/A)
    print(old_div(velocity,1.e5))
    
    p = 2.*np.pi * A / velocity

    p  /= (60*60*24.)
    return p
593695
def gravity(RAW_IMU, SENSOR_OFFSETS=None, ofs=None, mul=None, smooth=0.7):
    ''''''
    if hasattr(RAW_IMU, 'xacc'):
        rx = RAW_IMU.xacc * 9.81 / 1000.0
        ry = RAW_IMU.yacc * 9.81 / 1000.0
        rz = RAW_IMU.zacc * 9.81 / 1000.0
    else:
        rx = RAW_IMU.AccX
        ry = RAW_IMU.AccY
        rz = RAW_IMU.AccZ
    if SENSOR_OFFSETS is not None and ofs is not None:
        rx += SENSOR_OFFSETS.accel_cal_x
        ry += SENSOR_OFFSETS.accel_cal_y
        rz += SENSOR_OFFSETS.accel_cal_z
        rx -= ofs[0]
        ry -= ofs[1]
        rz -= ofs[2]
        if mul is not None:
            rx *= mul[0]
            ry *= mul[1]
            rz *= mul[2]
    return sqrt(rx**2+ry**2+rz**2)
1335445
def predmean(self, xp): 
        
        out = np.empty_like(xp)
        for i in range(3):
            out[:, i] = xp[:, i] + self.delta * xp[:, i+3]  # position
            out[:, i+3] = xp[:, i+3]                      # velocity 
        return out
230692
def mass_3d(self, R, Rs, rho0, r_core):
        
        b = r_core * Rs ** -1
        x = R * Rs ** -1

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

        return M_0 * (x * (1+x) ** -1 * (-1+b) ** -1 + (-1+b) ** -2 *
                      ((2*b-1)*np.log(1/(1+x)) + b **2 * np.log(x / b + 1)))
729968
def _mdens_deriv(self,m):
        """"""
        return -2.*self.a2**2*(self.a+2.*m)/m**3/(self.a+m)**3
606280
def MAU(self):
        ''''''
        if self._MAU is None:
            self._MAU = self.linear_system.M * self.AU
        return self._MAU