def scf_compute_coeffs(dens, N, L, a=1., radial_order=None, costheta_order=None, phi_order=None):
        """        
        NAME:

           scf_compute_coeffs

        PURPOSE:

           Numerically compute the expansion coefficients for a given triaxial density

        INPUT:

           dens - A density function that takes a parameter R, z and phi

           N - size of the Nth dimension of the expansion coefficients

           L - size of the Lth and Mth dimension of the expansion coefficients
           
           a - parameter used to shift the basis functions

           radial_order - Number of sample points of the radial integral. If None, radial_order=max(20, N + 3/2L + 1)

           costheta_order - Number of sample points of the costheta integral. If None, If costheta_order=max(20, L + 1)

           phi_order - Number of sample points of the phi integral. If None, If costheta_order=max(20, L + 1)

        OUTPUT:

           (Acos,Asin) - Expansion coefficients for density dens that can be given to SCFPotential.__init__

        HISTORY:

           2016-05-27 - Written - Aladdin 

        """
        def integrand(xi, costheta, phi):
            l = nu.arange(0, L)[nu.newaxis, :, nu.newaxis]
            m = nu.arange(0, L)[nu.newaxis,nu.newaxis,:]
            r = _xiToR(xi, a)
            R = r*nu.sqrt(1 - costheta**2.)
            z = r*costheta
            Legandre = lpmn(L - 1,L-1,costheta)[0].T[nu.newaxis,:,:]
            dV = (1. + xi)**2. * nu.power(1. - xi, -4.)
            
            
            phi_nl = - a**3*(1. + xi)**l * (1. - xi)**(l + 1.)*_C(xi, N, L)[:,:,nu.newaxis] * Legandre
            
            return dens(R,z, phi) * phi_nl[nu.newaxis, :,:,:]*nu.array([nu.cos(m*phi), nu.sin(m*phi)])*dV
            
               
        Acos = nu.zeros((N,L,L), float)
        Asin = nu.zeros((N,L,L), float)
        
       
        Ksample = [max(N + 3*L//2 + 1,20), max(L + 1,20 ), max(L + 1,20)]
        if radial_order != None:
            Ksample[0] = radial_order
        if costheta_order != None:
            Ksample[1] = costheta_order
        if phi_order != None:
            Ksample[2] = phi_order
        integrated = _gaussianQuadrature(integrand, [[-1., 1.], [-1., 1.], [0, 2*nu.pi]], Ksample = Ksample)
        n = nu.arange(0,N)[:,nu.newaxis, nu.newaxis]
        l = nu.arange(0,L)[nu.newaxis,:, nu.newaxis]
        m = nu.arange(0,L)[nu.newaxis,nu.newaxis,:]
        K = .5*n*(n + 4*l + 3) + (l + 1)*(2*l + 1)
        
        
        Nln = .5*gammaln(l - m + 1) - .5*gammaln(l + m + 1) - (2*l)*nu.log(2)
        NN = nu.e**(Nln)

        NN[nu.where(NN == nu.inf)] = 0 ## To account for the fact that m cant be bigger than l
            
        constants = NN*(2*l + 1.)**.5
        
        lnI = -(8*l + 6)*nu.log(2) + gammaln(n + 4*l + 3) - gammaln(n + 1) - nu.log(n + 2*l + 3./2) - 2*gammaln(2*l + 3./2)
        I = -K*(4*nu.pi) * nu.e**(lnI)
        Acos[:,:,:],Asin[:,:,:] = 2*(I**-1.)[nu.newaxis,:,:,:] * integrated * constants[nu.newaxis,:,:,:]
        
        return Acos, Asin
NAME:

           scf_compute_coeffs

        PURPOSE:

           Numerically compute the expansion coefficients for a given triaxial density

        INPUT:

           dens - A density function that takes a parameter R, z and phi

           N - size of the Nth dimension of the expansion coefficients

           L - size of the Lth and Mth dimension of the expansion coefficients
           
           a - parameter used to shift the basis functions

           radial_order - Number of sample points of the radial integral. If None, radial_order=max(20, N + 3/2L + 1)

           costheta_order - Number of sample points of the costheta integral. If None, If costheta_order=max(20, L + 1)

           phi_order - Number of sample points of the phi integral. If None, If costheta_order=max(20, L + 1)

        OUTPUT:

           (Acos,Asin) - Expansion coefficients for density dens that can be given to SCFPotential.__init__

        HISTORY:

           2016-05-27 - Written - Aladdin