def MEVG0(Rhwhm, RefE, RefV, Rrms) :
	print '... ... ... MEVG0: Rhwhm, RefE, RefV, Rrms', Rhwhm, RefE, RefV, Rrms

	if RefE !=0 and RefV !=0 :
		print "Ooops RE OR RV please"
		sys.exit(0)

	if RefE == 0 and RefV == 0:
		result = MG0(Rhwhm, Rrms)

	if RefV == 0 and RefE != 0 :
		result = tplquad(MEVG0_integrand_exp,\
						0.0, numpy.pi, \
						lambda th: 0.0,\
						lambda th: 3.0*Rrms,\
						lambda th,rt: 0.0,\
						lambda th,rt: 5.0*RefE,\
						args = (Rhwhm, RefE, RefV, Rrms),\
						)[0]

	if RefE == 0 and RefV != 0:
		result = tplquad(MEVG0_integrand_vau,\
						0, numpy.pi,
						lambda th: 0.0,\
						lambda th: 3.0*Rrms,\
						lambda th,rt: 0.0,\
						lambda th,rt: 5.0*RefV,\
						args = (Rhwhm, RefE, RefV, Rrms),\
						)[0]


	return result

def losvprojden(param = [2.0, -5.3, 2.5, 0.16, 1.5, -9.0, 6.9, 0.086, 21.0, 1.5, 1,3,1], dR = 0.05):
	a,d,e, Ec, rlim, b, q, Jb, Vmax, rmax, alpha, beta, gamma = param
	rs = rmax/2.16           # rmax=2.16*rs

	losvsiglist = []
        projdenlist = []
        Rlist = []
        Ri = 0.0000
        while Ri<rlim:  #for each Ri, integrate over r
		def losvI(theta, v,r):
			zz = r*r - Ri*Ri + 10**-10  #add 10**-10 to prevent divide by zero
                	dz = r / np.sqrt(zz)
			f  = ftheta( (r/rs), v, theta, param )
			Ja = f * zz*pow(v,4)*np.cos(theta)*np.cos(theta)*np.sin(theta)/(r*r)
			Jb = f * Ri*Ri*pow(v,4)*pow(np.sin(theta),3)/(r*r)
			print 'using Ri', Ri, theta, v, r
			#missing rhor again! WRONG
			rhor = ftheta( (r/rs), v, theta, param ) * v*v * np.sin(theta) 
			return (Ja+Jb)*dz

		def projI(theta,v,r):
			zz = r*r - Ri*Ri + 10**-10
                        dz = r / np.sqrt(zz)
			I  = ftheta( (r/rs), v, theta, param ) * v*v * np.sin(theta) * dz
			return I
	
		def bounds_theta(v,r):
			return [0, np.pi]
		def bounds_v(r):
			vmax = vesc( (r/rs), [rlim, Vmax, rmax, alpha, beta, gamma] )
	 		return [0, vmax]
		def bounds_r():
			return [Ri, rlim]
		
		def vmax(r):
                        return vesc( (r/rs), [rlim, Vmax, rmax, alpha, beta, gamma] )

		def vmin(r):
                        return 0.0
		def tmax(v,r):
			return np.pi
		def tmin(v,r):
			return 0.0

		#Ilos  = integrate.nquad( losvI, [bounds_theta, bounds_v, bounds_r], epsrel = 1e-6 )[0] 
		#Iproj = integrate.nquad( projI, [bounds_theta, bounds_v, bounds_r], epsrel = 1e-6 )[0]
		
		Ilos = integrate.tplquad(losvI, Ri, rlim, vmin,vmax, tmin,tmax, epsrel=1e-4, epsabs=0)[0]
		Iproj= integrate.tplquad(projI, Ri, rlim, vmin,vmax, tmin,tmax, epsrel=1e-4, epsabs=0)[0]

		losvsiglist.append( Ilos  )
		projdenlist.append( Iproj )
		Rlist.append(Ri)
		print "3quad ftheta: ", Ri, rlim
		Ri = Ri + dR

	Rarr	    = np.asarray(Rlist)
	projdenarr  = np.asarray(projdenlist)
        losvsigarr2 = np.asarray(losvsiglist) / (projdenarr+10**-8)
        return Rarr, losvsigarr2**0.5, projdenarr/(max(projdenarr) + 10**-8)

 def calc_ff(self,Q_VAL):
     for q in Q_VAL:
         r = self.r	
         a_int= integrate.tplquad(lambda x,y,z:np.cos(q*(x**2+y**2+z**2)**.5), 
                 -r, r, 
                 lambda x:-(r**2-x**2)**.5, lambda x: (r**2-x**2)**.5, 
                 lambda x,y: -(r**2-x**2-y**2)**.5, lambda x,y: (r**2-x**2-y**2)**.5)
         b_int = integrate.tplquad(lambda x,y,z:np.sin(q*(x**2+y**2+z**2)**.5), 
                 -r, r, 
                 lambda x:-(r**2-x**2)**.5 ,lambda x: (r**2-x**2)**.5, 
                 lambda x,y: -(r**2-x**2-y**2)**.5, lambda x,y: (r**2-x**2-y**2)**.5)	
         pq = (a_int[0]-1j*b_int[0])
         print(pq)	
         yield pq

