code: def scf_compute_coeffs_spherical(dens, N, a=1., radial_order=None):
          
          numOfParam = 0
          try:
              dens(0)
              numOfParam=1
          except:
              try:
                  dens(0,0)
                  numOfParam=2
              except:
                  numOfParam=3
          param = [0]*numOfParam;
          
          
          def integrand(xi):
              r = _xiToR(xi, a)
              R = r
              param[0] = R
              return a**3. * dens(*param)*(1 + xi)**2. * (1 - xi)**-3. * _C(xi, N, 1)[:,0]
                 
          Acos = nu.zeros((N,1,1), float)
          Asin = None
          
          Ksample = [max(N + 1, 20)]
          
          if radial_order != None:
              Ksample[0] = radial_order
              
          integrated = _gaussianQuadrature(integrand, [[-1., 1.]], Ksample=Ksample)
          n = nu.arange(0,N)
          K = 16*nu.pi*(n + 3./2)/((n + 2)*(n + 1)*(1 + n*(n + 3.)/2.))
          Acos[n,0,0] = 2*K*integrated
          return Acos, Asin
docs: NAME:
  
             scf_compute_coeffs_spherical
  
          PURPOSE:
  
             Numerically compute the expansion coefficients for a given spherical density
  
          INPUT:
  
             dens - A density function that takes a parameter R
  
             N - size of 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 + 1)
  
          OUTPUT:
  
             (Acos,Asin) - Expansion coefficients for density dens that can be given to SCFPotential.__init__
  
          HISTORY:
  
             2016-05-18 - Written - Aladdin