def simulateBackward(self, direction=1):
		'''
		Propagate the signal in backward direction using the population
		found in the previous forward iteration. Since N2 and N1 are constant
		we can solve each equations with a simple integration
		'''		

		# Get the initiale conditions
		Pp_ini = self.P_p_b[:,-1]
		Ps_ini = self.P_s_b[:,-1]
		Pase_ini = self.P_ase_b[:,-1]
		
		self.invSptProfil()

		for m in arange(self.nbrPump):
			integrant = sign(direction)*(-self.sigma_abs_p[m]*self.N1[::-1] - self.alpha_p) * self.dopedFiber.pumpOverlap(self.pumpWL[m])
			self.P_p_b[m,::-1] = r_[Pp_ini[m], Pp_ini[m]*exp(integrate.cumtrapz(integrant, self.z))]

		for l in arange(self.nbrSignal):
			integrant = sign(direction)*(self.sigma_em_s[l]*self.N2[::-1] - self.sigma_abs_s[l]*self.N1[::-1] - self.alpha_s)
			integrant *= self.dopedFiber.modeOverlap(self.signalWL[l], self.sigDC[l])
			self.P_s_b[l,::-1] = r_[Ps_ini[l], Ps_ini[l]*exp(integrate.cumtrapz(integrant, self.z))]

		for v in arange(self.nbrAse):
			integrant = sign(direction)*(self.sigma_em_ase[v]*self.N2[::-1] - self.sigma_abs_ase[v]*self.N1[::-1] - self.alpha_ase)
			integrant *= self.dopedFiber.modeOverlap(self.aseWL[v])
			integrant2 = sign(direction)*2*(h*c/(self.aseWL[v]*1E-6)) * self.delta_nu[v] * self.sigma_em_ase[v]*self.N2[::-1]
			integrant2 *= self.dopedFiber.modeOverlap(self.aseWL[v])

			sol = integrate.cumtrapz(integrant, self.z)
			solTerme1 = exp(sol)
			solTerme1b = r_[1.0, exp(-sol)]
			solTerme2 = solTerme1 * integrate.cumtrapz(integrant2*solTerme1b, self.z)
			self.P_ase_b[v,::-1] = r_[Pase_ini[v], Pase_ini[v]*solTerme1 + solTerme2]

 def amin_i( self, demsi, depsf0, depsf_smaller, a1, damage, idx ):
         if a1 > np.max( self.sorted_l_ez[self.sf_mask] ):
             amin = 1. / ( demsi + depsf0 ) * ( ( demsi + depsf0 ) * ( demsi + depsf_smaller ) ) ** ( .5 ) * a1
         else:
             l_ez_behind_a1 = self.sorted_l_ez[idx + 1:] 
             a = np.linspace( a1, 1.2 * self.amin_i_guess  , 500 )
             Kf = self.Kf
             a_shaped = a.reshape( len( a ), 1 )
             l_ez_behind_a1 = l_ez_behind_a1.reshape( 1, len( l_ez_behind_a1 ) )
             m_p = a_shaped / l_ez_behind_a1
             p = np.abs( H( 1 - m_p ) * ( 1 - m_p ) )
             ####
             muT = np.sum( self.sorted_depsf[idx + 1:] * ( 1 - damage[idx + 1:] ) * Kf[idx + 1:] * p, 1 )
             Kf_intact = np.sum( ( 1 - damage[idx + 1:] ) * Kf[idx + 1:] * ( 1 - p ) , 1 ) + np.sum( ( 1 - damage[:idx + 1] ) * Kf[:idx + 1] ) 
             Kf_broken = np.sum( Kf * damage )
             Emtrx = ( 1. - self.V_f_tot ) * self.E_m + Kf_broken + Kf_intact
             depsm = muT / Emtrx
             em = np.hstack( ( 0, cumtrapz( depsm, a ) ) )
             um = np.hstack( ( 0, cumtrapz( em , a ) ) )
             condI = a ** 2. / 2. * depsf0 + a * em - um - a1 ** 2 / 2.*depsf_smaller
             ip_f = interp1d( condI, a, bounds_error = False, fill_value = a[-1] )
             cut_at = np.nonzero( H( ip_f( 0 ) - a ) )
             amin = np.hstack( ( a[cut_at], ip_f( 0 ) ) )
             amin = amin[-1]
             self.amin_i_guess = amin
         return amin

