def Lookback(M1,L1,K1,function,p):
	if p == 1:
		file = open(function+'_'+str(M1)+'_'+str(L1)+'.txt','wt')
	LB=[]
	LB2=[]
	x=[]
	args = (M1,L1,K1)
	for z in drange(0,5,0.1):
	        x.append(z)
	        result, err = integrate.quad(E_z,0.0,z,args) # Lookback
	        result2, err = integrate.quad(E_z2,z,inf,args) # Age of Universe
	        LB.append(result)
	        LB2.append(result2)
		if p ==1:
	        	file.writelines(str(z)+" "+str(result)+" "+str(result2)+"\n")
	
	if p ==1:
		file.close
	plot(x,LB)
	plot(x,LB2)
	print LB[-1], LB2[-1]
	ylabel('Lookback Time $T_L/T_H$')
	xlabel('Redshift $z$')
	if p!= 1:
		show()

def foc(gov_policies, psi, sig, start=0, end=10):
    # Initialize local variables for government policies for better
    # readability.
    tax = gov_policies[0]
    trans = gov_policies[1]
    result = []

    # Compute different parts of first FOC (respect tax rate) and combine
    part_a = integrate.quad(
        lambda w: dell_u_tax(w, cons(w, tax, psi, trans),
                             hours(w, tax, psi), psi, tax) * lognorm.pdf(
                                 w, s=sig, scale=np.exp(-sig**2 / 2)),
        0, 10                   # set integration borders
    )[0]
    part_b = integrate.quad(
        lambda w: w * hours(w, tax, psi), start, end
    )[0]

    # Compute first part of the second FOC (respect transfers)
    part_c = integrate.quad(
        lambda w: lognorm.pdf(w, s=sig, scale=np.exp(-sig**2 / 2)) *
        dell_u_trans(cons(w, tax, psi, trans), hours(w, tax, psi), psi),
        start, end
    )[0]

    # Store first foc in results vector
    result.append(part_a + part_c * part_b)

    # Compute budget constraint
    bud_const = trans - integrate.quad(
        lambda w: tax * w * hours(w, tax, psi), start, end
    )[0]
    result.append(bud_const)

    return result

def _calculate_levy(x, alpha, beta, cdf=False):
    """ Calculation of Levy stable distribution via numerical integration.
        This is used in the creation of the lookup table. """
    # "0" parameterization as per http://academic2.american.edu/~jpnolan/stable/stable.html
    # Note: fails for alpha=1.0
    #       (so make sure alpha=1.0 isn't exactly on the interpolation grid)
    from scipy import integrate

    C = beta * N.tan(N.pi*0.5*alpha)
    
    def func_cos(u):
        ua = u ** alpha
        if ua > 700.0: return 0.0
        return N.exp(-ua)*N.cos(C*ua-C*u)

    def func_sin(u):
        ua = u ** alpha
        if ua > 700.0: return 0.0
        return N.exp(-ua)*N.sin(C*ua-C*u)

    if cdf:
        # Cumulative density function
        return (integrate.quad(lambda u: u and func_cos(u)/u or 0.0, 0.0, integrate.Inf, weight="sin", wvar=x, limlst=1000)[0]
               + integrate.quad(lambda u: u and func_sin(u)/u or 0.0, 0.0, integrate.Inf, weight="cos", wvar=x, limlst=1000)[0] 
               ) / N.pi + 0.5
    else:
        # Probability density function
        return ( integrate.quad(func_cos, 0.0, integrate.Inf, weight="cos", wvar=x, limlst=1000)[0]
               - integrate.quad(func_sin, 0.0, integrate.Inf, weight="sin", wvar=x, limlst=1000)[0] 
               ) / N.pi

    def radEinstein(self, zs):
        '''
        SDSSObject.radEinstein(z1, z2)
        ==============================
        Estimated Einstein Radius from Single Isothermal Sphere (SIS) model

        Parameters:
            zs: source redshift
        Returns:
            None
        '''
        # Necessary parameter
        H0 = 72e3
        c = 299792458.0
        OmegaM = 0.258
        OmegaL = 0.742
        vdisp = self.vdisp * 1000.0

        # Function needed for numerical computation of angular diameter distance
        def x12(z):
            return 1.0 / np.sqrt((1.0 - OmegaM - OmegaL) * (1.0 + z) *
                                 (1.0 + z) + OmegaM * (1.0 + z) ** 3.0 + OmegaL)
        # compute ang. diam. distances
        Ds = ((c / H0) * quad(x12, 0.0, zs)[0]) / (1.0 + zs)
        Dls = ((c / H0) * quad(x12, self.z, zs)[0]) / (1.0 + zs)
        # return value in arcsec
        coeff = 3600.0 * (180.0 / np.pi)
        return coeff * 4.0 * np.pi * vdisp * vdisp * Dls / (c * c * Ds)