def integrate_nd(f, domain, shape, dtype):

    if shape == () or shape == (1,):
        if dtype in continuous_types:
            return integrate.quad(f, domain.lower, domain.upper, epsabs=1e-8)[0]
        else:
            return np.sum(list(map(f, np.arange(domain.lower, domain.upper + 1))))
    elif shape == (2,):
        def f2(a, b):
            return f([a, b])

        return integrate.dblquad(f2,
                                 domain.lower[0],
                                 domain.upper[0],
                                 lambda a: domain.lower[1],
                                 lambda a: domain.upper[1])[0]

    elif shape == (3,):
        def f3(a, b, c):
            return f([a, b, c])

        return integrate.tplquad(f3,
                                 domain.lower[0], domain.upper[0],
                                 lambda a: domain.lower[1],
                                 lambda a: domain.upper[1],
                                 lambda a, b: domain.lower[2],
                                 lambda a, b: domain.upper[2])[0]
    else:
        raise ValueError("Dont know how to integrate shape: " + str(shape))

def integrate_that_shit():
  volume = tplquad(diff_volume, r1, r3, lambda r:   t1, lambda r:   t2,
                                        lambda r,t: p1, lambda r,t: p2)[0]
  print 'volume check (' + str(4./3.*np.pi*(r3)**3.) + '): ' + str(volume)

  print 'now the nasty integral...'
  Uabs_T = tplquad(new_diff_Uabs, r1, r3, lambda r:   t1, lambda r:   t2,
                                        lambda r,t: p1, lambda r,t: p2)[0]

  print '... no problem!'

  eps0 = 8.854187E-12
  mu0 = 4.*np.pi*1.E-7
  Ap = np.pi*radshell*radshell
  Qabs_corrected = Uabs_T/Ap/(1./2.*np.sqrt(eps0/mu0))
  print 'my calc divided by sqrt(eps0/mu0) = ' + str(Qabs_corrected*1.E-6)

 def test_norm(self):
     """
     Tests whether distribution function is normalized, and integrates to 1.
     """
     # converting vTh to unitless
     vTh = self.vTh.si.value
     # setting up integration from 0 to 10*vTh, which is close to Inf
     infApprox = (10 * vTh)
     # integral should be close to 1
     integ = spint.tplquad(Maxwellian_speed_3D,
                           0,
                           infApprox,
                           lambda z: 0,
                           lambda z: infApprox,
                           lambda z, y: 0,
                           lambda z, y: infApprox,
                           args=(self.T,
                                 self.particle,
                                 0,
                                 0,
                                 0,
                                 vTh,
                                 "unitless"),
                           epsabs=1e0,
                           epsrel=1e0,
                           )
     # value returned from tplquad is (integral, error), we just need
     # the 1st
     integVal = integ[0]
     exceptStr = ("Integral of distribution function should be 1 "
                  f"and not {integVal}")
     assert np.isclose(integVal,
                       1,
                       rtol=1e-3,
                       atol=0.0), exceptStr

