def _calc_or(self,Rmean,rperi,rap,E,L,fixed_quad,**kwargs):
     Tr= 0.
     if Rmean > rperi and not fixed_quad:
         Tr+= nu.array(integrate.quadrature(_TrSphericalIntegrandSmall,
                                            0.,m.sqrt(Rmean-rperi),
                                            args=(E,L,self._2dpot,
                                                  rperi),
                                            **kwargs))[0]
     elif Rmean > rperi and fixed_quad:
         Tr+= integrate.fixed_quad(_TrSphericalIntegrandSmall,
                                   0.,m.sqrt(Rmean-rperi),
                                   args=(E,L,self._2dpot,
                                         rperi),
                                   n=10,**kwargs)[0]
     if Rmean < rap and not fixed_quad:
         Tr+= nu.array(integrate.quadrature(_TrSphericalIntegrandLarge,
                                            0.,m.sqrt(rap-Rmean),
                                            args=(E,L,self._2dpot,
                                                  rap),
                                            **kwargs))[0]
     elif Rmean < rap and fixed_quad:
         Tr+= integrate.fixed_quad(_TrSphericalIntegrandLarge,
                                   0.,m.sqrt(rap-Rmean),
                                   args=(E,L,self._2dpot,
                                         rap),
                                   n=10,**kwargs)[0]
     Tr= 2.*Tr
     return 2.*nu.pi/Tr

def G(view, arch):
  '''The Geometry factor for a specific view or solar
  direction based on Myneni III.16. The projection of 1 unit 
  area of leaves within a unit volume onto the plane perpen-
  dicular to the view or solar direction.
  ------------------------------------------------------------
  Input: view - the view or solar zenith angle in radians, 
    arch - archetype, see gl function for description of each.
  Output: The integral of the Geometry function (G).
  '''
  g = lambda angle, view, arch: gl(angle, arch)\
      *psi(angle,view) # the G function as defined in Myneni III.16.
  if isinstance(view, np.ndarray):
    if arch == 's': # avoid integration in case of isometric distr.
      G = np.ones_like(view) * 0.5
      return G
    view = np.where(view > np.pi, 2.*np.pi - view, view)
    G = np.zeros_like(view)
    for j,v in enumerate(view):
      G[j] = fixed_quad(g, 0., np.pi/2., args=(v, arch),n=16)[0]
  else:
    if arch == 's': # avoid integration, see above...
      G = 0.5
      return G
    if view > np.pi: # symmetry of distribution about z axis
      view = np.pi - view
    G = fixed_quad(g, 0., np.pi/2., args=(view, arch),n=16)[0] 
    # integrate g function between 0 to pi/2.
  return G

