Total energy in protons above 1.22 GeV threshold (erg).
973525
def Wp(self):
        
        from scipy.integrate import quad

        Eth = 1.22e-3

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            Wp = quad(
                lambda x: x * self._particle_distribution(x), Eth, np.Inf
            )[0]

        return (Wp * u.TeV).to("erg")
1037164
def _calc_specpp_hiE(self, Egamma):
        
        # Fixed quad with n=40 is about 15 times faster and is always within
        # 0.5% of the result of adaptive quad for Egamma>0.1
        # WARNING: It also produces artifacts for steep distributions (e.g.
        # Maxwellian) at ~500 GeV. Reverting to adaptative quadrature
        # from scipy.integrate import fixed_quad
        # result=c*fixed_quad(self._photon_integrand, 0., 1., args = [Egamma,
        # ], n = 40)[0]
        from scipy.integrate import quad

        Egamma = Egamma.to("TeV").value
        specpp = (
            c.cgs.value
            * quad(
                self._photon_integrand,
                0.0,
                1.0,
                args=Egamma,
                epsrel=1e-3,
                epsabs=0,
            )[0]
        )

        return specpp * u.Unit("1/(s TeV)")
581534
def _finite_well_energy(P, n=1, atol=1e-6):
    '''
    
    '''
    assert n > 0 and n <= _finite_well_states(P)
    pi_2 = pi / 2.
    r = (1 / (P + pi_2)) * (n * pi_2)
    eta = n * pi_2 - arcsin(r) - r * P
    w = 1  # relaxation parameter (for succesive relaxation)
    while True:
        assert r <= 1
        if abs(eta) < atol:
            break
        r2 = r ** 2.
        sqrt_1mr2 = sqrt(1. - r2)
        denom = (1. + P * sqrt_1mr2)
        t1 = P * sqrt_1mr2 / denom * eta
#         t2 = -r * P / (2 * (1. + P * sqrt_1mr2) ** 3) * eta ** 2
        while True:
            next_r = (1 - w) * r + w * (r + t1)
#             next_r = (1 - w) * r + w * (r + t1 + t2)
            next_eta = n * pi_2 - arcsin(next_r) - next_r * P
            # decrease w until eta is converging
            if abs(next_eta / eta) < 1:
                r = next_r
                eta = next_eta
                break
            else:
                w *= 0.5

    alpha = P * r
    E = 2 * (alpha) ** 2  # hbar**2 / (m * L**2)
    return E
1337255
def remnant_lifetime_agb(self):
                '''
                        
                '''
                inim=[]
                remnm=[]
		time11=[]
		tottime=[]
		c_core=[]
		o_core=[]
		small_co_core=[]
		c_core_center=[]	
                for i in range(len(self.runs_H5_surf)):
                        m1p65_last=se(self.runs_H5_out[i])
                        mass_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'mass')
                        top_of_envelope=mass_dummy[len(mass_dummy)-1]
                        h_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'iso_massf','H-1')
			c_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'iso_massf','C-12')
                        o_dummy=m1p65_last.se.get(m1p65_last.se.cycles[len(m1p65_last.se.cycles)-2],'iso_massf','O-16')
			for j in range(len(mass_dummy)):
                                if h_dummy[j] > 1e-1:
                                        bottom_of_envelope = mass_dummy[j]
                                        break
                        inim.append(m1p65_last.get("mini"))
                        remnm.append(bottom_of_envelope)

			###Calculate the lifetime (MS)
                        sefiles=m1p65_last
                        cycs=[]
                        for k in range(5,len(sefiles.se.cycles),5):
                                cycs.append(int(sefiles.se.cycles[k]))
                        w=0
                        for cyc in cycs:
                                c12_center=sefiles.get(cyc,'C-12')[0]
                                #c12_center=c12[w][0]
                                w+=1
                                if c12_center>1e-1:
                                        time1=(sefiles.get(cyc,'age')*sefiles.get('age_unit'))/31557600.
                                        time11.append(time1)
                                        break
			tottime.append(sefiles.get(int(sefiles.se.cycles[-1]),'age')/31557600.)



                print "M_initial | M_remn/bottom of envelope | total lifetime"
                for i in range(len(inim)):
                        print inim[i],"|",'{:.3E}'.format(remnm[i]),"|",'{:.3E}'.format(tottime[i])