def phi(d, dp):
    """vdW-DF kernel."""
    from scipy.integrate import quad

    kwargs = dict(epsabs=1.0e-6, epsrel=1.0e-6, limit=400)
    cut = 35
    return quad(lambda y: quad(f, 0, cut, (y, d, dp), **kwargs)[0], 0, cut, **kwargs)[0]

 def getConsumerWelfare(self, equilibriumValue):
     """
     Calculate total consumer welfare
         W = \int_{A} u(t,A) dt +\int_{B} u(t,B) dt.
     Can break this up as:
     W = Network Effect Benefits + Heterogeneity Stuff
     W = t^{e} * \nu(t^{e}, A) + (1-t^{e}) * \nu(1-t^{e}) + Heterogeneity
     NB: we are assuming demand has been normalized to 1
     
     @param equilibriumValue: the equilibrium to use when calculating social welfare
     @return: total consumer welfare
     """
     networkAPart = equilibriumValue * self._nu_A(equilibriumValue)
     # Remember _nu_B defined a function of network A size ...
     networkBPart = (1.0 - equilibriumValue) * \
         self._nu_B(equilibriumValue)
     
     hetAPart = integrate.quad(self._het_A, 0, equilibriumValue)[0]
     hetBPart = integrate.quad(self._het_B, equilibriumValue, 1)[0]
     
     priceAPart = equilibriumValue * self._priceFunction_A(self._price_A)
     priceBPart = equilibriumValue * self._priceFunction_B(self._price_B)
     
     # logger.debug('Welfare: (netA, netB, hetA, hetB): ' + \
     #    str((networkAPart, networkBPart, hetAPart, hetBPart)))
     totalWelfare = (networkAPart + networkBPart +
                     hetAPart + hetBPart + priceAPart + priceBPart)
     return totalWelfare

