def funcbG(E,verbose = False):
     tolerance = 1.49e-8
     try:
         t = E.shape
         Gans = []
         problems = []
         for i in range(len(E)):
             print i+1, 'of', len(E)
             rapoval = rapo(E[i])
             try:
                 temp = intg.dblquad(bGinterior,0,rapoval,lambda r: 1e-4, lambda r: 1,args = (E[i],),epsabs = tolerance,epsrel = tolerance)
             except UserWarning as e:
                 if verbose == True:
                     print 'G, E = ', E[i], 'message = ', e
                     problems.append(i)
             Gans.append(temp[0])
         return array(Gans),problems
     except AttributeError:
         rapoval = rapo(E)
         problem = []
         try:
             temp = intg.dblquad(bGinterior,0,rapoval,lambda r: 0, lambda r: 1,args = (E,verbose))
         except UserWarning as e:
             if verbose == True:
                 print 'G, E = ', E, 'message = ', temp[3]
             problem = [E]
         return temp[0],problem

def test2():
    D = [[0,2], [2,4], [4,6], [6,8]]

    f = lambda x,y: np.exp(-x)*np.exp(-y)*(y-x)

    I1 = dblquad(f, 0, 8, lambda l:0, lambda l:l)[0]

    I3 = 0
    for i in np.arange(1,4):
        for j in np.arange(0,i):
            I3 += dblquad(f, D[i][0], D[i][1], lambda l:D[j][0], lambda l:D[j][1])[0]

    for i in range(4):
        I3 += dblquad(f, D[i][0], D[i][1], lambda l:D[i][0], lambda l:l)[0]

    I2 = 0
    for i in np.arange(1,4):
        for j in np.arange(0,i):
            tmp1 = gp.gintegral_seg(1,1,D[i][0],D[i][1])*gp.gintegral_seg(0,1,D[j][0],D[j][1])
            tmp2 = gp.gintegral_seg(0,1,D[i][0],D[i][1])*gp.gintegral_seg(1,1,D[j][0],D[j][1])
            I2 += tmp1-tmp2

    for i in range(4):
        I2 += dblquad(f, D[i][0], D[i][1], lambda l:D[i][0], lambda l:l)[0]

    print I1
    print I2
    print I3

    return

def over1(i,j,B,B_func,w_vec):
    """
    Calculates the first overlap integral. If it is found that it is with itself then the
    inverse effective area is returned otherwise the integrals are calculated. For the mode calculations
    the hermit-gaussian approximation is taken.
    Also the calculation is done in terms of microns^2 and is transformed in to m^2 in calc_overlaps
    Inputs::
        i,j (int,int): Integer on what whave the overlap is calculated for
        B(str vec shape[4]): Holding the mode for each wave. (lp01 or lp11)
        B_func( function vec shape[4]) : Points to what function is used to calculate each mode(field0 or field1)
        w_vec(float vec shape[2]) : The width of the lp01 or the lp11 modes. (calculated in other script)
    Local::
        fieldi,fieldj (function): Holds the ith and jth wave mode function calculator
        r(float): The radius of the fibre (there is no need to calculate infinities as the definition might give you)
        int1,int2,int3,top bottom (float vectors [4,4]): Integrals (look at Agrawal for the integrals themselves)
    Returns::
        The first overlap integrals
    """
    if i == j:

        if  B[i] == 'LP01':
            return 1/161
        elif B[i] == 'LP11':
            return 1/170
    r = 62.45
    fieldi = B_func[i]
    fieldj = B_func[j]
    int1 = lambda y,x : np.abs(fieldi(y,x,w_vec))**2 * np.abs(fieldj(y,x,w_vec))**2
    top = dblquad(int1,-r,r,lambda x : -r,lambda x: r)[0]

    int2 = lambda y,x : np.abs(fieldi(y,x,w_vec))**2
    int3 = lambda y,x : np.abs(fieldj(y,x,w_vec))**2
    bottom = dblquad(int2,-r,r,lambda x : -r,lambda x: r)[0]*\
            dblquad(int3,-r,r,lambda x : -r,lambda x: r)[0]
    return top/bottom