54418
private void processElementGamma(GeneratorSingleCluster cluster, Node cur) {
    double k = 1.0;
    double theta = 1.0;
    String kstr = ((Element) cur).getAttribute(ATTR_K);
    if(kstr != null && kstr.length() > 0) {
      k = ParseUtil.parseDouble(kstr);
    }
    String thetastr = ((Element) cur).getAttribute(ATTR_THETA);
    if(thetastr != null && thetastr.length() > 0) {
      theta = ParseUtil.parseDouble(thetastr);
    }

    // *** New normal distribution generator
    Random random = cluster.getNewRandomGenerator();
    Distribution generator = new GammaDistribution(k, theta, random);
    cluster.addGenerator(generator);

    // TODO: check for unknown attributes.
    XMLNodeIterator iter = new XMLNodeIterator(cur.getFirstChild());
    while(iter.hasNext()) {
      Node child = iter.next();
      if(child.getNodeType() == Node.ELEMENT_NODE) {
        LOG.warning("Unknown element in XML specification file: " + child.getNodeName());
      }
    }
  }
367513
def Flavio_to_Fierz_chrom(C, qq, parameters):
    """"""
    p = parameters
    V = ckmutil.ckm.ckm_tree(p["Vus"], p["Vub"], p["Vcb"], p["delta"])
    if qq == 'sb':
        xi = V[2, 2] * V[2, 1].conj()
    elif qq == 'db':
        xi = V[2, 2] * V[2, 0].conj()
    elif qq == 'ds':
        xi = V[2, 1] * V[2, 0].conj()
    elif qq == 'uc':
        xi = V[1, 2].conj() * V[0, 2]
    else:
        raise ValueError("Unexpected flavours: {}".format(qq))
    qqfl = qq[::-1]
    e = sqrt(4 * pi * parameters['alpha_e'])
    gs = sqrt(4 * pi * parameters['alpha_s'])
    if qq == 'sb' or qq == 'db':
        mq = parameters['m_b']
    elif qq == 'ds':
        mq = parameters['m_s']
    elif qq == 'uc':
        mq = parameters['m_c']
    else:
        KeyError("Not sure what to do with quark mass for flavour {}".format(qq))
    dic = {
        'F7gamma' + qq: C["C7_" + qqfl] / ((16 * pi**2) / e / mq),
        'F8g' + qq: C["C8_" + qqfl] / ((16 * pi**2) / gs / mq),
        'F7pgamma' + qq: C["C7p_" + qqfl] / ((16 * pi**2) / e / mq),
        'F8pg' + qq: C["C8p_" + qqfl] / ((16 * pi**2) / gs / mq)
            }
    prefactor = sqrt(2)/p['GF']/xi/4
    return {k: v / prefactor for k, v in dic.items()}
882588
def total_energy(self):
        
        return sum(sum(self._recip)) + sum(sum(self._real)) + sum(self._point) + self._charged_cell_energy
304911
def _psi2computations(dL_dpsi2, variance, Z, mu, S, gamma):
    
    # here are the "statistics" for psi1 and psi2
    # Produced intermediate results:
    # _psi2_dvariance      Q
    # _psi2_dZ             MxQ
    # _psi2_dgamma         NxQ
    # _psi2_dmu            NxQ
    # _psi2_dS             NxQ
    
    mu2 = np.square(mu)
    gamma2 = np.square(gamma)
    variance2 = np.square(variance)
    mu2S = mu2+S # NxQ
    gvm = gamma*mu*variance
    common_sum = gvm.dot(Z.T)
    Z_expect = (np.dot(dL_dpsi2,Z)*Z).sum(axis=0)
    Z_expect_var2 = Z_expect*variance2
    dL_dpsi2T = dL_dpsi2+dL_dpsi2.T
    common_expect = common_sum.dot(dL_dpsi2T).dot(Z)
    Z2_expect = common_sum.dot(dL_dpsi2T)
    Z1_expect = dL_dpsi2T.dot(Z)
    
    dL_dvar = variance*Z_expect*2.*(gamma*mu2S-gamma2*mu2).sum(axis=0)+(common_expect*gamma*mu).sum(axis=0)
        
    dL_dgamma = Z_expect_var2*(mu2S-2.*gamma*mu2)+common_expect*mu*variance
                
    dL_dmu = Z_expect_var2*mu*2.*(gamma-gamma2) + common_expect*gamma*variance

    dL_dS = gamma*Z_expect_var2
    
    dL_dZ = (gamma*(mu2S-gamma*mu2)).sum(axis=0)*variance2*Z1_expect+ Z2_expect.T.dot(gamma*mu)*variance

    return dL_dvar, dL_dgamma, dL_dmu, dL_dS, dL_dZ
745503
def electric_field_amplitude_intensity(s0,Omega=1.0e6):
	''''''
	e0=hbar*Omega/(e*a0) #This is the electric field scale.
	
	I0=2.50399 #mW/cm^2
	I0=1.66889451102868 #mW/cm^2
	
	I0=I0/1000*(100**2) #W/m^2
	r_ciclic=4.226983616875483 #a0
	gamma_D2=2*Pi*6.065e6/Omega # The decay frequency of the D2 line.
	E0_sat=gamma_D2/r_ciclic/sqrt(2.0)

	E0_sat=E0_sat*e0
	I0=E0_sat**2/2/c/mu0

    #return sqrt(c*mu0*s0*I0/2)/e0
	#return sqrt(c*mu0*s0*I0)/e0
	return sqrt(2*c*mu0*s0*I0)/e0
69830
def read_excitation_energies(self):
        

        transitions = list()

        # read in file
        with zopen(self.filename, "r") as f:
            line = f.readline()
            td = False
            while line != "":
                if re.search(r"^\sExcitation energies and oscillator strengths:", line):
                    td = True

                if td:
                    if re.search(r"^\sExcited State\s*\d", line):
                        val = [float(v) for v in float_patt.findall(line)]
                        transitions.append(tuple(val[0:3]))
                line = f.readline()
        return transitions
1211625
def _emiss_ep(self, Eph):
        
        if self.weight_ep == 0.0:
            return np.zeros_like(Eph)

        gam = np.vstack(self._gam)
        eps = (Eph / mec2).decompose().value
        # compute integral with electron distribution
        emiss = c.cgs * trapz_loglog(
            np.vstack(self._nelec) * self._sigma_ep(gam, eps),
            self._gam,
            axis=0,
        ).to(u.cm ** 2 / Eph.unit)
        return emiss
1241362
def energ_orb(m1,m2,r):
    ''' 
    
        
    '''
    epo = -grav_const * m1 * m2 * msun_g**2 / (r * rsun_cm)
    return epo
1060027
def radiance2tb(self, rad):
        
        #
        # Tb = C2 * νc/{α * log[C1*νc**3 / L + 1]} - β/α
        #
        # C1 = 2 * h * c**2 and C2 = hc/k
        #
        c_1 = 2 * H_PLANCK * C_SPEED ** 2
        c_2 = H_PLANCK * C_SPEED / K_BOLTZMANN

        vc_ = SEVIRI[self.bandname][self.platform_name][0]
        # Multiply by 100 to get SI units!
        vc_ *= 100.0
        alpha = SEVIRI[self.bandname][self.platform_name][1]
        beta = SEVIRI[self.bandname][self.platform_name][2]

        tb_ = c_2 * vc_ / \
            (alpha * np.log(c_1 * vc_ ** 3 / rad + 1)) - beta / alpha

        return tb_