def complex_quad(func, a, b, **kw_args):
    func_re = lambda t: np.real(func(t))
    func_im = lambda t: np.imag(func(t))
    I_re = quad(func_re, a, b, **kw_args)[0]
    I_im = quad(func_im, a, b, **kw_args)[0]
    
    return I_re + 1j*I_im

    def test_b_less_than_a_3(self):
        def f(x):
            return 1.0

        val_1, err_1 = quad(f, 0, 1, weight='alg', wvar=(0, 0))
        val_2, err_2 = quad(f, 1, 0, weight='alg', wvar=(0, 0))
        assert_allclose(val_1, -val_2, atol=max(err_1, err_2))

    def test_ctypes_variants(self):
        sin_0 = get_clib_test_routine('_sin_0', ctypes.c_double,
                                      ctypes.c_double, ctypes.c_void_p)

        sin_1 = get_clib_test_routine('_sin_1', ctypes.c_double,
                                      ctypes.c_int, ctypes.POINTER(ctypes.c_double),
                                      ctypes.c_void_p)

        sin_2 = get_clib_test_routine('_sin_2', ctypes.c_double,
                                      ctypes.c_double)

        sin_3 = get_clib_test_routine('_sin_3', ctypes.c_double,
                                      ctypes.c_int, ctypes.POINTER(ctypes.c_double))

        sin_4 = get_clib_test_routine('_sin_3', ctypes.c_double,
                                      ctypes.c_int, ctypes.c_double)

        all_sigs = [sin_0, sin_1, sin_2, sin_3, sin_4]
        legacy_sigs = [sin_2, sin_4]
        legacy_only_sigs = [sin_4]

        # LowLevelCallables work for new signatures
        for j, func in enumerate(all_sigs):
            callback = LowLevelCallable(func)
            if func in legacy_only_sigs:
                assert_raises(ValueError, quad, callback, 0, pi)
            else:
                assert_allclose(quad(callback, 0, pi)[0], 2.0)

        # Plain ctypes items work only for legacy signatures
        for j, func in enumerate(legacy_sigs):
            if func in legacy_sigs:
                assert_allclose(quad(func, 0, pi)[0], 2.0)
            else:
                assert_raises(ValueError, quad, func, 0, pi)

    def test_b_less_than_a(self):
        def f(x, p, q):
            return p * np.exp(-q*x)

        val_1, err_1 = quad(f, 0, np.inf, args=(2, 3))
        val_2, err_2 = quad(f, np.inf, 0, args=(2, 3))
        assert_allclose(val_1, -val_2, atol=max(err_1, err_2))

    def test_b_less_than_a_2(self):
        def f(x, s):
            return np.exp(-x**2 / 2 / s) / np.sqrt(2.*s)

        val_1, err_1 = quad(f, -np.inf, np.inf, args=(2,))
        val_2, err_2 = quad(f, np.inf, -np.inf, args=(2,))
        assert_allclose(val_1, -val_2, atol=max(err_1, err_2))

  def get_FO(self,x,y,z,Q2,qT2,muR2,muF2,charge,ps='dxdQ2dzdqT2',method='gauss'):
    D=self.D
    D['A1']=1+(2/y-1)**2
    D['A2']=-2
    w2=qT2/Q2*x*z
    w=w2**0.5
    xia_=lambda xib: w2/(xib-z)+x
    xib_=lambda xia: w2/(xia-x)+z

    integrand_xia=lambda xia: self.get_M(xia,xib_(xia),x/xia,z/xib_(xia),Q2,muF2,qT2,charge)
    integrand_xib=lambda xib: self.get_M(xia_(xib),xib,x/xia_(xib),z/xib,Q2,muF2,qT2,charge)

    if method=='quad':
      FO = quad(integrand_xia,x+w,1)[0] + quad(integrand_xib,z+w,1)[0]

    elif method=='gauss':
      integrand_xia=np.vectorize(integrand_xia)
      integrand_xib=np.vectorize(integrand_xib)
      FO = fixed_quad(integrand_xia,x+w,1,n=40)[0] + fixed_quad(integrand_xib,z+w,1,n=40)[0]

    if ps=='dxdQ2dzdqT2':
      s=x*y*Q2
      prefactor = D['alphaEM']**2 * self.SC.get_alphaS(muR2) 
      prefactor/= 2*s**2*Q2*x**2
      prefactor*= D['GeV**-2 -> nb'] 
      return prefactor * FO
    else: 
      print 'ps not inplemented'
      return None

def nc_of_norm():
    f1 = lambda x: x**4
    f2 = lambda x: x**2-x+2
    
    rv = norm(loc = 2 , scale=20)
    print rv.mean()
    print rv.var()
    print rv.moment(1)
    print rv.moment(4)
    #     moments的参数可为m(均值),v(方差值),s(偏度),k(峰度),默认为mv
    print rv.stats(moments = 'mvsk')
# 3）scipy.stats.norm.expect(func,loc=0,scale=1)函数返回func函数相对于正态分布的期望值，其中函数f(x)相对于分布dist的期望值定义为E[x]=Integral(f(x) * dist.pdf(x))
    print(norm.expect(f1, loc=1, scale=2))
    print(norm.expect(f2, loc=2, scale=5))
    
    
    