def compute_exact_solution(outdir='./_output',frame=0):
#------------------------------------
    from scipy.integrate import fixed_quad, quad
    from pyclaw.plotters.data import ClawPlotData
    plotdata = ClawPlotData()
    plotdata.outdir=outdir

    #Read in the solution
    dat = plotdata.getframe(frame)


    ap=ac.AcousticsProblem('setprob.data','sharpclaw.data')

    t = dat.t
    print t
    grid = dat.grids[0]
    xc = dat.center[0]
    dx=dat.d
    xe = dat.edge[0]
    
    print "Computing exact solution for mx = ",dat.n
    qq=zeros([2,size(xc)])
    for ii in range(size(xc)):
	qq[0,ii],dummy = fixed_quad(ap.pvec,xe[ii],xe[ii+1],args=(t,))
	qq[1,ii],dummy = fixed_quad(ap.uvec,xe[ii],xe[ii+1],args=(t,))
	#qq[0,ii],dummy= quad(ap.pvec,xe[ii],xe[ii+1],args=(t,),epsabs=1.e-11,epsrel=1.e-11)
	#qq[1,ii],dummy= quad(ap.uvec,xe[ii],xe[ii+1],args=(t,),epsabs=1.e-11,epsrel=1.e-11)
    qq/=dx
    return qq

	def FM(self,m,Q,M,alpha,r2bar):
		if self.presentation==True: #for the presentation period
			# self consistent equation for m
			if self.fam==True: #familiar stimuli
				fun=lambda x:self.distRates(x)*r2bar(x,m,Q,M,alpha)
				if self.adap_quad==True:
					var,err=integrate.quad(fun,0.,self.rmax)
				else:
					var,err=integrate.fixed_quad(fun,0.,self.rmax,n=self.nquad)
				return var
			else: # novel stimuli 
				fun=lambda x:self.distRates(x)*r2bar(x,Q,M,alpha)
				if self.adap_quad==True:
					var,err=integrate.quad(fun,0.,self.rmax)
				else:
					var,err=integrate.fixed_quad(fun,0.,self.rmax,n=self.nquad)
				return var

		else: #for the delay period
			if self.fam==True: #familiar stimuli
				fun=lambda x:self.distRates(x)*r2bar(x,m,M,alpha)	
				if self.adap_quad==True:
					var,err=integrate.quad(fun,0.,self.rmax)
				else:
					var,err=integrate.fixed_quad(fun,0.,self.rmax,n=self.nquad)
				return var
			else: #novel stimuli
				fun=lambda x:self.distRates(x)*r2bar(x,M,alpha)	
				if self.adap_quad==True:
					var,err=integrate.quad(fun,0.,self.rmax)
				else:
					var,err=integrate.fixed_quad(fun,0.,self.rmax,n=self.nquad)
				return var

 def _calc_op(self,Or,Rmean,rperi,rap,E,L,fixed_quad,**kwargs):
     #Azimuthal period
     I= 0.
     if Rmean > rperi and not fixed_quad:
         I+= nu.array(integrate.quadrature(_ISphericalIntegrandSmall,
                                           0.,m.sqrt(Rmean-rperi),
                                           args=(E,L,self._2dpot,
                                                 rperi),
                                           **kwargs))[0]
     elif Rmean > rperi and fixed_quad:
         I+= integrate.fixed_quad(_ISphericalIntegrandSmall,
                                  0.,m.sqrt(Rmean-rperi),
                                  args=(E,L,self._2dpot,rperi),
                                  n=10,**kwargs)[0]
     if Rmean < rap and not fixed_quad:
         I+= nu.array(integrate.quadrature(_ISphericalIntegrandLarge,
                                           0.,m.sqrt(rap-Rmean),
                                           args=(E,L,self._2dpot,
                                                 rap),
                                           **kwargs))[0]
     elif Rmean < rap and fixed_quad:
         I+= integrate.fixed_quad(_ISphericalIntegrandLarge,
                                  0.,m.sqrt(rap-Rmean),
                                  args=(E,L,self._2dpot,rap),
                                  n=10,**kwargs)[0]
     I*= 2*L
     return I*Or/2./nu.pi

