Find the total DOS shifted by the Fermi energy
967465
def get_dos(self):
        ''''''
        # find the dos file
        fildos = ''
        for f in self._files:
            with open(f, 'r') as fp:
                first_line = next(fp)
                if "E (eV)" in first_line and "Int dos(E)" in first_line:
                    fildos = f
                    ndoscol = len(next(fp).split())-2 # number of spin channels
                    fp.close()
                    break
                fp.close()
        if not fildos:
            return None # cannot find DOS

        # get the Fermi energy
        line = self._get_line('the Fermi energy is', self.outputf)
        efermi = float(line.split('is')[-1].split()[0])

        # grab the DOS
        energy = [] ; dos = []
        fp = open(fildos, 'r')
        next(fp) # comment line
        for line in fp:
            ls = line.split()
            energy.append(Scalar(value=float(ls[0])-efermi))
            dos.append(Scalar(value=sum([float(i) for i in ls[1:1+ndoscol]])))
        return Property(scalars=dos, units='number of states per unit cell', conditions=Value(name='energy', scalars=energy, units='eV'))
174649
def _diff(self, x, th, eps):
        
        foo = lambda y: self.igen(numpy.sum(self.gen(y, th), 0), th)

        out1 = out2 = 0.
        sign = 1 - 2*(x > .5).T
        for I in numpy.ndindex(*((2,)*(len(x)-1)+(1,))):

            eps_ = numpy.array(I)*eps
            x_ = (x.T + sign*eps_).T
            out1 += (-1)**sum(I)*foo(x_)

            x_[-1] = 1
            out2 += (-1)**sum(I)*foo(x_)

        out = out1/out2
        return out
769304
def _dirint_from_dni_ktprime(dni, kt_prime, solar_zenith, use_delta_kt_prime,
                             temp_dew):
    
    times = dni.index
    delta_kt_prime = _delta_kt_prime_dirint(kt_prime, use_delta_kt_prime,
                                            times)
    w = _temp_dew_dirint(temp_dew, times)
    dirint_coeffs = _dirint_coeffs(times, kt_prime, solar_zenith, w,
                                   delta_kt_prime)
    dni_dirint = dni * dirint_coeffs
    return dni_dirint
620566
private double calculateDistanceBetweenTwoAtoms(IAtom atom1, IAtom atom2) {
        double distance;
        Point3d firstPoint = atom1.getPoint3d();
        Point3d secondPoint = atom2.getPoint3d();
        distance = firstPoint.distance(secondPoint);
        return distance;
    }
852224
def rho2theta(self, rho0, gamma):
        
        fac = np.sqrt(np.pi) * special.gamma(1. / 2 * (-1 + gamma)) / special.gamma(gamma / 2.) * 2 / (3 - gamma) * rho0

        #fac = theta_E**(gamma - 1)
        theta_E = fac**(1. / (gamma - 1))
        return theta_E
643812
function getSetJ(h, lw, phi, dec, n, M, L) {

    var w = hourAngle(h, phi, dec),
        a = approxTransit(w, lw, n);
    return solarTransitJ(a, M, L);
}
325777
def StrikeDip(n, e, u):
    
    r2d = 180 / np.pi
    if u < 0:
        n = -n
        e = -e
        u = -u

    strike = np.arctan2(e, n) * r2d
    strike = strike - 90
    while strike >= 360:
        strike = strike - 360
    while strike < 0:
        strike = strike + 360
    x = np.sqrt(np.power(n, 2) + np.power(e, 2))
    dip = np.arctan2(x, u) * r2d
    return (strike, dip)
1326673
def get_element_spd_dos(self, el):
        
        el = get_el_sp(el)
        el_dos = {}
        for site, atom_dos in self.pdos.items():
            if site.specie == el:
                for orb, pdos in atom_dos.items():
                    orbital_type = _get_orb_type_lobster(orb)
                    if orbital_type not in el_dos:
                        el_dos[orbital_type] = pdos
                    else:
                        el_dos[orbital_type] = \
                            add_densities(el_dos[orbital_type], pdos)

        return {orb: Dos(self.efermi, self.energies, densities)
                for orb, densities in el_dos.items()}
690551
def bethe_findfill_zeroT(particles, orbital_e, hopping):
    """"""

    assert 0. <= particles <= len(orbital_e)

    zero = lambda e: np.sum([bethe_filling_zeroT(e-e_m, t) \
                        for t, e_m in zip(hopping, orbital_e)]) - particles
    return fsolve(zero, 0)