def prepare_mf(mpart, grid, mf_kwargs):
    M = np.linspace(np.log10(mpart), np.log10(grid.max()), 2000)
    mf_obj = MassFunction(M=M, **mf_kwargs)

    mf = mf_obj.dndm
    m_outside_range = mf_obj.mltm[0] + mf_obj.mgtm[-1]

    cumfunc = cumtrapz(10 ** M * mf, M, initial=0) * np.log(10)

    cdf = spline(M, cumfunc, k=3)
    icdf = spline(cumfunc, M, k=3)

    if MAKE_PLOTS:
        plt.clf()
        plt.plot(M, cumfunc)

        plt.plot(M, cdf(M))
        plt.savefig("cumfunc.pdf")

        plt.clf()
        mcumfunc = cumtrapz(10 ** (2 * M) * mf, dx=M[1] - M[0], initial=1e-20) * np.log(10)
        plt.plot(M, mcumfunc)
        plt.savefig("mcumfunc.pdf")

        # How much mass is above 10**12.5
        minvcumfunc = cumtrapz(10 ** (2 * M[::-1]) * mf[::-1], dx=M[1] - M[0]) * np.log(10)
        minvcumfunc = np.concatenate((np.array([minvcumfunc[0]]), minvcumfunc))
        minvcumfunc /= minvcumfunc[-1]
        plt.clf()
        plt.plot(M, minvcumfunc[::-1])
        plt.yscale('log')
        plt.grid(True)
        plt.savefig("minvcumfunc.pdf")

    return cdf, icdf, M, mf, m_outside_range

