trapz integrates over altitude axis, axis = -2
    concatenate over reaction dimension, axis = -1

    br: column integrated brightness
    lamb: wavelength [nm]
    ver: volume emission rate  [photons / cm^-3 s^-3 ...]
485361
def catvl(z, ver, vnew, lamb, lambnew, br):
    
    if ver is not None:
        br = np.concatenate((br, np.trapz(vnew, z, axis=-2)), axis=-1)  # must come first!
        ver = np.concatenate((ver, vnew), axis=-1)
        lamb = np.concatenate((lamb, lambnew))
    else:
        ver = vnew.copy(order='F')
        lamb = lambnew.copy()
        br = np.trapz(ver, z, axis=-2)

    return ver, lamb, br
666731
def densu(alt, dlb, tinf, tlb, xm, alpha, tz, zlb, s2, mn1, zn1, tn1, tgn1):
    '''

 '''
    rgas = 831.4
    #rgas = 831.44621    #maybe make this a global constant?
    densu_temp = 1.0

    xs = [0.0]*5
    ys = [0.0]*5
    y2out = [0.0]*5

    #/* joining altitudes of Bates and spline */
    za=zn1[0];
    if (alt>za):
        z=alt;
    else:
        z=za;

    #/* geopotential altitude difference from ZLB */
    zg2 = zeta(z, zlb);
    
    #/* Bates temperature */
    tt = tinf - (tinf - tlb) * exp(-s2*zg2);
    ta = tt;
    tz[0] = tt
    densu_temp = tz[0]

    if (alt<za):
        #/* calculate temperature below ZA
        # * temperature gradient at ZA from Bates profile */
        dta = (tinf - ta) * s2 * pow(((re[0]+zlb)/(re[0]+za)),2.0);
        tgn1[0]=dta;
        tn1[0]=ta;
        if (alt>zn1[mn1-1]):
            z=alt;
        else:
            z=zn1[mn1-1];
        mn=mn1;
        z1=zn1[0];
        z2=zn1[mn-1];
        t1=tn1[0];
        t2=tn1[mn-1];
        #/* geopotental difference from z1 */
        zg = zeta (z, z1);
        zgdif = zeta(z2, z1);
        #/* set up spline nodes */
        for k in range(mn):
            xs[k] = zeta(zn1[k], z1) / zgdif;
            ys[k] = 1.0 / tn1[k];
        
        #/* end node derivatives */
        yd1 = -tgn1[0] / (t1*t1) * zgdif;
        yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re[0]+z2)/(re[0]+z1)),2.0);
        #/* calculate spline coefficients */
        spline (xs, ys, mn, yd1, yd2, y2out);
        x = zg / zgdif;
        y = [0.0]
        splint (xs, ys, y2out, mn, x, y);
        #/* temperature at altitude */
        tz[0] = 1.0 / y[0];
        densu_temp = tz[0];

    if (xm==0):
        return densu_temp;

    #/* calculate density above za */
    glb = gsurf[0] / pow((1.0 + zlb/re[0]),2.0);
    gamma = xm * glb / (s2 * rgas * tinf);
    expl = exp(-s2 * gamma * zg2);
    if (expl>50.0): # pragma: no cover
        expl=50.0;
    if (tt<=0): # pragma: no cover
        expl=50.0;

    #/* density at altitude */
    densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl;
    densu_temp=densa;
    if (alt>=za):
        return densu_temp;

    #/* calculate density below za */
    glb = gsurf[0] / pow((1.0 + z1/re[0]),2.0);
    gamm = xm * glb * zgdif / rgas;

    #/* integrate spline temperatures */
    yi = [0]
    splini (xs, ys, y2out, mn, x, yi);
    expl = gamm * yi[0];
    if (expl>50.0): # pragma: no cover
        expl=50.0;
    if (tz[0]<=0): # pragma: no cover
        expl=50.0;

    #/* density at altitude */
    densu_temp = densu_temp * pow ((t1 / tz[0]),(1.0 + alpha)) * exp(-expl);
    return densu_temp;
650118
def _differential_to_integral(hpx_map):
        
        ebins = hpx_map.hpx.ebins
        ratio = ebins[1:] / ebins[0:-1]
        half_log_ratio = np.log(ratio) / 2.
        int_map = ((hpx_map.data[0:-1].T * ebins[0:-1]) +
                   (hpx_map.data[1:].T * ebins[1:])) * half_log_ratio
        return HpxMap(int_map.T, hpx_map.hpx)
