code: def _interpolate_v(p, r, v):
  
      
      # interpolate v (p should be in table)
      # ordinate: y**2
      # abcissa:  1./v
  
      # find the 3 closest v values
      # only p >= .9 have table values for 1 degree of freedom.
      # The boolean is used to index the tuple and append 1 when
      # p >= .9
      v0, v1, v2 = _select_vs(v, p)
  
      # y = f - 1.
      y0_sq = (_func(A[(p,v0)], p, r, v0) + 1.)**2.
      y1_sq = (_func(A[(p,v1)], p, r, v1) + 1.)**2.
      y2_sq = (_func(A[(p,v2)], p, r, v2) + 1.)**2.
  
      # if v2 is inf set to a big number so interpolation
      # calculations will work
      if v2 > 1e38: v2 = 1e38
  
      # transform v
      v_, v0_, v1_, v2_ = 1./v, 1./v0, 1./v1, 1./v2
  
      # calculate derivatives for quadratic interpolation
      d2 = 2.*((y2_sq-y1_sq)/(v2_-v1_) - \
               (y0_sq-y1_sq)/(v0_-v1_)) / (v2_-v0_)
      if (v2_ + v0_) >= (v1_ + v1_):
          d1 = (y2_sq-y1_sq) / (v2_-v1_) - 0.5*d2*(v2_-v1_)
      else:
          d1 = (y1_sq-y0_sq) / (v1_-v0_) + 0.5*d2*(v1_-v0_)
      d0 = y1_sq
  
      # calculate y
      y = math.sqrt((d2/2.)*(v_-v1_)**2. + d1*(v_-v1_)+ d0)
  
      return y
docs: interpolates v based on the values in the A table for the
      scalar value of r and th