def funcbG(E,verbose,prereqs):
    """
    functional form of mathcalG
    relies on Ginterior
    returns mathcalG(E)
    """
    model,psigood,ggood = prereqs
    tolerance = 1.49e-8
    try:
        t = E.shape
        Gans = []
        problems = []
        for i in range(len(E)):
            print i+1, 'of', len(E)
            rapoval = rapo(E[i],psigood)
            try:
                temp = intg.dblquad(bGinterior,0,rapoval,lambda r: 1e-4, lambda r: 1,args = (E[i],prereqs),epsabs = tolerance,epsrel = tolerance)
            except UserWarning as e:
                if verbose == True:
                    print 'G, E = ', E[i], 'message = ', e
                problems.append(i)
            Gans.append(temp[0])
        return array(Gans),problems
    except AttributeError:
        rapoval = rapo(E,psigood)
        problem = []
        try:
            temp = intg.dblquad(bGinterior,0,rapoval,lambda r: 0, lambda r: 1,args = (E,prereqs))
        except UserWarning as e:
            if verbose == True:
                print 'G, E = ', E, 'message = ', temp[3]
            problem = [E]
        return temp[0],problem

    def H(self):
        """Calculates the integrated beam just over the horizon (pi/2 - 10*pi/180 <theta<pi/2). Assumes Peak power is directly overhead"""

        [O1,err1] = integrate.dblquad(self.beam_pattern_integrand,0,2*np.pi,self.H_gfun,self.H_hfun,([0,0,1],0))
        [O2,err2] = integrate.dblquad(self.beam_pattern_integrand,0,2*np.pi,self.H_gfun,self.H_hfun,([0,1,0],0))

        Peak_power = self.beam_pattern(0,0,[0,0,1]) + self.beam_pattern(0,0,[0,1,0])

        return (O1+O2)/Peak_power

def bunchlength(bunch, cavity, sigma_dz):

    print 'Iterative evaluation of bunch length...'

    counter = 0
    eps = 1

    R = cavity.circumference / (2 * np.pi)
    eta = cavity.eta(bunch)
    Qs = cavity.Qs(bunch)

    zmax = np.pi * R / cavity.h
    Hmax = cavity.hamiltonian(zmax, 0, bunch)

    # Initial values
    z0 = sigma_dz
    p0 = z0 * Qs / eta / R            #Matching condition
    H0 = eta * bunch.beta * c * p0 ** 2

    z1 = z0
    while abs(eps)>1e-6:
        # cf1 = 2 * Qs ** 2 / (eta * h) ** 2
        # dplim = lambda dz: np.sqrt(cf1 * (1 + np.cos(h / R * dz) + (h / R * dz - np.pi) * np.sin(cavity.phi_s)))
        # dplim = lambda dz: np.sqrt(2) * Qs / (eta * h) * np.sqrt(np.cos(h / R * dz) - np.cos(h / R * zmax))
        # Stationary distribution
        # psi = lambda dz, dp: np.exp(cavity.hamiltonian(dz, dp, bunch) / H0) - np.exp(Hmax / H0)

        # zs = zmax / 2.

        psi = stationary_exponential(cavity.hamiltonian, Hmax, H0, bunch)
        dplim = cavity.separatrix.__get__(cavity)
        N = dblquad(lambda dp, dz: psi(dz, dp), -zmax, zmax,
                    lambda dz: -dplim(dz, bunch), lambda dz: dplim(dz, bunch))
        I = dblquad(lambda dp, dz: dz ** 2 * psi(dz, dp), -zmax, zmax,
                    lambda dz: -dplim(dz, bunch), lambda dz: dplim(dz, bunch))

        # Second moment
        z2 = np.sqrt(I[0] / N[0])
        eps = z2 - z0

        # print z1, z2, eps
        z1 -= eps

        p0 = z1 * Qs / eta / R
        H0 = eta * bunch.beta * c * p0 ** 2

        counter += 1
        if counter > 100:
            print "\n*** WARNING: too many interation steps! There are several possible reasons for that:"
            print "1. Is the Hamiltonian correct?"
            print "2. Is the stationary distribution function convex around zero?"
            print "3. Is the bunch too long to fit into the bucket?"
            print "4. Is this algorithm not qualified?"
            print "Aborting..."
            sys.exit(-1)

    return z1

def CavityIntegral(a,b,t):
    
    # Integrate over the top and bottom half of the hexagonal chunk we are
    # trying to remove. 
    
    TopIntegral = dblquad(integrand, b, b + (np.sqrt(3.0)/2.0)*t, lambda y: y/np.sqrt(3.0) - b/np.sqrt(3.0) + a - t, lambda y: -y/np.sqrt(3.0) + b/np.sqrt(3.0) + a +t)
    BottomIntegral = dblquad(integrand,  b - (np.sqrt(3.0)/2.0)*t, b, lambda y: -y/np.sqrt(3.0) + b/np.sqrt(3.0) + a - t, lambda y: y/np.sqrt(3.0) - b/np.sqrt(3.0) + a +t)
    
    return TopIntegral[0] + BottomIntegral[0]

 def effective_area(self,lim):
     """
     Computes the effective area of mode
     """
     integrand1 = dblquad(self.Eabs2, -lim, lim, lambda x: -lim,lambda x: lim)
     integrand2 = dblquad(lambda y,x: self.Eabs2(y,x)**2, -lim, lim, lambda x: -lim,lambda x: lim)
     
     self.Aeff =  integrand1[0]**2/integrand2[0]
     return None