# （2）lambda argument_list:expression用来表示匿名函数
# （3）numpy.exp(x)计算x的指数
# （4）numpy.inf表示正无穷大
# （5）scipy.integrate.quad(func,a,b)计算func函数从a至b的积分
# （6）scipy.stats.expon(scale)函数返回符合指数分布的参数为scale的随机变量rv
# （7）rv.moment(n,*args,*kwds)返回随机变量的n阶距
        #1st non-center moment of expon distribution whose lambda is 0.5
    E1 = lambda x: x*0.5*np.exp(-x/2)
    #2nd non-center moment of expon distribution whose lambda is 0.5
    E2 = lambda x: x**2*0.5*np.exp(-x/2)
    print(integrate.quad(E1, 0, np.inf))
    print(integrate.quad(E2, 0, np.inf))

 def computeColor(self, filtX, filtY, z):
     """Compute color (flux in filter X - filter Y) of SED at redshift z, return color in magnitudes
     
        @param filtX   lower wavelength filter
        @param filtY   higher wavelength filter
        @param z       redshift of SED
     """
     
     if filtX not in self.filterDict:
         emsg = "Filter " + filtX + " is not in the filter dictionary"
         raise LookupError(emsg)
     if filtY not in self.filterDict:
         emsg = "Filter " + filtY + " is not in the filter dictionary"
         raise LookupError(emsg)
     if filtX == filtY:
         raise ValueError("ERROR! cannot have color as difference between same filter")
     
     # define integral limits by filter extents
     aX, bX = self.filterDict[filtX].returnFilterRange()
     aY, bY = self.filterDict[filtY].returnFilterRange()
     
     int1 = integ.quad(self._integrand1, aX, bX, args=(z, filtX))[0]
     int2 = integ.quad(self._integrand1, aX, bX, args=(z, filtY))[0]
     
     # Do we need this zero-point term?
     int3 = integ.quad(self._integrand2, aX, bX, args=(filtX))[0]
     int4 = integ.quad(self._integrand2, aX, bX, args=(filtY))[0]
     zp = -2.5*math.log10(int4/int3)
     
     return -2.5*math.log10(int1/int2) + zp

 def funcf(E,verbose=False):
     epsilon = 0
     try:
         t = E.shape
         fans = []
         problems = []
         for i in range(len(E)):
             print i+1, ' of ', len(E)
             rapoval = rapo(E[i])
             prefactor = (1./(sqrt(8)*pi**2*(model.Mnorm + 10**Mencgood(log10(rapoval)))))
             temp = intg.quad(finterior,epsilon,1-epsilon,args=(E[i],rapoval),full_output = 1)
             t = temp[0]
             try:
                 if temp[3] != '':
                     if verbose == True:
                         print 'f, E = ',E[i],'message = ',temp[3]
                     problems.append(i)                  
             except IndexError:
                 pass
             fans.append((prefactor*t)[0])
         return array(fans),problems
     except AttributeError:
         rapoval = rapo(E)
         prefactor = (1./(sqrt(8)*pi**2*(model.Mnorm + 10**Mencgood(log10(rapoval)))))
         temp = intg.quad(finterior,epsilon,1-epsilon,args=(E,rapoval),full_output = 1)
         t = temp[0]
         problem = []
         try:
             if temp1[3] != '':
                 if verbose == True:
                     print 'f, E = ',E,'message = ',temp1[3]
                 problem = [E]
         except IndexError:
             pass
         return prefactor*t, problem

 def modelFunc(ti, p, tau, t0):
     p = 0.078
     z = abs(ti - t0) / tau
     def delta(p, r, z):
         if ((r >= (z + p)) or (r <= (z - p))):
            # print "boO"
             return 0
         elif ((r + z) <= p):
             #print "boo dos"
             return 1
         else: 
             #print "ay"
             return ((1/pi) * acos((z**2 - p**2 + r**2)/(2*z*r)))
         
     def I(r):
         mu = (1 - r**2)**0.5
         return 1 - (1 - (mu**0.5))
         
     def function1(r):
         deltaresult = delta(p, r, z) # stored delta value in deltaresult 
         return I(r) * (1 - deltaresult) * (2*r)
     
     def function2(r):
         return I(r) * 2*r    
     
     func1 = integrate.quad(function1, 0, 1)
     func2 = integrate.quad(function2, 0, 1)
     flux = func1[0] / func2[0]
     return flux

    def potential_from_particles(self, nbins = None, r_bins = None):
        """
        Uses the binned density profile from the particles to calculate 
        (approximately) what the total gravitational potential should be.
        This should be compared to the exact potential generated by the 
        N-body calculation. With enough particles, they should be the same.
        """

        r, dens = self.density_profile(nbins, r_bins)
        nbins = np.size(r)
        dens_function = interpolate.interp1d(r, dens)


        integrand_1 = lambda x : x * x * dens_function(x)
        integrand_2 = lambda x :     x * dens_function(x)
        
        rmin, rmax = np.min(r), np.max(r)
        
        pot = np.zeros(nbins)       
        for i in np.arange(nbins):
            

            A = integrate.quad(integrand_1, rmin,  r[i])[0]
            B = integrate.quad(integrand_2,  r[i], rmax)[0]

            pot[i] = A/r[i] + B

        pot = - 4.0 * np.pi * cgs.G * pot

        return r, pot

def complex_integral(func,a,b,args,intype='stupid'):
    real_func = lambda z: np.real(func(z,*args))
    imag_func = lambda z: np.imag(func(z,*args))
    if intype == 'quad':
        real_int = integ.quad(real_func,a,b)
        imag_int = integ.quad(imag_func,a,b)