907615
def luminosity_to_flux(lumErg_S, dist_Mpc):
    
    ################ > IMPORTS ################
    ## STANDARD LIB ##
    ## THIRD PARTY ##
    import numpy as np
    import math
    ## LOCAL APPLICATION ##

    ################ > VARIABLE SETTINGS ######

    ################ >ACTION(S) ################
    # Convert the distance to cm
    distCm = dist_Mpc * MPC_2_CMS
    fluxErg_cm2_S = lumErg_S / (4 * np.pi * distCm ** 2)

    return fluxErg_cm2_S
1270774
def perform_es_corr(self, lattice, q, step=1e-4):
        
        logger.info("Running Freysoldt 2011 PC calculation (should be " "equivalent to sxdefectalign)")
        logger.debug("defect lattice constants are (in angstroms)" + str(lattice.abc))

        [a1, a2, a3] = ang_to_bohr * np.array(lattice.get_cartesian_coords(1))
        logging.debug("In atomic units, lat consts are (in bohr):" + str([a1, a2, a3]))
        vol = np.dot(a1, np.cross(a2, a3))  # vol in bohr^3

        def e_iso(encut):
            gcut = eV_to_k(encut)  # gcut is in units of 1/A
            return scipy.integrate.quad(lambda g: self.q_model.rho_rec(g * g)**2, step, gcut)[0] * (q**2) / np.pi

        def e_per(encut):
            eper = 0
            for g2 in generate_reciprocal_vectors_squared(a1, a2, a3, encut):
                eper += (self.q_model.rho_rec(g2)**2) / g2
            eper *= (q**2) * 2 * round(np.pi, 6) / vol
            eper += (q**2) * 4 * round(np.pi, 6) \
                * self.q_model.rho_rec_limit0 / vol
            return eper

        eiso = converge(e_iso, 5, self.madetol, self.energy_cutoff)
        logger.debug("Eisolated : %f", round(eiso, 5))

        eper = converge(e_per, 5, self.madetol, self.energy_cutoff)

        logger.info("Eperiodic : %f hartree", round(eper, 5))
        logger.info("difference (periodic-iso) is %f hartree", round(eper - eiso, 6))
        logger.info("difference in (eV) is %f", round((eper - eiso) * hart_to_ev, 4))

        es_corr = round((eiso - eper) / self.dielectric * hart_to_ev, 6)
        logger.info("Defect Correction without alignment %f (eV): ", es_corr)
        return es_corr
1072706
def mzminus(df, minus=0, noise=10000):
    
    mol_ions = ((df.values > noise) * df.columns).max(axis=1) - minus
    mol_ions[np.abs(mol_ions) < 0] = 0
    d = np.abs(np.ones(df.shape) * df.columns -
               (mol_ions[np.newaxis].T * np.ones(df.shape))) < 1
    d = (df.values * d).sum(axis=1)
    return Trace(d, df.index, name='m-' + str(minus))