def prob_decay_KsKs(t1a,t1e,t2a,t2e):
	def f(t1,t2):
		return exp(-Gamma_S*t1-Gamma_S*t2)
	Int = dblquad(f, 0, np.inf, # limits t2
                  lambda x : 0, # limits t1
                  lambda x: np.inf)[0]
	Part = dblquad(f, t1a, t1e, # limits t2
                  lambda x : t2a, # limits t1
                  lambda x: t2e)[0]
	return (BrKs)**2*Part/Int

def integrate_circle(func, R, args = None):
    if args is None:
        result, _ = dblquad(func, -R, R,
                            lambda x: -np.sqrt(R**2 - x**2),
                            lambda x: np.sqrt(R**2 - x**2))
    else:
        result, _ = dblquad(func, -R, R,
                            lambda x: -np.sqrt(R**2 - x**2),
                            lambda x: np.sqrt(R**2 - x**2), args = args)
    return result

 def calc_ff(self,Q_VAL):
     for q in Q_VAL:
         muz = np.cos(self.theta)
         mul = np.cos(self.phi)
         fz = integrate.quad(lambda x: cos(q*muz*x)*self.rho(x), -self.L/2, self.L/2)
         fzi = integrate.quad(lambda x: sin(q*muz*x)*rho_L(x),-self.L/2, self.L/2)
         fl = integrate.dblquad(lambda x,y: x*cos(q*(1-mul**2)**.5*cos(y)*x)*self.rho_theta(y)*self.rho_R(R), 0, 2*3.14, 0, self.R)	
         fli = integrate.dblquad(lambda x,y: x*sin(q*(1-mul**2)**.5*cos(y)*x), 0,2*3.14, 0, self.R)
         pq = (fl[0]-1j*fli[0])*(fz[0]-1j*flz[0])*(1/(self.L*3.14*self.R**2))
         yield pq		 

    def Omega(self):
        """Calculates the integrated beam (0<theta<pi/2). Assumes Peak power is directly overhead"""

        [O1,err1] = integrate.dblquad(self.beam_pattern_integrand,0,2*np.pi,self.Omega_gfun,self.Omega_hfun,([0,0,1],0))
        [O2,err2] = integrate.dblquad(self.beam_pattern_integrand,0,2*np.pi,self.Omega_gfun,self.Omega_hfun,([0,1,0],0))

        Peak_power = self.beam_pattern(0,0,[0,0,1]) + self.beam_pattern(0,0,[0,1,0])
    #        print O1,err1
    #        print O2,err2
    #        print Peak_power
        return (O1+O2)/Peak_power

def effective(omega,xmin,xmax,ymin,ymax,l,m):
    omega1,omega2 = omega    
    n =2
    a = dblquad(psi,xmin,xmax,lambda x : ymin,lambda x: ymax,args = (omega1,0,0,n))[0]
    n = 4
    eff1 = a**2 / dblquad(psi,x.min(),x.max(),lambda x : ymin,lambda x: ymax,args = (omega1,0,0,n))[0]

    n =2
    a = dblquad(psi,xmin,xmax,lambda x : ymin,lambda x: ymax,args = (omega2,0,1,n))[0]
    n = 4
    eff2 = a**2 / dblquad(psi,x.min(),x.max(),lambda x : ymin,lambda x: ymax,args = (omega2,0,1,n))[0]
    
    return (eff1 - 1.61e-10,eff2 - 1.70e-10)

    def integrate(self, Q):

        rDispersion = RealDispersion(self.distribution, self.detuning, Q)
        iDispersion = ImaginaryDispersion(self.distribution, self.detuning, Q)

        realPart, realErr = dblquad(rDispersion,
                                    self.minJx, self.maxJx,
                                    self.minJy, self.maxJy)
        imagPart, imagErr = dblquad(iDispersion,
                                    self.minJx, self.maxJx,
                                    self.minJy, self.maxJy)

        return -1.0/complex(realPart, imagPart)

