Return the density of the PACl solution.

    From Stock Tank Mixing report Fall 2013:
    https://confluence.cornell.edu/download/attachments/137953883/20131213_Research_Report.pdf
350859
def dens_pacl_solution(ConcAluminum, temp):
    
    return ((0.492 * ConcAluminum * PACl.MolecWeight
             / (PACl.AluminumMPM * MOLEC_WEIGHT_ALUMINUM)
             ) + pc.density_water(temp).magnitude
            )
787346
def density_angle(rho0: Density, rho1: Density) -> bk.BKTensor:
    """"""
    return fubini_study_angle(rho0.vec, rho1.vec)
312585
def ent_tank_a(self):
        
        # first guess planview area
        a_new = 1 * u.m**2
        a_ratio = 2  # set to >1+tolerance to start while loop
        tolerance = 0.01
        a_floc_pv = (
            self.floc.vol /
            (self.floc.downstream_H + (self.floc.HL / 2))
        )
        while a_ratio > (1 + tolerance):
            a_et_pv = a_new
            a_etf_pv = a_et_pv + a_floc_pv
            w_tot = a_etf_pv / self.floc.max_L
            w_chan = w_tot / self.floc.channel_n

            a_new = self.floc.max_L * w_chan
            a_ratio = a_new / a_et_pv
        return a_new
1354427
def _saturation(self, T):
        """"""
        rhoc = self._constants.get("rhoref", self.rhoc)
        Tc = self._constants.get("Tref", self.Tc)

        if T > Tc:
            T = Tc
        tau = Tc/T

        rhoLo = self._Liquid_Density(T)
        rhoGo = self._Vapor_Density(T)

        def f(parr):
            rhol, rhog = parr
            deltaL = rhol/rhoc
            deltaG = rhog/rhoc
            phirL = _phir(tau, deltaL, self._constants)
            phirG = _phir(tau, deltaG, self._constants)
            phirdL = _phird(tau, deltaL, self._constants)
            phirdG = _phird(tau, deltaG, self._constants)
            Jl = deltaL*(1+deltaL*phirdL)
            Jv = deltaG*(1+deltaG*phirdG)
            Kl = deltaL*phirdL+phirL+log(deltaL)
            Kv = deltaG*phirdG+phirG+log(deltaG)
            return Kv-Kl, Jv-Jl

        rhoL, rhoG = fsolve(f, [rhoLo, rhoGo])
        if rhoL == rhoG:
            Ps = self.Pc
        else:
            deltaL = rhoL/self.rhoc
            deltaG = rhoG/self.rhoc
            firL = _phir(tau, deltaL, self._constants)
            firG = _phir(tau, deltaG, self._constants)

            Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(firL-firG+log(deltaL/deltaG))
        return rhoL, rhoG, Ps
