code: def _2ndDerivInt(x,y,z,dens,densDeriv,b2,c2,i,j,glx=None,glw=None):
      
      def integrand(s):
          t= 1/s**2.-1.
          m= numpy.sqrt(x**2./(1.+t)+y**2./(b2+t)+z**2./(c2+t))
          return (densDeriv(m)
                  *(x/(1.+t)*(i==0)+y/(b2+t)*(i==1)+z/(c2+t)*(i==2))
                  *(x/(1.+t)*(j==0)+y/(b2+t)*(j==1)+z/(c2+t)*(j==2))/m\
                      +dens(m)*(i==j)*((1./(1.+t)*(i==0)+1./(b2+t)*(i==1)+1./(c2+t)*(i==2))))\
                      /numpy.sqrt((1.+(b2-1.)*s**2.)*(1.+(c2-1.)*s**2.))
      if glx is None:
          return integrate.quad(integrand,0.,1.)[0]
      else:
          return numpy.sum(glw*integrand(glx))
docs: Integral that gives the 2nd derivative of the potential in x,y,z