Returns total flux of source at the given energy
        :param energies: energies (array or float)
        :return: differential flux at given energy
437081
def get_spatially_integrated_flux( self, energies):
    
        
       
        if not isinstance(energies, np.ndarray):
            energies = np.array(energies, ndmin=1)

        # Get the differential flux from the spectral components

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

        if isinstance(energies, u.Quantity):

            # Slow version with units

            # We need to sum like this (slower) because using np.sum will not preserve the units
            # (thanks astropy.units)

            differential_flux = sum(results)

        else:

            # Fast version without units, where x is supposed to be in the same units as currently defined in
            # units.get_units()

            differential_flux = np.sum(results, 0)

        return differential_flux
1146758
def get_potential_energy(self, a):
        """"""
        e = 0.0
        for c in self.calcs:
            e += c.get_potential_energy(a)
        return e
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)
1231179
def _energy_evaluation(self, operator):
        
        if self._quantum_state is not None:
            input_circuit = self._quantum_state
        else:
            input_circuit = [self.opt_circuit]
        if operator._paulis:
            mean_energy, std_energy = operator.evaluate_with_result(self._operator_mode, input_circuit,
                                                                    self._quantum_instance.backend, self.ret)
        else:
            mean_energy = 0.0
            std_energy = 0.0
        operator.disable_summarize_circuits()
        logger.debug('Energy evaluation {} returned {}'.format(self._eval_count, np.real(mean_energy)))
        return np.real(mean_energy), np.real(std_energy)
1371816
def estimate_flux(self):
        

        from .epsf import _interpolate_missing_data

        if np.any(self.mask):
            data_interp = _interpolate_missing_data(self.data, method='cubic',
                                                    mask=self.mask)
            data_interp = _interpolate_missing_data(data_interp,
                                                    method='nearest',
                                                    mask=self.mask)
            flux = np.sum(data_interp, dtype=np.float64)

        else:
            flux = np.sum(self.data, dtype=np.float64)

        return flux
873339
def total_power(self):
        
        power = self.average_current * self.voltage
        return round(power, self.sr)
284464
def get_potential_energy(self, a=None):
        """"""
        if a is None:
            a = self.a
        if self.mask is None:
            return -np.sum(np.dot(a.get_positions()[self.mask], self.force))
        else:
            return -np.sum(np.dot(a.get_positions(), self.force))
1349549
def smoothness(self):
        
        energies = self.energies
        try:
            sp = self.spline()
        except:
            print("Energy spline failed.")
            return None
        spline_energies = sp(range(len(energies)))
        diff = spline_energies - energies
        rms = np.sqrt(np.sum(np.square(diff)) / len(energies))
        return rms
1206852
def _generate_energy_edges_single(ene):
    """"""
    midene = np.sqrt((ene[1:] * ene[:-1]))
    elo, ehi = np.zeros(len(ene)) * ene.unit, np.zeros(len(ene)) * ene.unit
    elo[1:] = ene[1:] - midene
    ehi[:-1] = midene - ene[:-1]
    elo[0] = ene[0] * (1 - ene[0] / (ene[0] + ehi[0]))
    ehi[-1] = elo[-1]
    return u.Quantity([elo, ehi])
137217
def _subtract_path_flux(net_flux, path):
    

    net_flux = copy.copy(net_flux)

    net_flux[path[:-1], path[1:]] -= net_flux[path[:-1], path[1:]].min()

    # The above *should* make the bottleneck have zero flux, but
    # numerically that may not be the case, so just set it to zero
    # to be sure.
    bottleneck_ind = net_flux[path[:-1], path[1:]].argmin()
    net_flux[path[bottleneck_ind], path[bottleneck_ind + 1]] = 0.0

    return net_flux
844494
def get_form_energy(self, entry):
        
        c = entry.composition
        return entry.energy - sum([c[el] * self.el_refs[el].energy_per_atom
                                   for el in c.elements])
