code: def HernquistX(s):
      
      if(s<0.):
          raise ValueError("s must be positive in Hernquist X function")
      elif(s<1.):
          return numpy.log((1+numpy.sqrt(1-s*s))/s)/numpy.sqrt(1-s*s)
      elif(s==1.):
          return 1.
      else:
          return numpy.arccos(1./s)/numpy.sqrt(s*s-1)
docs: Computes X function from equations (33) & (34) of Hernquist (1990)