def exchange(a, b, grid, end_point=-1):
    abc, asc, al, aj = a  # radial and angular for a
    bbc, bsc, bl, bj = b  # radial and angular for b
    jmax, jmin = max(aj, bj), min(aj, bj)
    ro, w, h = grid  # radial grid and weights
    end_ind = end_point < 0 and len(ro) or (sc.where(ro > end_point)[0][0]+1)
    cabc = abc[:end_ind]
    casc = asc[:end_ind]
    cbbc = bbc[:end_ind]
    cbsc = bsc[:end_ind]
    cro = ro[:end_ind]
    cw = w[:end_ind]
    k = jmax - jmin
    if (k+al+bl) % 2 != 0:
        k += 1
    dens0 = cabc*cbbc+casc*cbsc

    dens_kg = dens0/cro**(k+1)
    dens_kl = dens0*cro**k
    res = 0e0
    while (k <= (jmax+jmin)):
        dens_g_int = cumtrapz(dens_kg*cw, initial=0e0)*h
        # внешнее интегрирование от r до бесконечности
        dens_g_int = dens_g_int[-1] - dens_g_int
        dens_l_int = cumtrapz(dens_kl*cw, initial=0e0)*h
        res_k = sc.trapz((dens_kg*dens_l_int+dens_kl*dens_g_int)*cw)*h
        res += res_k*w3js0[(jmax, jmin, k)]
        k += 2
        dens_kg = dens_kg/cro**2
        dens_kl = dens_kl*cro**2
    return res

 def test_kernels_stddev(self):
     """
     Test that the standard deviation calculated from the kernel (almost)
     equals the parameter sigma with which the kernel was constructed.
     """
     sigma = 0.5 * pq.s
     kernel_resolution = sigma / 50.0
     for invert in (False, True):
         kernel_list = [kernel_type(sigma, invert) for
                        kernel_type in self.kernel_types]
         for kernel in kernel_list:
             b = kernel.boundary_enclosing_area_fraction(self.fraction).magnitude
             restric_defdomain = \
                 np.linspace(-b, b, 2*b/kernel_resolution.magnitude) * \
                 sigma.units
             kern = kernel(restric_defdomain)
             av_integr = kern * restric_defdomain
             average = spint.cumtrapz(y=av_integr.magnitude,
                                      x=restric_defdomain.magnitude)[-1] * \
                       sigma.units
             var_integr = (restric_defdomain-average)**2 * kern
             variance = spint.cumtrapz(y=var_integr.magnitude,
                                       x=restric_defdomain.magnitude)[-1] * \
                        sigma.units**2
             stddev = np.sqrt(variance)
             self.assertAlmostEqual(stddev, sigma, delta=0.01*sigma)

 def str_str(self, f):
     """
     Gets str_str.out file and analyze it
     """
     lines = f.read()
     lines = lines.split('\n')
     EVM, SVM, E11, E22, E33, S11, S22 = [], [], [], [], [], [], []
     nline = len(lines) - 1
     for i in range(nline + 5):
         try:
             EVM.append(float(lines[i+1].split()[0]))
             SVM.append(float(lines[i+1].split()[1]))
             E11.append(float(lines[i+1].split()[2]))
             E22.append(float(lines[i+1].split()[3]))
             E33.append(float(lines[i+1].split()[4]))
             S11.append(float(lines[i+1].split()[8]))
             S22.append(float(lines[i+1].split()[9]))
         except IndexError:
             break
     workx = integrate.cumtrapz(y=S11, x=E11)
     worky = integrate.cumtrapz(y=S22, x=E22)
     workTotal = []
     for k in range(len(workx)):
         workTotal.append(workx[k]+worky[k])
     return S11, S22, E11, E22, workTotal, workx, worky

    def plot(self, wlrange=[1.2e-6, 1.6e-6]):
        wls = np.linspace(wlrange[0], wlrange[1], 100)
        omega = 2 * pi * 3e8 / (wls)

        betas = np.array([self.beta(wl, [1, 2, 3]) for wl in wls])
        gammas = np.array([self.gamma(wl) for wl in wls])

        pl.subplot(311)
        pl.plot(wls / u_, betas[:, 0], "b--", label=r"$\beta_1$")
        pl.subplot(312)
        pl.plot(wls / u_, betas[:, 1], "g--", label=r"$\beta_2$")
        pl.subplot(313)
        pl.plot(wls / u_, betas[:, 2], "r--", label=r"$\beta_2$")

        # Integrate:
        D = 0.25 * self.S0 * (wls - self.lambda0 ** 4 / wls ** 3)
        b1_calc = integrate.cumtrapz(D, wls)
        b1_calc = np.append(b1_calc, b1_calc[-1])

        b_calc = integrate.cumtrapz(b1_calc, omega)
        b_calc = np.append(b_calc, b_calc[-1])
        print b1_calc + betas[0, 0]

        # betaw = betas[:,0]*omega+betas[:,1]*omega**2+betas[:,2]*omega**3
        pl.subplot(311)
        pl.plot(wls / u_, b1_calc + betas[0, 0], "r-", label=r"$\beta_1 calc$")

 def geo_amin( self, lf, phi, depsmax, Kc, Ef, Vf, damage ):
         # print self.amin_it
         a = np.linspace( 0, self.amin_it, 300 )
         Kf = Vf * self.sorted_nu_r * self.sorted_stats_weights * Ef
         a_shaped = a.reshape( len( a ), 1 )
         phi = phi.reshape( 1, len( phi ) )
         m_p = a_shaped * 2 / np.cos( phi ) / lf
         mask1 = m_p > 0
         m_p = m_p * mask1
         p = np.abs( H( 1 - m_p ) * ( 1 - m_p ) )
         muT = np.sum( self.sorted_depsf * ( 1 - damage ) * Kf * p , 1 )
         Kf_intact = np.sum( ( 1 - damage ) * Kf * p , 1 )[::-1]
         Kf_broken = np.sum( Kf * damage )
         Emtrx = ( 1. - self.V_f_tot ) * self.E_m + Kf_broken + Kf_intact
         depsm = muT / Emtrx
         # print depsm
         em = np.hstack( ( 0, cumtrapz( depsm, a ) ) )
         um = np.hstack( ( 0, cumtrapz( em , a ) ) )
         plt.plot( a, em )
         plt.plot( a, um )
         plt.show()
         ind = np.argmin( np.abs( self.w - depsmax * a ** 2. / 2. + em * a - um ) )
         amin = a[:ind + 1]
         self.amin_it = 1.2 * a[ind]
         # print amin
         return amin, em[:ind + 1]

