r"""
        The potential energy *per unit mass*:

        .. math::

            E_\Phi = \Phi(\boldsymbol{q})

        Returns
        -------
        E : :class:`~astropy.units.Quantity`
            The potential energy.
427254
def potential_energy(self, potential=None):
        
        """
        if self.hamiltonian is None and potential is None:
            raise ValueError("To compute the potential energy, a potential"
                             " object must be provided!")
        if potential is None:
            potential = self.hamiltonian.potential

        return super(Orbit,self).potential_energy(potential)
1084182
def flux_units(self, units):
        
        # Check that the units are valid
        dtypes = (q.core.PrefixUnit, q.quantity.Quantity, q.core.CompositeUnit)
        if not isinstance(units, dtypes):
            raise ValueError(units, "units not understood.")

        # Check that the units changed
        if units != self.flux_units:

            # Convert to new units
            sfd = q.spectral_density(self.wave_eff)
            self.zp = self.zp.to(units, equivalencies=sfd)

            # Store new units
            self._flux_units = units
156971
def accel_fl(q: np.ndarray):
    """"""
    # Infer number of dimensions from q
    dims: int = len(q)
    # Number of celestial bodies
    B: int = dims // 3

    # The force given the positions q of the bodies
    f = force(q)

    # The accelerations from this force
    a = np.zeros(dims)
    for i in range(B):
        a[slices[i]] = f[slices[i]] / mass[i]
    
    return a
659109
def est_entropy(self):
        
        """
        nz_weights = self.particle_weights[self.particle_weights > 0]
        return -np.sum(np.log(nz_weights) * nz_weights)
443585
def U_ij(q_vars: List[fl.Var], mass: np.ndarray, i: int, j: int):
    """"""
    # Check that the lengths are consistent
    assert len(q_vars) == 3 * len(mass)
    
    # Masses of the bodies i and j
    mi = mass[i]
    mj = mass[j]
    
    # Gravitational potential is -G * m1 * m2 / r
    U = -(G * mi * mj) / flux_r(q_vars, i, j)
    return U
1198636
def _get_energy(self, x):
        
        return self.pd.get_hull_energy(self.comp1 * x + self.comp2 * (1-x)) - \
            self.e1 * x - self.e2 * (1-x)
895174
def ellipticity2phi_q(e1, e2):
    
    phi = np.arctan2(e2, e1)/2
    c = np.sqrt(e1**2+e2**2)
    if c > 0.999:
        c = 0.999
    q = (1-c)/(1+c)
    return phi, q
480117
def _R2deriv(self,R,phi=0.,t=0.):
        
        return self._Pot.R2deriv(R,0.,phi=phi,t=t,use_physical=False)
73971
def calculate_energy(self, energies):
        
        return sum([amt * energies[c] for amt, c in zip(self._coeffs,
                                                        self._all_comp)])
280905
function(phi, c) {
  this.phi = phi;
  this.c = c;

  this.psi = 1.0 + phi * phi / 2;
  this.zeta = 1.0 + phi * phi;

  this.mu = {};
  this.sigma = {};
  return this;
}
1149702
def get_flux(self, energies):

        """"""

        results = [component.shape(energies) for component in self.components.values()]

        return numpy.sum(results, 0)
589794
def energy(self):
        
        return self.uncorrected_energy + np.sum(list(self.corrections.values()))
297626
def _get_phi(self, C, mag):
        
        if mag < 5.5:
            return C["phi1"]
        elif mag < 5.75:
            return C["phi1"] + (C["phi2"] - C["phi1"]) * ((mag - 5.5) / 0.25)
        else:
            return C["phi2"]
401020
def get_gammadot(F, mc, q, e):
    

    # chirp mass
    mc *= SOLAR2S

    #total mass
    m = (((1+q)**2)/q)**(3/5) * mc

    dgdt = 6*np.pi*F * (2*np.pi*F*m)**(2/3) / (1-e**2) * \
        (1 + 0.25*(2*np.pi*F*m)**(2/3)/(1-e**2)*(26-15*e**2))

    return dgdt
17384
def unit(self):
        """"""
        s = self.x * self.x + self.y * self.y
        if s < 1e-5:
            return Point(0,0)
        s = math.sqrt(s)
        return Point(self.x / s, self.y / s)
761968
def density(self):
        
        m = Mass(self.composition.weight, "amu")
        return m.to("g") / (self.volume * Length(1, "ang").to("cm") ** 3)
1222532
def make_force(q_vars, mass):
    """"""
    # Number of bodies
    B: int = len(mass)
    
    # Build the potential energy fluxion by iterating over distinct pairs of bodies
    U = fl.Const(0.0)
    for i in range(B):
        for j in range(i+1, B):
            U +=  U_ij(q_vars, mass, i, j)
    
    # Varname arrays for both the coordinate system and U
    vn_q = np.array([q.var_name for q in q_vars])
    vn_fl = np.array(sorted(U.var_names))
    # Permutation array for putting variables in q in the order expected by U (alphabetical)
    q2fl = np.array([np.argmax((vn_q == v)) for v in vn_fl])
    # Permutation array for putting results of U.diff() in order of q_vars
    fl2q = np.array([np.argmax((vn_fl == v)) for v in vn_q])
    # Return a force function from this potential
    force_func = lambda q: -U.diff(q[q2fl]).squeeze()[fl2q]
    return force_func