def pcs_numeric_quad_int_iso_cone(theta_max=None, sigma_max=None, c=None, r_pivot_atom=None, r_ln_pivot=None, A=None, R_eigen=None, RT_eigen=None, Ri_prime=None):
    """Determine the averaged PCS value via numerical integration.

    @keyword theta_max:     The half cone angle.
    @type theta_max:        float
    @keyword sigma_max:     The maximum torsion angle.
    @type sigma_max:        float
    @keyword c:             The PCS constant (without the interatomic distance and in Angstrom units).
    @type c:                float
    @keyword r_pivot_atom:  The pivot point to atom vector.
    @type r_pivot_atom:     numpy rank-1, 3D array
    @keyword r_ln_pivot:    The lanthanide position to pivot point vector.
    @type r_ln_pivot:       numpy rank-1, 3D array
    @keyword A:             The full alignment tensor of the non-moving domain.
    @type A:                numpy rank-2, 3D array
    @keyword R_eigen:       The eigenframe rotation matrix.
    @type R_eigen:          numpy rank-2, 3D array
    @keyword RT_eigen:      The transpose of the eigenframe rotation matrix (for faster calculations).
    @type RT_eigen:         numpy rank-2, 3D array
    @keyword Ri_prime:      The empty rotation matrix for the in-frame isotropic cone motion, used to calculate the PCS for each state i in the numerical integration.
    @type Ri_prime:         numpy rank-2, 3D array
    @return:                The averaged PCS value.
    @rtype:                 float
    """

    # Perform numerical integration.
    result = tplquad(pcs_pivot_motion_full_quad_int, -sigma_max, sigma_max, lambda phi: -pi, lambda phi: pi, lambda theta, phi: 0.0, lambda theta, phi: theta_max, args=(r_pivot_atom, r_ln_pivot, A, R_eigen, RT_eigen, Ri_prime))

    # The surface area normalisation factor.
    SA = 4.0 * pi * sigma_max * (1.0 - cos(theta_max))

    # Return the value.
    return c * result[0] / SA

 def __init__(self, r0, basis, mean, covar):
     self.r0 = r0
     self.basis = basis
     self.mean = mean
     self.covar = covar
     self.C = 1
     self.C = integrate.tplquad(lambda psi, phi, theta: sin(psi)**2*sin(phi)*self.pdfSphere(Quaternion(cos(psi), sin(psi)*cos(phi), sin(psi)*sin(phi)*cos(theta), sin(psi)*sin(phi)*sin(theta))), 0, 2*pi, lambda theta: 0, lambda theta: pi, lambda phi, theta: 0, lambda phi,theta: pi)[0]

def _scipy_integrate(f,lb,ub):  
  """
Returns the integral of an n-dimensional function f from lb to ub
(where n = 1, 2, or 3), using scipy.integrate

Inputs:
    f -- a function that takes a list and returns a number.
    lb -- a list of lower bounds
    ub -- a list of upper bounds
"""
  from scipy.integrate import quad, dblquad, tplquad
  if len(lb) == 3:
    def func(z,y,x): return f([x,y,z])
    def qf(x,y): return lb[2]
    def rf(x,y): return ub[2]
    def gf(x): return lb[1]
    def hf(x): return ub[1]
    expectation,confidence = tplquad(func,lb[0],ub[0],gf,hf,qf,rf)
    return expectation
  if len(lb) == 2:
    def func(y,x): return f([x,y])
    def gf(x): return lb[1]
    def hf(x): return ub[1]
    expectation,confidence = dblquad(func,lb[0],ub[0],gf,hf)
    return expectation 
  if len(lb) == 1:
    expectation,confidence = quad(f,lb[0],ub[0])
    return expectation 

    def spherical_vol_avg(self,center,r_max):
        """
        Averages the density over the volume of a sphere centered at the point
        center, with a radis r_max.
        """
        # Physics notation for spherical coordinates are use:
        # rho - radial distance
        # theta - inclination
        # phi - azimuthal angle
        from scipy.integrate import tplquad
        def func(theta,phi,rho,center,dens):
            """
            Function to find the value of the density at a point equal to
            center+Vector(rho).
            """
            # Determine direct rectilinear coordinates relative to center from
            # spherical coordinates away from center.
            rho_x = rho*np.cos(phi)*np.sin(theta)
            rho_y = rho*np.sin(phi)*np.sin(theta)
            rho_z = rho*np.cos(theta)
            # Determine absolute direct coordinates from relative coordinates
            # and center.
            r_x = rho_x + center[0]
            r_y = rho_y + center[1]
            r_z = rho_z + center[2]
            # Find the value of the density at the point r.
            val = dens.interpolate([r_x,r_y,r_z])
            # Weight the value by the spherical coordinate Jacobian.
            val = val*rho*rho*np.sin(theta)
            return val

        integral,error = tplquad(func,0,r_max,lambda a:0,lambda b:2*np.pi,
                                 lambda c,d:0,lambda e,f: np.pi,args=(center,self))
        return integral

    def test_matching_tplquad(self):
        def func3d(x0, x1, x2, c0, c1):
            return x0 ** 2 + c0 * x1 ** 3 - x0 * x1 + 1 + c1 * np.sin(x2)

        res = tplquad(func3d, -1, 2, lambda x: -2, lambda x: 2, lambda x, y: -np.pi, lambda x, y: np.pi, args=(2, 3))
        res2 = nquad(func3d, [[-np.pi, np.pi], [-2, 2], (-1, 2)], args=(2, 3))
        assert_almost_equal(res, res2)

