docs: def power_law_anisotropy(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
          """
          equation (19) in Suyu+ 2010
          :param r:
          :return:
          """
          # first term
          theta_E = kwargs_profile['theta_E']
          gamma = kwargs_profile['gamma']
          r_ani = kwargs_anisotropy['r_ani']
          a = 0.551 * kwargs_light['r_eff']
          rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
          prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
          prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
          hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
          hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
          fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
          sigma2_dim_less = prefac1 * prefac2 * fac
          return sigma2_dim_less * (self.cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000)**2