Calculates the elastic energy density due to a strain
193140
def energy_density(self, strain, convert_GPa_to_eV=True):
        
        return sum([c.energy_density(strain, convert_GPa_to_eV)
                    for c in self])
1119665
def density(self, *args):
        
        M = self.mass(*args) * MSUN
        V = 4./3 * np.pi * (self.radius(*args) * RSUN)**3
        return  M/V
867246
def get_hessian(self):
        """"""
        force_const = self.fields.get("Cartesian Force Constants")
        if force_const is None:
            return None
        N = len(self.molecule.numbers)
        result = np.zeros((3*N, 3*N), float)
        counter = 0
        for row in range(3*N):
            result[row, :row+1] = force_const[counter:counter+row+1]
            result[:row+1, row] = force_const[counter:counter+row+1]
            counter += row + 1
        return result
254440
def _get_anelastic_attenuation_term(self, C, rrup):
        
        f_atn = np.zeros(len(rrup))
        idx = rrup >= 80.0
        f_atn[idx] = (C["c20"] + C["Dc20"]) * (rrup[idx] - 80.0)
        return f_atn
685693
def _update(self):
        """"""

        if self.strains.size and self.strains.size == self.values.size:
            x = np.log(self.strains)
            y = self.values

            if x.size < 4:
                self._interpolater = interp1d(
                    x,
                    y,
                    'linear',
                    bounds_error=False,
                    fill_value=(y[0], y[-1]))
            else:
                self._interpolater = interp1d(
                    x,
                    y,
                    'cubic',
                    bounds_error=False,
                    fill_value=(y[0], y[-1]))
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
482596
def unit_tangent(self, t):
        """"""
        dseg = self.derivative(t)
        return dseg/abs(dseg)
325985
def containment_angle(self, energies=None, fraction=0.68, scale_fn=None):
        """"""

        if energies is None:
            energies = self.energies

        vals = self.interp(energies[np.newaxis, :], self.dtheta[:, np.newaxis],
                           scale_fn=scale_fn)
        dtheta = np.radians(self.dtheta[:, np.newaxis] * np.ones(vals.shape))
        return self._calc_containment(dtheta, vals, fraction)
743844
def density(self):
        
        r = float(self.sp_index.npoints) / float(self.sp_index.length)
        return r
705375
def spectral_rolloff(data, coeff):
    """"""
    data = np.mean(data, axis=1)
    nFFT = len(data) // 2
    X = FFT(data, nFFT)

    totalEnergy = np.sum(X ** 2)
    fftLength = len(X)
    Thres = coeff * totalEnergy

    # Find the spectral rolloff as the frequency position where the
    # respective spectral energy is equal to c*totalEnergy
    CumSum = np.cumsum(X ** 2) + EPSILON
    [a, ] = np.nonzero(CumSum > Thres)
    if len(a) > 0:
        mC = np.float64(a[0]) / (float(fftLength))
    else:
        mC = 0.0

    return mC
687273
def get_energy(self, circuit, sampler):
        """"""
        val = 0.0
        for meas in self.hamiltonian:
            c = circuit.copy()
            for op in meas.ops:
                if op.op == "X":
                    c.h[op.n]
                elif op.op == "Y":
                    c.rx(-np.pi / 2)[op.n]
            measured = sampler(c, meas.n_iter())
            for bits, prob in measured.items():
                if sum(bits) % 2:
                    val -= prob * meas.coeff
                else:
                    val += prob * meas.coeff
        return val.real
276058
def get_forces(self, a):
        """"""
        f = np.zeros( [ len(a), 3 ], dtype=float )
        for c in self.calcs:
            f += c.get_forces(a)
        return f
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
389129
def stSpectralRollOff(X, c, fs):
    """"""
    totalEnergy = numpy.sum(X ** 2)
    fftLength = len(X)
    Thres = c*totalEnergy
    # Ffind the spectral rolloff as the frequency position 
    # where the respective spectral energy is equal to c*totalEnergy
    CumSum = numpy.cumsum(X ** 2) + eps
    [a, ] = numpy.nonzero(CumSum > Thres)
    if len(a) > 0:
        mC = numpy.float64(a[0]) / (float(fftLength))
    else:
        mC = 0.0
    return (mC)
973910
def make_force_ES(q_vars, mass):
    """"""
    # Build the potential energy fluxion; just one pair of bodies
    U =  U_ij(q_vars, mass, 0, 1)
    # 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
630648
def density(self):
        
        r = self.radius * _Rsun
        m = self.mass * _Msun
        return 0.75 * m / (np.pi * r * r * r)
761968
def density(self):
        
        m = Mass(self.composition.weight, "amu")
        return m.to("g") / (self.volume * Length(1, "ang").to("cm") ** 3)