#        print(real_int)
#        print(imag_int)
        return real_int[0] + 1j * imag_int[0]
    elif intype == 'quadrature':
        real_int = integ.quadrature(real_func,a,b)
        imag_int = integ.quadrature(imag_func,a,b)
#        print(real_int)
#        print(imag_int)
        return real_int[0] + 1j * imag_int[0]
    elif intype == 'romberg':
        real_int = integ.romberg(real_func,a,b)
        imag_int = integ.romberg(imag_func,a,b)
#        print(real_int)
#        print(imag_int)
        return real_int + 1j * imag_int
    elif intype == 'stupid':
        Npoints = 500
        z,dz = np.linspace(a,b,Npoints,retstep=True)
        real_int = np.sum(real_func(z))*dz
        imag_int = np.sum(imag_func(z))*dz
#        print(real_int)
#        print(imag_int)
        return real_int + 1j*imag_int

def C_l(freq1, freq2, Pk_interp, l_vals): #freqs are entered in GHz
    
    f1 = freq1*10**9 #[Hz]
    f2 = freq2*10**9 #[Hz]
    z1 = (nu21/f1)-1
    z2 = (nu21/f2)-1

    one_over_Ez = lambda zprime:1/numpy.sqrt(omg_m*(1+zprime)**3+omg_lambda)
    Dc1 = (c/H0)*integrate.quad(one_over_Ez,0,z1)[0]
    Dc2 = (c/H0)*integrate.quad(one_over_Ez,0,z2)[0]

    C_ls = []

    for i in range(len(l_vals)):

        ans2=0
        for kk in range(len(k_data)):
            #val = (2/numpy.pi)*Pk_interp(k_data[kk])*(k_data[kk]**2)*scipy.special.sph_jn(l_vals[i],k_data[kk]*Dc1)[0][l_vals[i]]*scipy.special.sph_jn(l_vals[i],k_data[kk]*Dc2)[0][l_vals[i]]
            #ans2+=val

            J_sub = (1/2.)*(2*l_vals[i]+1)
            x1 = k_data[kk]*Dc1
            x2 = k_data[kk]*Dc2
            bessel1 = numpy.sqrt(numpy.pi/2)*scipy.special.jn(J_sub,x1)/(numpy.sqrt(x1))
            bessel2 = numpy.sqrt(numpy.pi/2)*scipy.special.jn(J_sub,x2)/(numpy.sqrt(x2))
            val = (2/numpy.pi)*Pk_interp(k_data[kk])*(k_data[kk]**2)*bessel1*bessel2
            ans2 += val
        ans2*=(k_data[1]-k_data[0])

        C_ls.append(ans2)

    return C_ls

def BG(params,t,rho):
  K = params['K'].value
  K2 = params['K2'].value
  Dt = params['Dt'].value 
  rho_0 = params['rho_0'].value
  j1 = params['j1'].value
  Tk = params['Tk'].value
  
  a=numpy.ones(len(t))
  b=numpy.ones(len(t))
  c=numpy.ones(len(t))
  for i in range(len(t)):
    func_ph = lambda x:(x**5)/((numpy.exp(x)-1)*(1-numpy.exp(-x)))#((numpy.sinh(x))**2)
    func_sd = lambda x:(x**3)/((numpy.exp(x)-1)*(1-numpy.exp(-x)))
    func_ee = lambda x:(x**2)/((numpy.exp(x)-1)*(1-numpy.exp(-x)))
    ph = quad(func_ph,0,(Dt/t[i]))
    sd = quad(func_sd,0,(Dt/t[i]))
    ee = quad(func_ee,0,(Dt/t[i]))
    a[i]=ph[0]  #Phonon scattering
    b[i]=sd[0]  # s-d scattering
    c[i]=ee[0]  # e-e scattering
    
  model3 =  rho_0 + K * ((t/Dt)**5) * a + K2 * ((t/Dt)**3) * b + K * ((t/Dt)**2) * c
  model2 =  rho_0 + K * ((t/Dt)**5) * a + K2 * ((t/Dt)**3) * b 
  model1 =  rho_0 + K * ((t/Dt)**5) * a 
  
  x1 = numpy.log(t/Tk)
  p = numpy.pi
  x2 = ((x1**2) + p**2)**0.5
  Py = model1 + j1*(1-(x1/x2))
  return Py-rho