def cumtrapzmid(x, y, c):
    """
    cumulative trapezoidal numerical integration taken from midpoint

    :param x: vector of size N describing the time samples
    :param y: vector of size N describing the function
    :param c: midpoint

    :rtype: vector
    :return fa: cumulative integration

    """
    a = x.shape[0]
    mid = int(round(a / 2.))

    # case < mid
    fa = zeros(a)
    tmpx = x[0:mid]
    tmpy = y[0:mid]
    tmp = c + cumtrapz(tmpy[::-1], tmpx[::-1], initial=0)
    fa[0:mid] = tmp[::-1]

    # case >= mid
    fa[mid:a] = c + cumtrapz(y[mid - 1:a - 1], x[mid - 1:a - 1], initial=0)

    return fa

    def calc_stream(self):
        """computes and returns the stream function

        only make sense with vp fields
        """
        # should add test if vp fields or not
        vphi = self.fields['v'][:, :, 0]
        # interpolate to the same phi
        vph2 = -0.5 * (vphi + np.roll(vphi, 1, 1))
        v_r = self.fields['w'][:, :, 0]
        n_r, nph = np.shape(v_r)
        stream = np.zeros(np.shape(vphi))
        # integrate first on phi
        stream[0, 1:nph - 1] = self.rcmb * \
            integrate.cumtrapz(v_r[0, 0:nph - 1], self._ph_coord)
        stream[0, 0] = 0
        # use r coordinates where vphi is defined
        rcoord = self.rcmb + np.array(
            self.rgeom[0:np.shape(self.rgeom)[0] - 1:2])
        for iph in range(0, np.shape(vph2)[1] - 1):
            stream[1:n_r, iph] = stream[0, iph] + \
                integrate.cumtrapz(vph2[:, iph], rcoord)  # integrate on r
        stream = stream - np.mean(stream[n_r / 2, :])
        # remove some typical value.
        # Would be better to compute the global average
        # taking into account variable grid spacing
        return stream