671188
def _process_input_seed(seed_photon_fields):
        

        Tcmb = 2.72548 * u.K  # 0.00057 K
        Tfir = 30 * u.K
        ufir = 0.5 * u.eV / u.cm ** 3
        Tnir = 3000 * u.K
        unir = 1.0 * u.eV / u.cm ** 3

        # Allow for seed_photon_fields definitions of the type 'CMB-NIR-FIR' or
        # 'CMB'
        if type(seed_photon_fields) != list:
            seed_photon_fields = seed_photon_fields.split("-")

        result = OrderedDict()

        for idx, inseed in enumerate(seed_photon_fields):
            seed = {}
            if isinstance(inseed, six.string_types):
                name = inseed
                seed["type"] = "thermal"
                if inseed == "CMB":
                    seed["T"] = Tcmb
                    seed["u"] = ar * Tcmb ** 4
                    seed["isotropic"] = True
                elif inseed == "FIR":
                    seed["T"] = Tfir
                    seed["u"] = ufir
                    seed["isotropic"] = True
                elif inseed == "NIR":
                    seed["T"] = Tnir
                    seed["u"] = unir
                    seed["isotropic"] = True
                else:
                    log.warning(
                        "Will not use seed {0} because it is not "
                        "CMB, FIR or NIR".format(inseed)
                    )
                    raise TypeError
            elif type(inseed) == list and (
                len(inseed) == 3 or len(inseed) == 4
            ):
                isotropic = len(inseed) == 3

                if isotropic:
                    name, T, uu = inseed
                    seed["isotropic"] = True
                else:
                    name, T, uu, theta = inseed
                    seed["isotropic"] = False
                    seed["theta"] = validate_scalar(
                        "{0}-theta".format(name), theta, physical_type="angle"
                    )

                thermal = T.unit.physical_type == "temperature"

                if thermal:
                    seed["type"] = "thermal"
                    validate_scalar(
                        "{0}-T".format(name),
                        T,
                        domain="positive",
                        physical_type="temperature",
                    )
                    seed["T"] = T
                    if uu == 0:
                        seed["u"] = ar * T ** 4
                    else:
                        # pressure has same physical type as energy density
                        validate_scalar(
                            "{0}-u".format(name),
                            uu,
                            domain="positive",
                            physical_type="pressure",
                        )
                        seed["u"] = uu
                else:
                    seed["type"] = "array"
                    # Ensure everything is in arrays
                    T = u.Quantity((T,)).flatten()
                    uu = u.Quantity((uu,)).flatten()

                    seed["energy"] = validate_array(
                        "{0}-energy".format(name),
                        T,
                        domain="positive",
                        physical_type="energy",
                    )

                    if np.isscalar(seed["energy"]) or seed["energy"].size == 1:
                        seed["photon_density"] = validate_scalar(
                            "{0}-density".format(name),
                            uu,
                            domain="positive",
                            physical_type="pressure",
                        )
                    else:
                        if uu.unit.physical_type == "pressure":
                            uu /= seed["energy"] ** 2
                        seed["photon_density"] = validate_array(
                            "{0}-density".format(name),
                            uu,
                            domain="positive",
                            physical_type="differential number density",
                        )
            else:
                raise TypeError(
                    "Unable to process seed photon"
                    " field: {0}".format(inseed)
                )

            result[name] = seed

        return result
270140
def main():
    
    dataset = 'RD'

    observables = ['log(<\sigma v>)', '\Omega_{\chi}h^2', 'log(\sigma_p^{SI})']
    var  = ['log(m_{\chi})']
    var += ['log(C_1)', 'log(C_2)', 'log(C_3)', 'log(C_4)', 'log(C_5)', 'log(C_6)']
    var += observables

    plot_vs_mass(dataset, observables, 'mass_vs_observables.png')
    plot_oneD(dataset, var, 'oneD.png')
    pairplot(dataset, var, 'pairplot.png')
1217436
def _iso_ic_on_monochromatic(
        electron_energy, seed_energy, seed_edensity, gamma_energy
    ):
        
        photE0 = (seed_energy / mec2).decompose().value
        phn = seed_edensity

        # electron_energy = electron_energy[:, None]
        gamma_energy = gamma_energy[:, None]
        photE0 = photE0[:, None, None]
        phn = phn[:, None, None]

        b = 4 * photE0 * electron_energy
        w = gamma_energy / electron_energy
        q = w / (b * (1 - w))
        fic = (
            2 * q * np.log(q)
            + (1 + 2 * q) * (1 - q)
            + (1.0 / 2.0) * (b * q) ** 2 * (1 - q) / (1 + b * q)
        )

        gamint = (
            fic
            * heaviside(1 - q)
            * heaviside(q - 1.0 / (4 * electron_energy ** 2))
        )
        gamint[np.isnan(gamint)] = 0.0

        if phn.size > 1:
            phn = phn.to(1 / (mec2_unit * u.cm ** 3)).value
            gamint = trapz_loglog(gamint * phn / photE0, photE0, axis=0)  # 1/s
        else:
            phn = phn.to(mec2_unit / u.cm ** 3).value
            gamint *= phn / photE0 ** 2
            gamint = gamint.squeeze()

        # gamint /= mec2.to('erg').value

        # r0 = (e**2 / m_e / c**2).to('cm')
        # sigt = ((8 * np.pi) / 3 * r0**2).cgs
        sigt = 6.652458734983284e-25
        c = 29979245800.0

        gamint *= (3.0 / 4.0) * sigt * c / electron_energy ** 2

        return gamint
43802
def _photon_integrand(self, x, Egamma):
        
        try:
            return (
                self._sigma_inel(Egamma / x)
                * self._particle_distribution((Egamma / x))
                * self._Fgamma(x, Egamma / x)
                / x
            )
        except ZeroDivisionError:
            return np.nan