def compute_errors(odir='.',frame=0):
#------------------------------------
    from scipy.integrate import fixed_quad
    # change to output directory and read in solution from fort.q000N
    # for frame N.

    # read in clawpack solution for this frame:
    ap=ac.AcousticsProblem('setprob.data','claw.data')
    os.chdir(odir)
    clawframe = read_clawframe(frame)
    
    t = clawframe.t
    
    # Assume there is only one grid (AMR not used):
    grid = clawframe.grids[0]
    
    print "Computing errors for mx = ",grid.mx
    
    xc = grid.xcenter
    
    qq=zeros([2,size(xc)])
    for ii in range(size(xc)):
	qq[0,ii],dummy = fixed_quad(ap.pvec,grid.xedge[ii],grid.xedge[ii+1],args=(t,))
	qq[1,ii],dummy = fixed_quad(ap.uvec,grid.xedge[ii],grid.xedge[ii+1],args=(t,))
    dx=grid.xedge[1]-grid.xedge[0]
    qq/=dx
    errors = qq[0,:] - grid.q[:,0]
    ion()
    clf()
    plot(xc,qq[0,:],'k',xc,grid.q[:,0],'+b')
    draw()
    ioff()
    return errors

 def _calc_angler(self,Or,r,Rmean,rperi,rap,E,L,vr,fixed_quad,**kwargs):
     if r < Rmean:
         if r > rperi and not fixed_quad:
             wr= Or*integrate.quadrature(_TrSphericalIntegrandSmall,
                                         0.,m.sqrt(r-rperi),
                                         args=(E,L,self._2dpot,rperi),
                                         **kwargs)[0]
         elif r > rperi and fixed_quad:
             wr= Or*integrate.fixed_quad(_TrSphericalIntegrandSmall,
                                         0.,m.sqrt(r-rperi),
                                         args=(E,L,self._2dpot,rperi),
                                         n=10,**kwargs)[0]
         else:
             wr= 0.
         if vr < 0.: wr= 2*m.pi-wr
     else:
         if r < rap and not fixed_quad:
             wr= Or*integrate.quadrature(_TrSphericalIntegrandLarge,
                                         0.,m.sqrt(rap-r),
                                         args=(E,L,self._2dpot,rap),
                                         **kwargs)[0]
         elif r < rap and fixed_quad:
             wr= Or*integrate.fixed_quad(_TrSphericalIntegrandLarge,
                                         0.,m.sqrt(rap-r),
                                         args=(E,L,self._2dpot,rap),
                                         n=10,**kwargs)[0]
         else:
             wr= m.pi
         if vr < 0.:
             wr= m.pi+wr
         else:
             wr= m.pi-wr
     return wr

def Gamma(view, sun, arch, refl, trans):
  '''The one angle Area Scattering Phase Function based on 
  Myneni V.18 and Shultis (17) isotropic scattering assumption. 
  This is the phase function of the scattering in a particular 
  direction based also on the amount of interception in the direction.
  -------------------------------------------------------------
  Input: view - view zenith angle, sun - the solar zenith angle, 
    arch - archetype, see gl function for description, 
    refl - fraction reflected, trans - fraction transmitted.
  Output: Area Scattering Phase function value.
  '''
  '''
  # uncomment/comment the code below for bi-lambetian Gamma.
  B = sun - view # uncomment these lines to run test plot
  gam = (refl + trans)/np.pi/3.*(np.sin(B) - B*np.cos(B)) +\
      trans/np.pi*np.cos(B) # Myneni V.15
  '''
  func = lambda leaf, view, sun, arch, refl, trans: gl(leaf, arch)\
      *(refl*Big_psi(view,sun,leaf,'r') + (trans*Big_psi(view,sun,leaf,'t')))
      # the integral as defined in Myneni V.18.
  if isinstance(sun, np.ndarray):
    # to remove singularity at sun==0.
    sun = np.where(sun==0.,1.0e-10,sun) 
    gam = np.zeros_like(sun)
    for j,s in enumerate(sun):
      gam[j] = fixed_quad(func, 0., np.pi/2.,\
          args=(view,s,arch,refl,trans),n=16)[0]
  else:
    if sun==0.:
      sun = 1.0e-10 # to remove singularity at sun==0.
    gam = fixed_quad(func, 0., np.pi/2.,\
        args=(view,sun,arch,refl,trans),n=16)[0]
    # integrate leaf angles between 0 to pi/2.
  return gam 

  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 f_dmags(self, dmag, s):
     """Calculates the joint probability density of dMag and projected
     separation
     
     Args:
         dmag (float):
             Value of dMag
         s (float):
             Value of projected separation (AU)
     
     Returns:
         f (float):
             Value of joint probability density
     
     """
     
     if (dmag < self.mindmag(s)) or (dmag > self.maxdmag(s)):
         f = 0.0
     else:
         ztest = (s/self.x)**2*10.**(-0.4*dmag)/self.val
         if ztest >= self.zmax:
             f = 0.0
         elif (self.pconst & self.Rconst):
             f = self.f_dmagsz(self.zmin,dmag,s)
         else:
             if ztest < self.zmin:
                 f = integrate.fixed_quad(self.f_dmagsz, self.zmin, self.zmax, args=(dmag, s), n=61)[0]
             else:
                 f = integrate.fixed_quad(self.f_dmagsz, ztest, self.zmax, args=(dmag, s), n=61)[0]
     return f