def calc_theta_norm_2(theta_min, theta_max, dist_min, dist_max, dist_KDE_prior, region, alpha_1, alpha_2, s_crit):

    # args = dist_min, dist_max, dist_KDE_prior, region, alpha_1, alpha_2, s_crit
    args = dist_KDE_prior, region, alpha_1, alpha_2, s_crit

    if region == 1:
        norm = dblquad(calc_integrand, theta_min, theta_max, d_min_region_1, d_max_region_1, args=args)
    else:
        norm = dblquad(calc_integrand, theta_min, theta_max, d_min_region_2, d_max_region_2, args=args)

    # norm = quad(calc_integral, theta_min, theta_max, args=args)

    return norm[0]

def prob_decay_signal(t1a,t1e,t2a,t2e):
	Constant = BrKl*BrKs#propPiPiKs**2*eta**2
	def intensity(t1,t2):
		Part1 = exp(-Gamma_L*t1-Gamma_S*t2)
		Part2 = exp(-Gamma_L*t2-Gamma_S*t1)
		Part3 = -2*(1-zeta)*exp(-.5*(Gamma_L+Gamma_S)*(t1+t2))*cos(Delta_m*(t1-t2))
		return Part1+Part2+Part3
	Int = dblquad(intensity, 0, np.inf, # limits t2
                  lambda x : 0, # limits t1
                  lambda x: np.inf)[0]
	Part = dblquad(intensity, t1a, t1e, # limits t2
                  lambda x : t2a, # limits t1
                  lambda x: t2e)[0]
	return Constant*Part/Int

def checkdf(df,*args,**keywords):
	plt.clf()
	x = np.linspace(0,2,200)
	y = np.linspace(0,2,200)
	X,Y = np.meshgrid(x,y)
	plt.subplot(aspect='equal')
	plt.imshow(df(X,Y, (keywords)), origin='lower', extent=[0,2,0,2])
	plt.plot(x, np.sqrt(2-x**2),c='k', linewidth=3)
	plt.colorbar()
	if df==vel_func_benson: 
		print 'integral in bound region', dblquad(df, 0.,np.sqrt(2), lambda x:0, lambda x:np.sqrt(2.-x**2), args = (keywords['redshift'],))
	elif df==vel_func_gaussian:
		print 'integral in bound region',   dblquad(df, 0.,np.sqrt(2), lambda x:0, lambda x:np.sqrt(2.-x**2), args = (keywords,))
	else: print 'integral in bound region', dblquad(df, 0.,np.sqrt(2), lambda x:0, lambda x:np.sqrt(2.-x**2))

 def test_call(self):
     from scipy.integrate import dblquad
     m = MultiGauss2D(sigmas=[1, 2], norms=[3, 4])
     xy_max = 5 * m.max_sigma  # integration range
     integral = dblquad(m, -xy_max, xy_max,
                        lambda _: -xy_max, lambda _: xy_max)[0]
     assert_almost_equal(integral, 7, decimal=5)

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 test():
    fr = 0.5
    rhor = 400
    qitot = 1e-3
    nitot = 100000
    qifrac = qitot/nitot

    s = p3lookup(rhor, fr, qifrac, 0, 0.01, 'p3lookup_new.npy')

    i1 = gp.gintegral_seg(s.mu, s.l, s.dgr, s.dcr)
    i2, err = quad(lambda x: x**s.mu*np.exp(-s.l*x), s.dgr, s.dcr)

    print '---------SINGLE INTEGRATION---------'
    print 'result from gamma function:', i1
    print 'result from numerical int:', i2
    print 'difference i2-i1:', i2-i1

    cr1 = collision_rate_pyt(s.dgr, s.dcr, s.rhog, fr, s.mu, s.l)
    cr2 = collision_rate(s.dgr, s.dcr, s.rhog, fr, s.mu, s.l)

    print '---------COLLISION RATE---------'
    print 'result from gamma function:', cr1
    print 'result from numerical int:', cr2
    print 'difference i2-i1:', cr2-cr1

    mu2 = 2
    di1 = gp.gintegral_seg(s.mu, s.l, s.dgr, s.dcr)*gp.gintegral_seg(s.mu+1, s.l, s.dgr, s.dcr)
    di2, err = dblquad(lambda x,y: x**s.mu*np.exp(-s.l*x)*y**(s.mu+1)*np.exp(-s.l*y), s.dgr, s.dcr, lambda m:s.dgr, lambda m: s.dcr)

    print '---------DOUBLE INTEGRAL---------'
    print 'result from gamma function:', di1
    print 'result from numerical int:', di2
    print 'difference i2-i1:', di2-di1

    return