def properdistance(z,omegam=0.3,omegax=0.7,w0=-1,w1=0,wz=None):
    """
    Gives the proper distance in the defined cosmology
    The c/Ho factor is ommited
    Returns dist(z), w(z), omegax(z), H(z), curvature
    """
    # if no wz on input the calculate it from w0 and w1
    if wz is None: wz=w0+(z*1./(1.+z))*w1
    # calculate evolution of omegax accounting for its equation of state
    omegaxz=zeros(z.size)
    omegaxz[0]=omegax
    omegaxz[1:z.size]=omegax*exp(3*integrate.cumtrapz((1.+wz)/(1.+z),x=z))

    # curvature
    omega=omegam+omegax
    
    # calculation of H(z)
    hz=sqrt((1.-omega)*(1+z)**2+omegaxz+omegam*(1+z)**3)

    # calculate chi
    chi=zeros(z.size)
    chi[1:z.size]=integrate.cumtrapz(1./hz,x=z)

    #calculate proper distance
    if omega>1: curv=1
    if omega<1: curv=-1
    if omega==1: curv=0
    kk=abs(1.-omega)
    if curv==1: dist=sin(sqrt(kk)*chi)/sqrt(kk)
    if curv==-1: dist=sinh(sqrt(kk)*chi)/sqrt(kk)
    if curv==0: dist=chi

    return dist,wz,omegaxz,hz,curv,chi

 def geo_amin_lmin(self, damage, depsf, umLmin, emLmin, Lmin, idx):
     phi = self.sorted_phi
     Vf = self.sorted_V_f
     depsmax = depsf
     Ef = self.sorted_E_f
     al = self.sorted_l_ez
     a = np.linspace(Lmin, Lmin * 1.5, self.discr_amin)
     Kf = Vf * self.sorted_nu_r * self.sorted_stats_weights * Ef * self.soc
     a_shaped = a.reshape(len(a), 1)
     al = al.reshape(1, len(al))
     m_p = a_shaped / al
     m_p = m_p
     p = np.abs(H(1 - m_p) * (1 - m_p))
     ####
     muT = np.sum(self.sorted_depsf[idx:] * (1 - damage[idx:]) * Kf[idx:] * p[:, idx:], 1)
     Kf_intact = np.sum(((1 - damage) * Kf * (1 - p))[:, idx:], 1)
     Kf_broken = np.sum(Kf * damage)
     Emtrx = (1.0 - self.V_f_tot) * self.E_m + Kf_broken + Kf_intact
     depsm = muT / Emtrx
     emr = np.hstack((emLmin, cumtrapz(depsm, a)))
     umr = np.hstack((umLmin, cumtrapz(emr, a)))
     condI = self.w - (
         Lmin * (emr - emLmin + (a - Lmin) * depsmax)
         + depsmax * Lmin ** 2.0 / 2.0
         + emLmin * Lmin
         - umLmin
         + (depsmax * a ** 2.0 / 2.0 + emr * a - umr)
     )
     ip_f = interp1d(condI[::-1], a[::-1], bounds_error=False, fill_value=a[-1])
     cut_at = np.nonzero(H(ip_f(0) - a))
     amin = np.hstack((a[cut_at], ip_f(0)))
     return amin[-1]

 def geo_amin( self, damage, Lmax ):
         l_ez = self.sorted_l_ez
         depsmax = self.sorted_depsf[0]
         a = np.linspace( 0, np.min( ( self.amin_it[-1] * 1.2, Lmax + 1e-6 ) ) , self.discr_amin )
         Kf = self.Kf
         a_shaped = a.reshape( len( a ), 1 )
         l_ez2D = l_ez.reshape( 1, len( l_ez ) )
         m_p = a_shaped / l_ez2D
         p = np.abs( H( 1 - m_p ) * ( 1 - m_p ) )
         ####
         muT = np.sum( self.sorted_depsf * ( 1 - damage ) * Kf * p , 1 )
         Kf_intact = np.sum( ( 1 - damage ) * Kf * ( 1 - p ) , 1 ) 
         Kf_broken = np.sum( Kf * damage )
         Emtrx = ( 1. - self.V_f_tot ) * self.E_m + Kf_broken + Kf_intact
         depsm = muT / Emtrx
         em = np.hstack( ( 0, cumtrapz( depsm, a ) ) )
         um = np.hstack( ( 0, cumtrapz( em , a ) ) )
         ######
         # Aphi = pi * 0.03 ** 2 / np.cos( self.sorted_phi )
         # llambda = ( 0.03 - 0.03 / np.cos( self.sorted_phi) ) / ( 0.03 + 0.03/ np.cos( self.sorted_phi ) )
         # Uphi = pi * ( 0.03 + 0.03/ np.cos( self.sorted_phi ) ) * ( 1. + 3.*llambda ** 2. / ( 10. + ( 4 - 3. * llambda ** 2. ) ) )
         ######
         condI = self.w / 2. - ( depsmax * a ** 2. / 2. + em * a - um )
         ip_f = interp1d( condI[::-1], a[::-1], bounds_error = False, fill_value = a[-1] )
         cut_at = np.nonzero( H( ip_f( 0 ) - a ) )
         
         a_new = np.hstack( ( a[cut_at], ip_f( 0 ) ) )
         amin = a_new
         self.amin_it = a_new
         # interp depsm and em
         ind = len( amin )
         interp_data = ( a_new[-1] - a[ind - 2] ) / ( a[ind - 1 ] - a[ind - 2] )
         diff = depsm[ind - 1 ] - depsm[ind - 2]
         depsm[ind - 1] = depsm[ind - 2] + interp_data * diff
         return amin, depsm[:len( amin ) ], Emtrx[:len( amin ) ]

 def geo_amin(self, damage, Lmax):
     phi = self.sorted_phi
     Vf = self.sorted_V_f
     lf = self.sorted_lf
     al = self.sorted_l_ez
     depsmax = self.sorted_depsf[0]
     Ef = self.sorted_E_f
     a = np.linspace(0, np.min((self.amin_it[-1] * 1.2, Lmax + 1e-6)), self.discr_amin)
     Kf = Vf * self.sorted_nu_r * self.sorted_stats_weights * Ef * self.soc
     a_shaped = a.reshape(len(a), 1)
     al = al.reshape(1, len(al))
     m_p = a_shaped / al
     m_p = m_p
     p = np.abs(H(1 - m_p) * (1 - m_p))
     ####
     muT = np.sum(self.sorted_depsf * (1 - damage) * Kf * p, 1)
     Kf_intact = np.sum((1 - damage) * Kf * (1 - p), 1)
     Kf_broken = np.sum(Kf * damage)
     Emtrx = (1.0 - self.V_f_tot) * self.E_m + Kf_broken + Kf_intact
     depsm = muT / Emtrx
     em = np.hstack((0, cumtrapz(depsm, a)))
     um = np.hstack((0, cumtrapz(em, a)))
     condI = self.w / 2.0 - (depsmax * a ** 2.0 / 2.0 + em * a - um)
     ip_f = interp1d(condI[::-1], a[::-1], bounds_error=False, fill_value=a[-1])
     cut_at = np.nonzero(H(ip_f(0) - a))
     a_new = np.hstack((a[cut_at], ip_f(0)))
     amin = a_new
     self.amin_it = a_new
     # interp depsm and em
     ind = len(amin)
     interp_data = (a_new[-1] - a[ind - 2]) / (a[ind - 1] - a[ind - 2])
     diff = depsm[ind - 1] - depsm[ind - 2]
     depsm[ind - 1] = depsm[ind - 2] + interp_data * diff
     return amin, depsm[: len(amin)]

    def calculate_pressure_perturbation(self):
        """Perturbation pressure divided by density.
        Assumes hydrostatic balance.
        See: Kunze et. al. 2002 JPO.
        See: Nash et. al. 2005 JPO."""

        self.Pprime = self.P.copy()

        for i in xrange(len(self.hpid)):

            nans = np.isnan(self.P[:, i])

            z = self.z[~nans, i]
            b = self.b[~nans, i]

            # z should be increasing.
            if z[0] > z[-1]:
                z = np.flipud(z)
                b = np.flipud(b)

                bi = cumtrapz(b, z, initial=0.)
                bii = cumtrapz(bi, z, initial=0.)

                Pprime = bi + (bii[0] - bii[-1])/(-z[0])

                self.Pprime[~nans, i] = np.flipud(Pprime)

            else:
                bi = cumtrapz(b, z, initial=0.)
                bii = cumtrapz(bi, z, initial=0.)

                self.Pprime[~nans, i] = bi + (bii[0] - bii[-1])/(-z[0])

 def amin_i(self, demsi, depsf0, depsf_smaller, a1, damage, idx):
     Vf = self.sorted_V_f
     Ef = self.sorted_E_f
     al = self.sorted_l_ez[idx + 1 :]
     # print al
     a = np.linspace(a1, self._a_long[idx + 1] * 1.2, self.discr_amin)
     Kf = Vf * self.sorted_nu_r * self.sorted_stats_weights * Ef * self.soc
     a_shaped = a.reshape(len(a), 1)
     al = al.reshape(1, len(al))
     m_p = a_shaped / al
     p = np.abs(H(1 - m_p) * (1 - m_p))
     # print p
     ####
     muT = np.sum(self.sorted_depsf[idx + 1 :] * (1 - damage[idx + 1 :]) * Kf[idx + 1 :] * p, 1)
     # print muT
     Kf_intact = np.sum((1 - damage[idx + 1 :]) * Kf[idx + 1 :] * (1 - p), 1) + np.sum(
         (1 - damage[: idx + 1]) * Kf[: idx + 1]
     )
     Kf_broken = np.sum(Kf * damage)
     Emtrx = (1.0 - self.V_f_tot) * self.E_m + Kf_broken + Kf_intact
     depsm = muT / Emtrx
     em = np.hstack((0, cumtrapz(depsm, a)))
     um = np.hstack((0, cumtrapz(em, a)))
     condI = a ** 2.0 / 2.0 * depsf0 + a * em - um - a1 ** 2 / 2.0 * depsf_smaller
     ip_f = interp1d(condI, a, bounds_error=False, fill_value=a[-1])
     cut_at = np.nonzero(H(ip_f(0) - a))
     amin = np.hstack((a[cut_at], ip_f(0)))
     # print self.sorted_depsf[self.sf_mask][0]
     # print self.sorted_depsf[self.c_mask][0]
     return amin[-1]