def onef_dmags(dmag, s, ranges, val, pdfs, funcs, x):
    """Returns joint probability density of s and dmag
    
    Args:
        dmag (float):
            Value of difference in brightness magnitude
        s (float):
            Value of planet-star separation
        ranges (tuple):
            pmin (float): minimum value of geometric albedo
            pmax (float): maximum value of geometric albedo
            Rmin (float): minimum value of planetary radius (km)
            Rmax (float): maximum value of planetary radius (km)
            rmin (float): minimum value of orbital radius (AU)
            rmax (float): maximum value of orbital radius (AU)
            zmin (float): minimum value of pR^2 (km^2)
            zmax (float): maximum value of pR^2 (km^2)
        val (float):
            Value of sin(bstar)**2*Phi(bstar)
        pdfs (tuple):
            Probability density functions for pR^2 and orbital radius
        funcs (tuple):
            Inverse functions of sin(b)**2*Phi(b)
        x (float):
            Conversion factor between AU and km
            
    Returns:
        f (float):
            Joint probability density of s and dmag
            
    """
            
    pmin, pmax, Rmin, Rmax, rmin, rmax, zmin, zmax = ranges
    minrange = (pmax, Rmax, rmin, rmax)
    maxrange = (pmin, Rmin, rmax)
    pconst = pmin == pmax
    Rconst = Rmin == Rmax
    if (dmag < util.mindmag(s, minrange, x)) or (dmag > util.maxdmag(s, maxrange, x)):
        f = 0.0
    elif (pconst & Rconst):
        ranges2 = (rmin, rmax)
        f = onef_dmagsz(pmin*Rmin**2, dmag, s, val, pdfs, ranges2, funcs, x)
    else:
        ztest = (s/x)**2*10.0**(-0.4*dmag)/val
        ranges2 = (rmin, rmax)
        if ztest >= zmax:
            f = 0.0
        else:
            if ztest < zmin:
                f = integrate.fixed_quad(onef_dmagsz, zmin, zmax, args=(dmag, s, val, pdfs, ranges2, funcs, x), n=61)[0]
            else:
                f = integrate.fixed_quad(onef_dmagsz, ztest, zmax, args=(dmag, s, val, pdfs, ranges2, funcs, x), n=61)[0]
            
    return f            

def calc_stream(k, nk, beta, sigma, flag):
  stream = k*0.0+1.0

  if(flag=='exp'):
    for cc in range(0,nk):
      tempi = integrate.fixed_quad(pmu_exp,-1,1,args=(k[cc],sigma,beta),n=8)
      stream[cc] = 0.5*tempi[0] #0.5* to normalize Legendre L_0
  elif(flag=='Gauss'):
    for cc in range(0,nk):
      tempi = integrate.fixed_quad(pmu_gauss,-1,1,args=(k[cc],sigma,beta),n=8)
      stream[cc] = 0.5*tempi[0] #0.5* to normalize Legendre L_0
  
  return stream

def Expectation(model,dist,distparam1,distparam2,Aproj,distparamextra = 0):
    from scipy import integrate
    import numpy as np