def heroes_effective_area_gaussian(energy_range=(20,30), fwhm=3, radius=9.5,
                                   offaxis=0):
    """
    Calculates the average effective area for a Gaussian exposure.  Off-axis
    exposures take *much* longer to calculate, by >~2 orders of magnitude in
    time!
    
    energy_range is in keV
    fwhm of source is in arcmin
    radius of integration area is in arcmin
    offaxis is in arcmin
    """
    f2d = heroes_effective_area_fit()
    sigma = fwhm/2.355
    if offaxis == 0:
        area = integrate.dblquad(lambda e,r: f2d(e,r)*r*np.exp(-(r/sigma)**2/2)/sigma**2,
                                 0, radius,
                                 lambda r:energy_range[0], lambda r: energy_range[1])[0]
    else:
        r2 = lambda x,y: x**2+y**2
#        area = integrate.dblquad(lambda y,x: np.exp(-(r(x-offaxis,y)/sigma)**2/2)/(2*np.pi*sigma**2),
#                                 -radius, radius,
#                                 lambda x: -np.sqrt(radius**2-x**2), lambda x: np.sqrt(radius**2-x**2))[0]
        area = integrate.tplquad(lambda e,y,x: f2d(e,np.sqrt(r2(x,y)))*np.exp(-r2(x-offaxis,y)/sigma**2/2)/(2*np.pi*sigma**2),
                                 -radius, radius,
                                 lambda x: -np.sqrt(radius**2-x**2), lambda x: np.sqrt(radius**2-x**2),
                                 lambda x,y: energy_range[0], lambda x,y: energy_range[1])[0]
    norm_area = 1.
    area /= norm_area*(energy_range[1]-energy_range[0])
    return area

def boxIntegrate(cx, cy, r, zMin, zMax, f):
    yBot = lambda x : cy - r
    yTop = lambda x : cy + r
    zBot = lambda x, y : zMin
    zTop = lambda x, y : zMax
    
    res, err = integrate.tplquad(f, cx - r, cx + r, yBot, yTop, zBot, zTop)
    return (res,err)

def cylinderIntegrate(cx, cy, r, zMin, zMax, f):
    yTop = lambda x : cy + (r*r - (x - cx)*(x - cx))**0.5
    yBot = lambda x : cy - (r*r - (x - cx)*(x - cx))**0.5
    zTop = lambda x, y : zMax
    zBot = lambda x, y : zMin
    
    res, err = integrate.tplquad(f, cx - r, cx + r, yBot, yTop, zBot, zTop)
    return res, err

    def test_triple_integral(self):
        # 9) Triple Integral test
        def simpfunc(z,y,x):      # Note order of arguments.
            return x+y+z

        a, b = 1.0, 2.0
        assert_quad(tplquad(simpfunc,a,b,
                            lambda x: x, lambda x: 2*x,
                            lambda x,y: x-y, lambda x,y: x+y),
                    8/3.0 * (b**4.0 - a**4.0))