684217
def detect_plates_vzcheck(step, seuil_memz):
    """"""
    v_z = step.fields['v3'][0, :, :, 0]
    v_x = step.fields['v2'][0, :, :, 0]
    tcell = step.fields['T'][0, :, :, 0]
    n_z = step.geom.nztot
    nphi = step.geom.nptot  # -1? should be OK, ghost not included
    rcmb = max(0, step.geom.rcmb)
    radius = step.geom.r_coord + rcmb
    radiusgrid = step.geom.rgeom[:, 0] + rcmb
    dphi = 1 / nphi

    # calculing temperature on the grid and vz_mean
    vz_mean = 0
    tgrid = np.zeros((nphi, n_z + 1))
    tgrid[:, 0] = 1
    for i_z in range(1, n_z):
        for phi in range(nphi):
            tgrid[phi, i_z] = (
                tcell[phi, i_z - 1] *
                (radiusgrid[i_z] - radius[i_z - 1]) + tcell[phi, i_z] *
                (-radiusgrid[i_z] + radius[i_z])) / (radius[i_z] -
                                                     radius[i_z - 1])
            vz_mean += abs(v_z[phi, i_z]) / (nphi * n_z)

    flux_c = n_z * [0]
    for i_z in range(1, n_z - 1):
        for phi in range(nphi):
            flux_c[i_z] += (tgrid[phi, i_z] - step.timeinfo.loc['Tmean']) * \
                v_z[phi, i_z] * radiusgrid[i_z] * dphi

    # checking stagnant lid
    if all(abs(flux_c[i_z]) <= np.max(flux_c) / 50
           for i_z in range(n_z - n_z // 20, n_z)):
        raise error.StagnantLidError(step.sdat)
    else:
        # verifying horizontal plate speed and closeness of plates
        dvphi = nphi * [0]
        dvx_thres = 16 * step.timeinfo.loc['vrms']

        for phi in range(0, nphi):
            dvphi[phi] = (v_x[phi, n_z - 1] -
                          v_x[phi - 1, n_z - 1]) / ((1 + rcmb) * dphi)
        limits = []
        for phi in range(0, nphi - nphi // 33):
            mark = all(abs(dvphi[i]) <= abs(dvphi[phi])
                       for i in range(phi - nphi // 33, phi + nphi // 33))
            if mark and abs(dvphi[phi]) >= dvx_thres:
                limits.append(phi)
        for phi in range(nphi - nphi // 33 + 1, nphi):
            mark = all(abs(dvphi[i]) <= abs(dvphi[phi])
                       for i in range(phi - nphi // 33 - nphi,
                                      phi + nphi // 33 - nphi))
            if mark and abs(dvphi[phi]) >= dvx_thres:
                limits.append(phi)

        # verifying vertical speed
        k = 0
        for i in range(len(limits)):
            vzm = 0
            phi = limits[i - k]
            for i_z in range(1 if phi == nphi - 1 else 0, n_z):
                vzm += (abs(v_z[phi, i_z]) +
                        abs(v_z[phi - 1, i_z]) +
                        abs(v_z[(phi + 1) % nphi, i_z])) / (n_z * 3)

            vz_thres = vz_mean * 0.1 + seuil_memz / 2 if seuil_memz else 0
            if vzm < vz_thres:
                limits.remove(phi)
                k += 1

    return limits, nphi, dvphi, vz_thres, v_x[:, n_z - 1]
1164955
def ccor2(alt, r, h1, zh, h2):
    '''

 '''
    e1 = (alt - zh) / h1;
    e2 = (alt - zh) / h2;
    if ((e1 > 70.0) or (e2 > 70)): # pragma: no cover
        return 1.0 # exp(0)
    if ((e1 < -70) and (e2 < -70)): # pragma: no cover
        return exp(r)
    ex1 = exp(e1);
    ex2 = exp(e2);
    ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
    return exp(ccor2v);
1161060
def ccor(alt, r, h1, zh):
    
    e = (alt - zh) / h1
    if(e>70.0):
        return 1.0 # exp(0) # pragma: no cover
    elif (e < -70.0):
        return exp(r)
    ex = exp(e)
    e = r / (1.0 + ex)
    return exp(e)
963895
def get_SZ_orient(self):
        

        tm_outdated = self._tm_signature != (self.radius, self.radius_type, 
            self.wavelength, self.m, self.axis_ratio, self.shape, self.ddelt, 
            self.ndgs)
        scatter_outdated = self._scatter_signature != (self.thet0, self.thet, 
            self.phi0, self.phi, self.alpha, self.beta, self.orient)

        orient_outdated = self._orient_signature != \
            (self.orient, self.or_pdf, self.n_alpha, self.n_beta)
        if orient_outdated:
            self._init_orient()
        
        outdated = tm_outdated or scatter_outdated or orient_outdated

        if outdated:
            (self._S_orient, self._Z_orient) = self.orient(self)
            self._set_scatter_signature()

        return (self._S_orient, self._Z_orient)
496286
def _get_v(self, rho):
        
        v_mol = self.mass / rho  # in cm^3
        v = vol_mol2uc(v_mol * 1.e-6, self.z)
        print(v)
        return v
1335550
def example_ylm(m=0, n=2, shape=128, limits=[-4, 4], draw=True, show=True, **kwargs):
    """"""
    import ipyvolume.pylab as p3

    __, __, __, r, theta, phi = xyz(shape=shape, limits=limits, spherical=True)
    radial = np.exp(-(r - 2) ** 2)
    data = np.abs(scipy.special.sph_harm(m, n, theta, phi) ** 2) * radial  # pylint: disable=no-member
    if draw:
        vol = p3.volshow(data=data, **kwargs)
        if show:
            p3.show()
        return vol
    else:
        return data
1130852
def int_flux_threshold(self, skydir, fn, ts_thresh, min_counts):
        

        ebins = 10**np.linspace(np.log10(self.ebins[0]),
                                np.log10(self.ebins[-1]), 33)
        ectr = np.sqrt(ebins[0] * ebins[-1])

        sig, bkg, bkg_fit = self.compute_counts(skydir, fn, ebins)

        norms = irfs.compute_norm(sig, bkg, ts_thresh,
                                  min_counts, sum_axes=[1, 2, 3], bkg_fit=bkg_fit,
                                  rebin_axes=[4, 10, 1])

        npred = np.squeeze(np.apply_over_axes(np.sum, norms * sig, [1, 2, 3]))
        npred = np.array(npred, ndmin=1)
        flux = np.squeeze(norms) * fn.flux(ebins[0], ebins[-1])
        eflux = np.squeeze(norms) * fn.eflux(ebins[0], ebins[-1])
        dnde = np.squeeze(norms) * fn.dnde(ectr)
        e2dnde = ectr**2 * dnde

        o = dict(e_min=self.ebins[0], e_max=self.ebins[-1], e_ref=ectr,
                 npred=npred, flux=flux, eflux=eflux,
                 dnde=dnde, e2dnde=e2dnde)

        sig, bkg, bkg_fit = self.compute_counts(skydir, fn)

        npred = np.squeeze(np.apply_over_axes(np.sum, norms * sig,
                                              [2, 3]))
        flux = np.squeeze(np.squeeze(norms, axis=(1, 2, 3))[:, None] *
                          fn.flux(self.ebins[:-1], self.ebins[1:]))
        eflux = np.squeeze(np.squeeze(norms, axis=(1, 2, 3))[:, None] *
                           fn.eflux(self.ebins[:-1], self.ebins[1:]))
        dnde = np.squeeze(np.squeeze(norms, axis=(1, 2, 3))
                          [:, None] * fn.dnde(self.ectr))
        e2dnde = ectr**2 * dnde

        o['bins'] = dict(npred=npred,
                         flux=flux,
                         eflux=eflux,
                         dnde=dnde,
                         e2dnde=e2dnde,
                         e_min=self.ebins[:-1], e_max=self.ebins[1:],
                         e_ref=self.ectr)

        return o
685628
def _evaluate_xyz(self,x,y,z):
        """"""
        return 2.*numpy.pi*self._b*self._c\
            *_potInt(x,y,z,self._psi,
                     self._b2,self._c2,glx=self._glx,glw=self._glw)
1156615
def emitResetAxisSignal(self, axisNumber):
        
        assert axisNumber in (VALID_AXES_NUMBERS), \
            "Axis Nr should be one of {}, got {}".format(VALID_AXES_NUMBERS, axisNumber)

        # Hide 'auto-scale (A)' button
        logger.debug("ArgosPgPlotItem.autoBtnClicked, mode:{}".format(self.autoBtn.mode))
        if self.autoBtn.mode == 'auto':
            self.autoBtn.hide()
        else:
            # Does this occur?
            msg = "Unexpected autobutton mode: {}".format(self.autoBtn.mode)
            if DEBUGGING:
                raise ValueError(msg)
            else:
                logger.warn(msg)

        logger.debug("Emitting sigAxisReset({}) for {!r}".format(axisNumber, self))
        self.sigAxisReset.emit(axisNumber)
1262433
def emissive_part_3x(self, tb=True):
        """"""
        try:
            # Emissive part:
            self._e3x = self._rad3x_t11 * (1 - self._r3x)
            # Unsure how much sense it makes to apply the co2 correction term here!?
            # FIXME!
            # self._e3x *= self._rad3x_correction

        except TypeError:
            LOG.warning(
                "Couldn't derive the emissive part \n" +
                "Please derive the relfectance prior to requesting the emissive part")

        if tb:
            return self.radiance2tb(self._e3x)
        else:
            return self._e3x
428609
def herz_me(val):
    """"""
    result = 0
    if val.endswith("MHz"):
        stripped = val.replace("MHz", "")
        strip_fl = float(stripped)
        result = strip_fl * 1000000
    elif val.endswith("kHz"):
        stripped = val.replace("kHz", "")
        strip_fl = float(stripped)
        result = strip_fl * 1000
    elif val.endswith("Hz"):
        stripped = val.replace("Hz", "")
        result = float(stripped)
    return(result)