# this routine takes in the name of the roof model desired, the distibution of mass to be integrated with respect to,
# and parameters needed to define the distribution
    #fwood,fsteel,fconcrete,fcomposite,flightMetal,ftile= getModel()
    functions = getModel()
    if model.lower()=='all':
    # HERE COMES AN ASSUMPTION TO GIVE ROOF CASUALTY AREA TO ROOF TYPES FOR WHICH DATA IS NOT AVAILABLE:
        # assuming that all roof of similar types have the same casualty areas. e.g all wood roofs are the same
    # the order of models is assumed to be as follows
    # Wood Roof, Wood 1st, Wood 2nd, Steel Roof, Steel 1st, Steel 2nd, Concrete, Concrete 1st, Concrete 2nd, Composite
    # Light Metal, Tile Roof, Tile 1st, Tile 2nd, Car, Open
        
        if dist.lower() =='uniform':
            lowbound = distparam1
            highbound = distparam2
            denc = highbound-lowbound
            if denc>0.:
                c = 1./denc
            elif denc<0.:
                print 'Check bounds in upper and lower bound. roofPenetration.Expectation'
                exit(1)

        elif dist.lower() =='gaussian' or dist.lower() == 'normal' or dist.lower()=='gauss':
            mu = distparam1
            sigma = distparam2
            lowbound = max(-1e3,mu-5.*sigma)
            highbound = min(1e6,mu+5.*sigma)
            fgauss = lambda x: np.exp(-.5*((x-mu)/sigma)**2)/(sigma*(2.*np.pi)**.5)
        
        EretMain = np.zeros((len(functions)))
        for index in range(len(functions)):
            froof = functions[index]
        
            if dist.lower() =='uniform':
                if denc ==0.0:
                    retval = froof(lowbound)
                else:
                    invals = integrate.fixed_quad(froof,lowbound,highbound,n=50)
                    retval = invals[0]*c
            elif dist.lower() =='gaussian' or dist.lower() == 'normal' or dist.lower()=='gauss':
                fvals = lambda x: fgauss(x)*froof(x)
                invals = integrate.fixed_quad(fvals,lowbound,highbound,n=50)
                retval = invals[0]
                
            EretMain[index] = retval
        Eret = np.zeros((16))
        Eret = np.concatenate((3*[EretMain[0]],3*[EretMain[1]],3*[EretMain[2]],[EretMain[3]],[EretMain[4]],3*[EretMain[5]],[EretMain[4]],[Aproj]),1)
        return Eret

 def _calc_anglez(self,Or,Op,ar,z,r,Rmean,rperi,rap,E,L,Lz,vr,axivz,
                  fixed_quad,**kwargs):
     #First calculate psi
     i= nu.arccos(Lz/L)
     sinpsi= z/r/nu.sin(i)
     if sinpsi > 1. and sinpsi < (1.+10.**-7.):
         sinpsi= 1.
     if sinpsi < -1. and sinpsi > (-1.-10.**-7.):
         sinpsi= -1.
     psi= nu.arcsin(sinpsi)
     if axivz > 0.: psi= nu.pi-psi
     psi= psi % (2.*nu.pi)
     #Calculate dSr/dL
     dpsi= Op/Or*2.*nu.pi #this is the full I integral
     if r < Rmean:
         if not fixed_quad:
             wz= L*integrate.quadrature(_ISphericalIntegrandSmall,
                                        0.,m.sqrt(r-rperi),
                                        args=(E,L,self._2dpot,
                                              rperi),
                                        **kwargs)[0]
         elif fixed_quad:
             wz= L*integrate.fixed_quad(_ISphericalIntegrandSmall,
                                        0.,m.sqrt(r-rperi),
                                        args=(E,L,self._2dpot,
                                              rperi),
                                        n=10,**kwargs)[0]
         if vr < 0.: wz= dpsi-wz
     else:
         if not fixed_quad:
             wz= L*integrate.quadrature(_ISphericalIntegrandLarge,
                                        0.,m.sqrt(rap-r),
                                        args=(E,L,self._2dpot,
                                              rap),
                                        **kwargs)[0]
         elif fixed_quad:
             wz= L*integrate.fixed_quad(_ISphericalIntegrandLarge,
                                        0.,m.sqrt(rap-r),
                                        args=(E,L,self._2dpot,
                                              rap),
                                        n=10,**kwargs)[0]
         if vr < 0.:
             wz= dpsi/2.+wz
         else:
             wz= dpsi/2.-wz
     #Add everything
     wz= -wz+psi+Op/Or*ar
     return wz

	def Eg2(self):# mean of g**2
		fun=lambda x:self.distRates2(x)*self.g(x)*self.g(x)
		if self.adap_quad==True:
			var,err=integrate.quad(fun,0.,self.rmaxSig)
		else:
			var,err=integrate.fixed_quad(fun,0.,self.rmaxSig,n=self.nquad)
		return var

	def Ef2(self):# mean of f**2
		fun=lambda x:self.distRates(x)*self.f(x)*self.f(x)
		if self.adap_quad==True:
			var,err=integrate.quad(fun,0.,self.rmax)
		else:
			var,err=integrate.fixed_quad(fun,0.,self.rmax,n=self.nquad)
		return var

  def get_A(self,x,mub2,zetaD,bT2,flav):

    # Eqs. A10 & A11 of PRD83 114042

    D=self.D

    factor1 = 0.5*np.log(4/(mub2*bT2))-D['gamma']
    factor2 =-0.5*(np.log(bT2*mub2)-2*(np.log(2)-D['gamma']))**2\
             -(np.log(bT2*mub2) - 2*(np.log(2)-D['gamma'])) * np.log(zetaD/mub2)

    if flav.startswith('g'):

      C0 = 0
      C1 = lambda z: 0
      C2 = lambda z: D['TF']*(2*(1-2*z*(1-z))*factor1 + 2*z*(1-z))
      C3 = factor2

    else: 

      C0 = 1
      C1 = lambda z: D['CF']*2*factor1*2
      C2 = lambda z: D['CF']*(2*factor1*(-1-z) + 1-z)
      C3 = factor2

    # get alpha strong
    alphaS=self.SC.get_alphaS(mub2)

    D['DI']  = lambda z: C0 +alphaS/(2*np.pi)*(-(1-x)*C1(z)/z**2/(1-z)\
                         + C1(1)*np.log(1-x) + C3)
    D['DII'] = lambda z: alphaS/(2*np.pi)*(1-x)*(C1(z)/z**2/(1-z) + C2(z)/z**2)
    D['ff']  = lambda x: self.CFF.get_FF(x,D['muF2'],flav,D['charge'])
    D['x']   = x

    func=np.vectorize(self.integrand4A)
    return fixed_quad(func,0,1,n=40)[0]

 def comp_s(self, smin, smax, dMag):
     """Calculates completeness by first integrating over dMag and then
     projected separation.
     
     Args:
         smin (ndarray):
             Values of minimum projected separation (AU) from instrument
         smax (ndarray):
             Value of maximum projected separation (AU) from instrument
         dMag (ndarray):
             Planet delta magnitude
     
     Returns:
         comp (ndarray):
             Completeness values
     
     """
     # cast to arrays
     smin = np.array(smin, ndmin=1, copy=False)
     smax = np.array(smax, ndmin=1, copy=False)
     dMag = np.array(dMag, ndmin=1, copy=False)
     
     comp = np.zeros(smin.shape)
     for i in xrange(len(smin)):
         comp[i] = integrate.fixed_quad(self.f_sv, smin[i], smax[i], args=(dMag[i],), n=50)[0]
     # ensure completeness values are between 0 and 1
     comp = np.clip(comp, 0., 1.)
     
     return comp

    def bellman_operator(self, v):
        """
        The Bellman operator.  Including for comparison. Value function
        iteration is not recommended for this problem.  See the reservation
        wage operator below.

        Parameters
        ==============
            v : An approximate value function represented as a
                one-dimensional array.

        """
        # == Simplify names == #
        f, g, beta, c, q = self.f, self.g, self.beta, self.c, self.q 

        vf = LinearNDInterpolator(self.grid_points, v)
        N = len(v)
        new_v = np.empty(N)
        for i in range(N):
            w, pi = self.grid_points[i,:]
            v1 = w / (1 - beta)
            integrand = lambda m: vf(m, q(m, pi)) * (pi * f(m) \
                    + (1 - pi) * g(m))
            integral, error = fixed_quad(integrand, 0, self.w_max)
            v2 = c + beta * integral
            new_v[i] = max(v1, v2)
        return new_v

 def test_scalar(self):
     n = 4
     func = lambda x: x**(2*n - 1)
     expected = 1/(2*n)
     got, _ = fixed_quad(func, 0, 1, n=n)
     # quadrature exact for this input
     assert_allclose(got, expected, rtol=1e-12)