723993
def density_2d(self, x, y, rho0, gamma, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        sigma = np.sqrt(np.pi) * special.gamma(1./2*(-1+gamma))/special.gamma(gamma/2.) * r**(1-gamma) * rho0
        return sigma
653511
def kpoints(self):
        
        return Kpoints.automatic_density_by_vol(self.structure,
                                                self.reciprocal_density,
                                                force_gamma=True)
362481
def __equalize_densities(self,nominal_bounds,nominal_density):
        
        left,bottom,right,top = nominal_bounds.lbrt()
        width = right-left; height = top-bottom
        center_y = bottom + height/2.0
        # True density is not equal to the nominal_density when
        # nominal_density*(right-left) is not an integer.
        true_density = int(nominal_density*(width))/float(width)

        n_cells = round(height*true_density,0)
        adjusted_half_height = n_cells/true_density/2.0

        return (BoundingBox(points=((left, center_y-adjusted_half_height),
                                    (right, center_y+adjusted_half_height))),
                true_density)
761786
def sep_dist_floc(ConcAluminum, ConcClay, coag, material,
                  DIM_FRACTAL, DiamTarget):
    """"""
    return (material.Diameter
            * (np.pi/(6
                      * frac_vol_floc_initial(ConcAluminum, ConcClay,
                                              coag, material)
                      ))**(1/3)
            * (DiamTarget / material.Diameter)**(DIM_FRACTAL / 3)
            )
579223
def Q_stock(self):
        
        return self._Q_sys * (self._C_sys / self._C_stock).to(u.dimensionless)
1321067
def drain_D(self):
        
        tank_A = 2 * self.channel_L * self.channel_W
        drain_D = (np.sqrt(8 * tank_A / (np.pi * self.drain_t) * np.sqrt(
           self.downstream_H * self.drain_K / (2 * u.standard_gravity)))).to_base_units()
        return drain_D
1073091
def compute_density(start, end, length, time_unit='us'):
    
    if isinstance(start, int): start = float(start)
    if isinstance(end, int): end = float(end)
    diff = end-start
    if isinstance(diff, timedelta_types):
        if isinstance(diff, np.timedelta64):
            diff = np.timedelta64(diff, time_unit).tolist()
        tscale = 1./np.timedelta64(1, time_unit).tolist().total_seconds()
        return (length/(diff.total_seconds()*tscale))
    else:
        return length/diff
115440
def sd(self):
        ''''''
        v = self.var()
        if len(v):
            return np.sqrt(v)
        else:
            return None
1040833
def _get_rho(self, v):
        
        v_mol = vol_uc2mol(v, self.z)  # in m^3
        rho = self.mass / v_mol * 1.e-6  # in g/cm^3
        return rho
788953
def _get_charge_distribution_df(self):
        
        # Fraction coordinates and corresponding indices
        axis_grid = np.array([np.array(self.chgcar.get_axis_grid(i)) /
                              self.structure.lattice.abc[i] for i in range(3)])
        axis_index = np.array([range(len(axis_grid[i])) for i in range(3)])

        data = {}

        for index in itertools.product(*axis_index):
            a, b, c = index
            f_coords = (axis_grid[0][a], axis_grid[1][b], axis_grid[2][c])
            data[f_coords] = self.chgcar.data["total"][a][b][c]

        # Fraction coordinates - charge density table
        df = pd.Series(data).reset_index()
        df.columns = ['a', 'b', 'c', 'Charge Density']
        self._charge_distribution_df = df

        return df
671310
def dense_output(t_current, t_old, h_current, rcont):
    
    # initialization
    s = (t_current - t_old) / h_current
    s1 = 1.0 - s

    return rcont[0] + s * (rcont[1] + s1 * (
            rcont[2] + s * (rcont[3] + s1 * (rcont[4] + s * (rcont[5] + s1 * (rcont[6] + s * rcont[7]))))))
802404
def _density_par_approx(self,dangle,tdisrupt,_return_array=False,
                            higherorder=False):
        """"""
        # First construct the breakpoints for this dangle
        Oparb= (dangle-self._kick_interpdOpar_poly.x)/self._timpact
        # Find the lower limit of the integration in the pw-linear-kick approx.
        lowbindx,lowx= self.minOpar(dangle,tdisrupt,_return_raw=True)
        lowbindx= numpy.arange(len(Oparb)-1)[lowbindx]
        Oparb[lowbindx+1]= Oparb[lowbindx]-lowx
        # Now integrate between breakpoints
        out= (0.5/(1.+self._kick_interpdOpar_poly.c[-2]*self._timpact)\
                           *(special.erf(1./numpy.sqrt(2.*self._sortedSigOEig[2])\
                                             *(Oparb[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO))\
                                 -special.erf(1./numpy.sqrt(2.*self._sortedSigOEig[2])*(numpy.roll(Oparb,-1)[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO\
                                                                                             -self._kick_interpdOpar_poly.c[-2]*self._timpact*(Oparb-numpy.roll(Oparb,-1))[:-1]))))
        if _return_array:
            return out
        out= numpy.sum(out[:lowbindx+1])
        if higherorder:
            # Add higher-order contribution
            out+= self._density_par_approx_higherorder(Oparb,lowbindx)
        # Add integration to infinity
        out+= 0.5*(1.+special.erf((self._meandO-Oparb[0])\
                                  /numpy.sqrt(2.*self._sortedSigOEig[2])))
        return out
630648
def density(self):
        
        r = self.radius * _Rsun
        m = self.mass * _Msun
        return 0.75 * m / (np.pi * r * r * r)