def tripleRectIntegrate(bounds, f):
    a1 = bounds[0][0]
    b1 = bounds[0][1]
    a2 = lambda x : bounds[1][0]
    b2 = lambda x : bounds[1][1]
    a3 = lambda x, y : bounds[2][0]
    b3 = lambda x, y : bounds[2][1]
    
    res, err = integrate.tplquad(f, a1,b1,a2,b2,a3,b3)
    return (res,err)

def b_pot_3d_cont(x, R, h, sigma, basis_func):
    """
    Returns the value of the potential at point (x,y,0) generated
    by a basis source located at (0,0,0)
    """
    pot, err = integrate.tplquad(int_pot_3D, -R, R,
                                 lambda x: -R, lambda x: R,
                                 lambda x, y: -R, lambda x, y: R,
                                 args=(x, R, h, basis_func))
    pot *= 1./(2.0*pi*sigma)  # TODO check the constant
    return pot

def heroes_atmospheric_attenuation(energy_range = (20, 30),
                                   altitude = (40, 40.1),
                                   elevation = (45, 45.1)):
    """First attempt, not fully tested and *very* slow"""
    f = lambda en, al, el: xray_transmission_in_atmosphere(en, al, view_angle=el)
    attenuation = integrate.tplquad(f, elevation[0], elevation[1],
                                    lambda el: altitude[0], lambda el: altitude[1],
                                    lambda el, al: energy_range[0], lambda el, al: energy_range[1])[0]
    attenuation /= energy_range[1]-energy_range[0]
    attenuation /= altitude[1]-altitude[0]
    attenuation /= elevation[1]-elevation[0]

    return attenuation

 def forward_model(self, x, R, h, sigma, src_type):
     """FWD model functions
     Evaluates potential at point (x,0) by a basis source located at (0,0)
     Utlizies sk monaco monte carlo method if available, otherwise defaults
     to scipy integrate
     Parameters
     ----------
     x : float
     R : float
     h : float
     sigma : float
     src_type : basis_3D.key
     Returns
     -------
     pot : float
         value of potential at specified distance from the source
     """
     if src_type.__name__ == "gauss_3D":
         if x == 0: x=0.0001
         pot = special.erf(x/(np.sqrt(2)*R/3.0)) / x
     elif src_type.__name__ == "gauss_lim_3D":
         if x == 0: x=0.0001
         d = R/3.
         if x < R:
             e = np.exp(-(x / (np.sqrt(2)*d))**2)
             erf = special.erf(x / (np.sqrt(2)*d))
             pot = 4*np.pi * ((d**2)*(e - np.exp(-4.5)) +
                              (1/x)*((np.sqrt(np.pi/2)*(d**3)*erf) - x*(d**2)*e))
         else:
             pot = 15.28828*(d)**3 / x
         pot /= (np.sqrt(2*np.pi)*d)**3
     elif src_type.__name__ == "step_3D":
         Q = 4.*np.pi*(R**3)/3.
         if x < R:
             pot = (Q * (3 - (x/R)**2)) / (2.*R)
         else:
             pot = Q / x
         pot *= 3/(4*np.pi*R**3)
     else:
         pot, err = integrate.tplquad(self.int_pot_3D,
                                      -R,
                                      R,
                                      lambda x: -R,
                                      lambda x: R,
                                      lambda x, y: -R,
                                      lambda x, y: R,
                                      args=(x, R, h, src_type))
     pot *= 1./(4.0*np.pi*sigma)
     return pot

def cyl_integral(function, cyl_radius, height, args):
    """
    Try to do a cylindrical integral of a function
    """
    from scipy.integrate import tplquad

    # z bounds
    qfun = lambda y,z: -height #0
    rfun = lambda y,z:  height #2*np.pi

    # y bounds
    gfun = lambda y: 0
    hfun = lambda y: 2*np.pi

    # x bounds
    llim_a = 0
    ulim_b = cyl_radius

    return tplquad(function, a=llim_a, b=ulim_b, gfun=gfun, hfun=hfun,
                   qfun=qfun, rfun=rfun, args=args)