code: def impulse_deltav_hernquist_curvedstream(v,x,b,w,x0,v0,GM,rs):
      
      if len(v.shape) == 1: v= numpy.reshape(v,(1,3))
      if len(x.shape) == 1: x= numpy.reshape(x,(1,3))
      b0 = numpy.cross(w,v0)
      b0 *= b/numpy.sqrt(numpy.sum(b0**2))
      b_ = b0+x-x0
      w = w-v
      wmag = numpy.sqrt(numpy.sum(w**2,axis=1))
      bdotw=numpy.sum(b_*w,axis=1)/wmag
      B = numpy.sqrt(numpy.sum(b_**2,axis=1)-bdotw**2)
      denom= wmag*(B**2-rs**2)
      denom=1./denom
      s = numpy.sqrt(2.*B/(rs+B))
      HernquistXv=numpy.vectorize(HernquistX)
      Xfac = 1.-2.*rs/(rs+B)*HernquistXv(s)
      return -2.0*GM*((b_.T-bdotw*w.T/wmag)*Xfac*denom).T
docs: NAME:
  
         impulse_deltav_plummer_hernquist
  
      PURPOSE:
  
         calculate the delta velocity to due an encounter with a Hernquist sphere in the impulse approximation; allows for arbitrary velocity vectors, and arbitrary position along the stream
  
      INPUT:
  
         v - velocity of the stream (nstar,3)
  
         x - position along the stream (nstar,3)
  
         b - impact parameter
  
         w - velocity of the Hernquist sphere (3)
  
         x0 - point of closest approach
  
         v0 - velocity of point of closest approach
  
         GM - mass of the Hernquist sphere (in natural units)
  
         rs - size of the Hernquist sphere
  
      OUTPUT:
  
         deltav (nstar,3)
  
      HISTORY:
  
         2015-08-13 - SANDERS, using Wyn Evans calculation