230709
def get_Fdot(F, mc, e):
    

    # chirp mass
    mc *= SOLAR2S

    dFdt = 48 / (5*np.pi*mc**2) * (2*np.pi*mc*F)**(11/3) * \
        (1 + 73/24*e**2 + 37/96*e**4) / ((1-e**2)**(7/2))

    return dFdt
719157
def ellipticity2phi_gamma(e1, e2):
    
    phi = np.arctan2(e2, e1)/2
    gamma = np.sqrt(e1**2+e2**2)
    return phi, gamma
1196558
def find_and_read_ebins(hdulist):
    
    from fermipy import utils
    ebins = None
    if 'ENERGIES' in hdulist:
        hdu = hdulist['ENERGIES']
        ectr = hdu.data.field(hdu.columns[0].name)
        ebins = np.exp(utils.center_to_edge(np.log(ectr)))
    elif 'EBOUNDS' in hdulist:
        hdu = hdulist['EBOUNDS']
        emin = hdu.data.field('E_MIN') / 1E3
        emax = hdu.data.field('E_MAX') / 1E3
        ebins = np.append(emin, emax[-1])
    return ebins
124981
def compute_partial_energy(self, removed_indices):
        
        total_energy_matrix = self.total_energy_matrix.copy()
        for i in removed_indices:
            total_energy_matrix[i, :] = 0
            total_energy_matrix[:, i] = 0
        return sum(sum(total_energy_matrix))
500190
def _J(self):
        
        pd = self.particle_distribution(self._Ep * u.GeV)
        return pd.to("1/GeV").value
756248
def calc_z(bands, filling, interaction, hund_cu, name):
    """"""

    while True:
        try:
            data = np.load(name+'.npz')
            break
        except IOError:
            dopout = []
            for dop in filling:
                slsp = Spinon(slaves=2*bands, orbitals=bands, \
                             hopping=[0.5]*2*bands, populations=[dop]*2*bands)
                dopout.append(solve_loop(slsp, interaction, hund_cu)[0][0])
            np.savez(name, zeta=dopout, u_int=interaction, doping=filling, hund=hund_cu)

    return data
537825
def ξ_bounds(self):
        
        λ, Δμ = self._dichalcogenide.λ, self.Δμ
        return (abs(Δμ) - 2 * λ, abs(Δμ))
699384
double getBetheFreeEnergy() {
        // 
        // G_{Bethe} = \sum_a \sum_{x_a} - b(x_a) ln \chi(x_a)
        //              + \sum_a \sum_{x_a} b(x_a) ln b(x_a)
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //           = \sum_a \sum_{x_a} b(x_a) ln (b(x_a) / \chi(x_a))
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //
        //     where n_i is the number of neighbors of the variable x_i,
        //     b(x_a) and b(x_i) are normalized distributions and x_a is 
        //     the set of variables participating in factor a. 
        //
                
        double bethe = 0.0;
        for (int a=0; a<fg.getFactors().size(); a++) {
            Factor f = fg.getFactors().get(a);
            if (!(f instanceof GlobalFactor)) {
                int numConfigs = f.getVars().calcNumConfigs();
                VarTensor beliefs = getFactorBeliefs(a);
                for (int c=0; c<numConfigs; c++) {                
                    // Since we want multiplication by 0 to always give 0 (not the case for Double.POSITIVE_INFINITY or Double.NaN.
                    double b_c = beliefs.getValue(c);
                    if (b_c != s.zero()) {
                        double r_b_c = s.toReal(b_c);
                        double log_b_c = s.toLogProb(b_c);
                        double log_chi_c = f.getLogUnormalizedScore(c);
                        bethe += r_b_c * (log_b_c - log_chi_c);
                    }
                }
            } else {
                VarTensor[] inMsgs = getMsgs(f.getId(), msgs, IN_MSG);
                bethe += ((GlobalFactor) f).getExpectedLogBelief(inMsgs);
            }
        }
        for (int v=0; v<fg.getVars().size(); v++) {
            Var var = fg.getVars().get(v);
            int numNeighbors = bg.numNbsT1(v);
            VarTensor beliefs = getVarBeliefs(v);
            double sum = 0.0;
            for (int c=0; c<var.getNumStates(); c++) {
                double b_c = beliefs.getValue(c);
                if (b_c != s.zero()) {
                    double r_b_c = s.toReal(b_c);
                    double log_b_c = s.toLogProb(b_c);
                    sum += r_b_c * log_b_c;
                }
            }
            bethe -= (numNeighbors - 1) * sum;
        }
        
        assert !Double.isNaN(bethe);        
        return bethe;
    }