SIE not valid at q=1

[ConnorStoneAstro]

SIE returns nan when evaluated with q=1. Numerically unstable due to 1/(1-q^2) term, which analytically should cancel. @AlexandreAdam has taylor expansion which solves this issue near q=1, just need a where statement to flip between them

def approximate_deflection_angles(self, r_ein, q, phi, x0, y0):
        b, s = self._param_conv(q, r_ein)
        # rotate to major/minor axis coordinates
        theta1, theta2 = self.rotated_and_shifted_coords(x0, y0, phi)
        psi = torch.sqrt(q ** 2 * (s ** 2 + theta1 ** 2) + theta2 ** 2)
        alpha1 = b * theta1 / (psi + s)
        alpha2 = b * theta2 / (psi + q**2 * s)
        # # rotate back to original orientation of coordinate system
        alpha1, alpha2 = self._rotate(alpha1, alpha2, -phi)
        return torch.concat([alpha1, alpha2], dim=1)  # stack alphas into tensor of shape [batch_size, 2, pix, pix]

    def analytical_deflection_angles(self, r_ein, q, phi, x0, y0):
        b, s = self._param_conv(q, r_ein)
        # rotate to major/minor axis coordinates
        theta1, theta2 = self.rotated_and_shifted_coords(x0, y0, phi)
        psi = torch.sqrt(q ** 2 * (s ** 2 + theta1 ** 2) + theta2 ** 2)
        alpha1 = b / np.sqrt(1. - q ** 2) * torch.atan(np.sqrt(1. - q ** 2) * theta1 / (psi + s))
        alpha2 = b / np.sqrt(1. - q ** 2) * torch.atanh(np.sqrt(1. - q ** 2) * theta2 / (psi + s * q ** 2))
        # # rotate back
        alpha1, alpha2 = self._rotate(alpha1, alpha2, -phi)
        return torch.concat([alpha1, alpha2], dim=1)

    def _param_conv(self, q, r_ein):
        r_ein_conv = 2. * q * r_ein / np.sqrt(1. + q ** 2)
        b = r_ein_conv * np.sqrt((1 + q ** 2) / 2)
        s = self.s_scale * np.sqrt((1 + q ** 2) / (2 * q ** 2))
        return b, s

[ConnorStoneAstro]

On further investigation it seems our SIE is not simply an elliptical version of our SIS, though both match with lenstronomy.

[ConnorStoneAstro]

Ok after further investigation the SIE is fine. And a new PR #301 closes this issue