def compute_jdot_grav_profile(snapshot,rad_list,code="AREPO", alpha = 0.1, h0 = 0.1):
    
    X0, Y0 = 0.5 * snapshot.header.boxsize, 0.5 * snapshot.header.boxsize
    vol = snapshot.gas.MASS/snapshot.gas.RHO
    # add gradients
    jdotdens = -snapshot.gas.RHO * ((snapshot.gas.POS[:,0] - X0) * (-snapshot.gas.ACCE[:,1]) - \
                                    (snapshot.gas.POS[:,1] - Y0) * (-snapshot.gas.ACCE[:,0])) 
    ind = snapshot.gas.ID > -2
    jdotdens_av=compute_profiles(jdotdens[ind],snapshot.gas.R[ind],vol[ind],rad_list)
    jdot_av = 2 * np.pi *  cumtrapz((jdotdens_av[:] * rad_list[:])[::-1],x = -rad_list[::-1],initial=0)[::-1]
    #jdot_av = 2 * np.pi *  cumtrapz((jdotdens_av[:] * rad_list[:]),x = rad_list,initial=0)
    
    return np.array(jdot_av)



    
    # Compute the cell-centered quantities
    jdotdens_per_cell = -snapshot.gas.RHO * ((snapshot.gas.POS[:,0] - X0) * (-snapshot.gas.ACCE[:,1]) - \
                                             (snapshot.gas.POS[:,1] - Y0) * (-snapshot.gas.ACCE[:,0])) 
    snapshot.add_data(jdotdens_per_cell,'TORQUEDENS')
    # interpolate onto the grid
    jdotdens_interp = disk_interpolate_primitive_quantities(snapshot,[gridX,gridY],\
                                                            quantities=['TORQUEDENS'],method = 'nearest')[0]

    del snapshot.gas.TORQUEDENS
    
    # In the case of gravity, we need to carry out an additional integration step
    gridR = grid.R.mean(axis=0)
    jdot_interp = cumtrapz((jdotdens_interp * grid.R)[:,::-1],x = -gridR[::-1],initial=0,axis=1)[:,::-1] / grid.R

    
    return jdot_interp

    def test_1d(self):
        x = np.linspace(-2, 2, num=5)
        y = x
        y_int = cumtrapz(y, x, initial=0)
        y_expected = [0., -1.5, -2., -1.5, 0.]
        assert_allclose(y_int, y_expected)

        y_int = cumtrapz(y, x, initial=None)
        assert_allclose(y_int, y_expected[1:])

