From a214980fc2fae4e5bda906ac5f77d5f1a047c383 Mon Sep 17 00:00:00 2001 From: "Sergey E. Koposov" Date: Wed, 7 Feb 2024 23:04:08 +0000 Subject: [PATCH] improve numerical accuracy of resolution matrix --- py/rvspecfit/read_grid.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/py/rvspecfit/read_grid.py b/py/rvspecfit/read_grid.py index a1b1609..4140da5 100644 --- a/py/rvspecfit/read_grid.py +++ b/py/rvspecfit/read_grid.py @@ -40,13 +40,13 @@ def gau_integrator(A, B, x1, x2, l1, l2, s): Erf = scipy.special.erf Pi = np.pi return (-((Sqrt(2 / Pi) * s * - (2 * B + A * (l1 + x1))) / E**((l1 - x1)**2 / (2. * s**2))) + - (Sqrt(2 / Pi) * s * (2 * B + A * (l2 + x1))) / E**((l2 - x1)**2 / + (2 * B + A * (l1 + x1))) * E**(-(l1 - x1)**2 / (2. * s**2))) + + (Sqrt(2 / Pi) * s * (2 * B + A * (l2 + x1))) * E**(-(l2 - x1)**2 / (2. * s**2)) + - (Sqrt(2 / Pi) * s * (2 * B + A * (l1 + x2))) / E**((l1 - x2)**2 / + (Sqrt(2 / Pi) * s * (2 * B + A * (l1 + x2))) * E**(-(l1 - x2)**2 / (2. * s**2)) - (Sqrt(2 / Pi) * s * - (2 * B + A * (l2 + x2))) / E**((l2 - x2)**2 / (2. * s**2)) + x1 * + (2 * B + A * (l2 + x2))) * E**(-(l2 - x2)**2 / (2. * s**2)) + x1 * (2 * B + A * x1) * Erf( (l1 - x1) / (Sqrt(2) * s)) - x1 * (2 * B + A * x1) * Erf( (l2 - x1) / (Sqrt(2) * s)) + 2 * B * l1 * Erf( @@ -98,8 +98,15 @@ def pix_integrator(x1, x2, l1, l2, s): # the reason for two integrations because we have a function # (f1*(-x+x2)+ f2*(x-x1))/(x2-x1) # - coeff1 = gau_integrator(-1. / (x2 - x1), x2 / (x2 - x1), x1, x2, l1, l2, s) - coeff2 = gau_integrator(1. / (x2 - x1), -x1 / (x2 - x1), x1, x2, l1, l2, s) + offset = x1 + # I put in offset, as integration is ambivalent to offsets, but + # gau_integrator has numerical issues when x1,x2 is large + coeff1 = gau_integrator(-1. / (x2 - x1), (x2 - offset) / (x2 - x1), + x1 - offset, x2 - offset, l1 - offset, l2 - offset, + s) + coeff2 = gau_integrator(1. / (x2 - x1), -(x1 - offset) / (x2 - x1), + x1 - offset, x2 - offset, l1 - offset, l2 - offset, + s) return coeff1, coeff2