1347429
def energy(self):
        
        """
        s, b, W, N = self.state, self.b, self.W, self.N
        self.E = - sum(s * b) - sum([s[i] * s[j] * W[i, j] for (i, j) in product(range(N), range(N)) if i < j])
        self.low_energies[-1] = self.E
        self.low_energies.sort()
        self.high_energies[-1] = self.E
        self.high_energies.sort()
        self.high_energies = self.high_energies[::-1]
        return self.E
179920
def read_spectral_data(hdu):
    
    ebins = read_energy_bounds(hdu)
    fluxes = np.ndarray((len(ebins)))
    try:
        fluxes[0:-1] = hdu.data.field("E_MIN_FL")
        fluxes[-1] = hdu.data.field("E_MAX_FL")[-1]
        npreds = hdu.data.field("NPRED")
    except:
        fluxes = np.ones((len(ebins)))
        npreds = np.ones((len(ebins)))
    return ebins, fluxes, npreds
113204
def calculate_energy(self):
        
        clibrebound.reb_tools_energy.restype = c_double
        return clibrebound.reb_tools_energy(byref(self))
138519
def energy(self):
        """"""
        result = 0.0
        for index1 in range(self.numc):
            for index2 in range(index1):
                if self.scaling[index1, index2] > 0:
                    for se, ve in self.yield_pair_energies(index1, index2):
                        result += se*ve*self.scaling[index1, index2]
        return result
615585
def determine_coefficient
      # Array of photon energies (in units of MeV).
      @energies = [
        0.001,
        0.0015,
        0.002,
        0.003,
        0.004,
        0.005,
        0.006,
        0.008,
        0.01,
        0.015,
        0.02,
        0.03,
        0.04,
        0.05,
        0.06,
        0.08,
        0.1,
        0.15,
        0.2,
        0.3,
        0.4,
        0.5,
        0.6,
        0.8,
        1.0,
        1.25,
        1.5,
        2.0,
        3.0,
        4.0,
        5.0,
        6.0,
        8.0,
        10.0,
        15.0,
        20.0
      ]
      # Array of mass attenuation coefficients for the above energies, for liquid water (in units of cm^2/g):
      @att_coeffs = [
        4078,
        1376,
        617.3,
        192.9,
        82.78,
        42.58,
        24.64,
        10.37,
        5.329,
        1.673,
        0.8096,
        0.3756,
        0.2683,
        0.2269,
        0.2059,
        0.1837,
        0.1707,
        0.1505,
        0.137,
        0.1186,
        0.1061,
        0.09687,
        0.08956,
        0.07865,
        0.07072,
        0.06323,
        0.05754,
        0.04942,
        0.03969,
        0.03403,
        0.03031,
        0.0277,
        0.02429,
        0.02219,
        0.01941,
        0.01813
      ]
      # Determine the coefficient:
      if @energy >= 20.0
        # When the energy is above 20, we use the coefficient for 20 MeV:
        @ac_water = @att_coeffs.last
      else
        if i = @energies.index(@energy)
          # When it exactly matches one of the listed energies, use its corresponding coefficient:
          @ac_water = @att_coeffs[i]
        else
          # When the given energy is between two of the listed values, interpolate:
          i_after = @energies.index {|x| x > @energy}
          if i_after
            e_high = @energies[i_after]
            e_low = @energies[i_after - 1]
            ac_high = @att_coeffs[i_after]
            ac_low = @att_coeffs[i_after - 1]
            @ac_water = (ac_high - ac_low ) / (e_high - e_low) * (@energy - e_low)
          else
            raise "Unexpected behaviour with index in the energy interpolation logic."
          end
        end
      end
    end
1272364
def energy(mesh, uniform_density=False):
    
    # E = 1/(d+1) sum_i ||x_i||^2 |omega_i| - int_Omega_i ||x||^2
    dim = mesh.cells["nodes"].shape[1] - 1

    star_volume = numpy.zeros(mesh.node_coords.shape[0])
    for i in range(3):
        idx = mesh.cells["nodes"][:, i]
        if uniform_density:
            # rho = 1,
            # int_{star} phi_i * rho = 1/(d+1) sum_{triangles in star} |triangle|
            fastfunc.add.at(star_volume, idx, mesh.cell_volumes)
        else:
            # rho = 1 / tau_j,
            # int_{star} phi_i * rho = 1/(d+1) |num triangles in star|
            fastfunc.add.at(star_volume, idx, numpy.ones(idx.shape, dtype=float))
    x2 = numpy.einsum("ij,ij->i", mesh.node_coords, mesh.node_coords)
    out = 1 / (dim + 1) * numpy.dot(star_volume, x2)

    # could be cached
    assert dim == 2
    x = mesh.node_coords[:, :2]
    triangles = numpy.moveaxis(x[mesh.cells["nodes"]], 0, 1)
    val = quadpy.triangle.integrate(
        lambda x: x[0] ** 2 + x[1] ** 2,
        triangles,
        # Take any scheme with order 2
        quadpy.triangle.Dunavant(2),
    )
    if uniform_density:
        val = numpy.sum(val)
    else:
        rho = 1.0 / mesh.cell_volumes
        val = numpy.dot(val, rho)

    assert out >= val

    return out - val