def numerical_potential(r,dens):
    '''Integrates the density profile to solve for the potential profile. This is the 2 integral method.
    Returned units are in Mpc^2/s^2
    '''
    deriv1 = dens*r**2
    deriv2 = dens*r
    inner = cumtrapz(deriv1,r)
    outer = -cumtrapz(deriv2[::-1],r[::-1])
    return -4*np.pi*G*(1.0/r[1:-1]*inner[:-1] + outer[::-1][1:])

    def __add__(self,fld):
        newfld = copy.deepcopy(self)
        newfld.vx += fld.vx
        newfld.vy += fld.vy
        newfld.vz += fld.vz
        newfld.vx2 += fld.vx2
        newfld.vy2 += fld.vy2
        newfld.vz2 += fld.vz2
        newfld.KE += fld.KE
        newfld.dTdz += fld.dTdz
        newfld.T += fld.T
        newfld.T2 += fld.T2
        newfld.Fk += fld.Fk
        newfld.Fc += fld.Fc
        newfld.Fe += fld.Fe
        newfld.kap += fld.kap
        newfld.vort += fld.vort
        newfld.area += fld.area
        newfld.vortx += fld.vortx
        newfld.vorty += fld.vorty
        newfld.visc3 += fld.visc3
        newfld.therm += fld.therm
        newfld.P += fld.P
        newfld.dzP += fld.dzP
        newfld.Pflux += fld.Pflux
        newfld.TdzP += fld.TdzP
        newfld.trip1 += fld.trip1
        newfld.trip2 += fld.trip2
        newfld.rey12 += fld.rey12
        newfld.rey13 += fld.rey13
        newfld.rey23 += fld.rey23
        newfld.vflux += fld.vflux
        newfld.keflux += fld.keflux
        newfld.tflux += fld.tflux
        newfld.thermflux += fld.thermflux
        newfld.Fcflux1 += fld.Fcflux1
        newfld.Fcflux2 += fld.Fcflux2
        newfld.visc1 += fld.visc1
        newfld.visc2 += fld.visc2



        newfld.thermflux *= newfld.kap[:,-1][:,np.newaxis]
        newfld.Fb = -newfld.vz*newfld.delad*newfld.z[:,np.newaxis]
        newfld.Nu = newfld.Fc/newfld.Fk

        newfld.Ftot = newfld.Fk + newfld.Fc
        newfld.vrms = np.sqrt(newfld.vz2)

        newfld.s = newfld.T - newfld.delad*newfld.z[:,np.newaxis]
        newfld.dsdz = newfld.dTdz - newfld.delad
        newfld.sk = cumtrapz(-1./newfld.kap[:,-1]-newfld.delad,x=newfld.z,initial=0)
        tmp = cumtrapz((-1./newfld.kap[:,-1]-newfld.delad)[::-1],x=newfld.z[::-1],initial=0)[::-1]
        newfld.sk_t = newfld.s[-1,:] + tmp[:,np.newaxis]
        return newfld