code: def theta2rho(self, theta_E, gamma):
          
          fac1 = np.sqrt(np.pi) * special.gamma(1. / 2 * (-1 + gamma)) / special.gamma(gamma / 2.) * 2 / (3 - gamma)
          fac2 = theta_E**(gamma - 1)
          rho0 = fac2 / fac1
          return rho0
docs: converts projected density parameter (in units of deflection) into 3d density parameter
          :param theta_E:
          :param gamma:
          :return: