From 1fc9ad6bcaabfdbf57ee13e201077ce7ec689968 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Tue, 26 Nov 2024 11:01:44 +0100 Subject: [PATCH 01/15] for PSF functions with two components, separate the x_c and y_c parameters (functions and jacobians) --- spectractor/extractor/psf.py | 188 +++++++++++++++++++++-------------- 1 file changed, 112 insertions(+), 76 deletions(-) diff --git a/spectractor/extractor/psf.py b/spectractor/extractor/psf.py index 4297cff0b..d8b289e68 100644 --- a/spectractor/extractor/psf.py +++ b/spectractor/extractor/psf.py @@ -230,15 +230,15 @@ def evaluate_moffat1d_jacobian(y, amplitude, y_c, gamma, alpha, norm, dnormda, f return J -@njit(["float32[:](int64[:], float32, float32, float32, float32, float32, float32, float32)", - "float32[:](float32[:], float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) -def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, norm_moffat): # pragma: no cover +@njit(["float32[:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32)", + "float32[:](float32[:], float32, float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) +def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, y_g, sigma, norm_moffat): # pragma: no cover r"""Compute a 1D Moffat-Gaussian function, whose integral is normalised to unity. .. math :: f(y) \propto A \left\lbrace - \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}+ \eta e^{-(y-y_c)^2/(2\sigma^2)}\right\rbrace + \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}+ \eta e^{-(y-y_g)^2/(2\sigma^2)}\right\rbrace \quad\text{ and } \quad \eta < 0, \alpha > 1/2 Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor for the Moffat+Gauss @@ -259,6 +259,8 @@ def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, no Exponent :math:`\alpha` of the Moffat function. eta_gauss: float Relative negative amplitude of the Gaussian function. + y_g: float + Center :math:`y_g` of the Gaussian function. sigma: float Width :math:`\sigma` of the Gaussian function. norm_moffat: float @@ -280,7 +282,7 @@ def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, no >>> eta_gauss = -0.1 >>> sigma = 1 >>> norm = evaluate_moffat1d_normalisation(gamma, alpha) - >>> a = evaluate_moffatgauss1d(y, amplitude=amplitude, y_c=Ny/2, gamma=gamma, alpha=alpha, eta_gauss=eta_gauss, sigma=sigma, norm_moffat=norm) + >>> a = evaluate_moffatgauss1d(y, amplitude=amplitude, y_c=Ny/2, gamma=gamma, alpha=alpha, eta_gauss=eta_gauss, y_g=Ny/2+3, sigma=sigma, norm_moffat=norm) >>> print(f"{np.sum(a):.6f}") 9.966492 >>> a.dtype @@ -301,7 +303,7 @@ def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, no y = np.arange(Ny) amplitude = 10 norm = evaluate_moffat1d_normalisation(gamma, alpha) - a = evaluate_moffatgauss1d(y, amplitude=amplitude, y_c=Ny/2, gamma=5, alpha=2, eta_gauss=-0.1, sigma=1, norm_moffat=norm) + a = evaluate_moffatgauss1d(y, amplitude=amplitude, y_c=Ny/2, gamma=5, alpha=2, eta_gauss=-0.1, y_c=Ny/2+3, sigma=1, norm_moffat=norm) plt.plot(a) plt.grid() plt.xlabel("y") @@ -312,28 +314,24 @@ def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, no yc = y - y_c rr = yc * yc rr_gg = rr / (gamma * gamma) - rr_ss = rr / (sigma * sigma) + rr_ss = (y - y_g) * (y - y_g) / (sigma * sigma) norm = (1. / norm_moffat) + eta_gauss * np.sqrt(2 * np.pi) * sigma a = (1 + rr_gg) ** -alpha + eta_gauss * np.exp(-(rr_ss / 2)) a *= (amplitude / norm) return a -@njit(["float32[:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32, float32)", - "float32[:](float32[:], float32, float32, float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) -def evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alpha2, norm1, norm2): # pragma: no cover +@njit(["float32[:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32, float32, float32)", + "float32[:](float32[:], float32, float32, float32, float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) +def evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, y_c2, gamma2, alpha2, norm1, norm2): # pragma: no cover r"""Compute a 1D DoubleMoffat function, whose integral is normalised to unity. .. math :: f(y) \propto A \left\lbrace - \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}+ \eta e^{-(y-y_c)^2/(2\sigma^2)}\right\rbrace + \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}+ \eta \frac{1}{\left[ 1 +\left(\frac{y-y_{c,2}}{\gamma_2}\right)^2 \right]^{\alpha_2}}\right\rbrace \quad\text{ and } \quad \eta < 0, \alpha > 1/2 - Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor for the Moffat+Gauss - :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} + \eta \sqrt{2\pi} \sigma` is not included as special functions - are not supproted by the numba library. - Parameters ---------- y: array_like @@ -348,6 +346,8 @@ def evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alph Exponent :math:`\alpha_1` of the Moffat function. eta: float Relative amplitude of the second Moffat function. + y_c2: float + Center :math:`y_{c,2}` of the second Moffat function. gamma2: float Width :math:`\gamma_2` of the second Moffat function. alpha2: float @@ -370,14 +370,14 @@ def evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alph >>> amplitude = 10 >>> gamma1 = 5 >>> alpha1 = 2 - >>> eta = 2 + >>> eta = 0.1 >>> gamma2 = 4 >>> alpha2 = 3 >>> norm1 = evaluate_moffat1d_normalisation(gamma1, alpha1) >>> norm2 = evaluate_moffat1d_normalisation(gamma2, alpha2) - >>> a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) + >>> a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, y_c2=Ny/2+3, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) >>> print(f"{np.sum(a):.6f}") - 9.985071 + 9.969373 >>> a.dtype dtype('float32') @@ -402,7 +402,7 @@ def evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alph alpha2 = 3 norm1 = evaluate_moffat1d_normalisation(gamma1, alpha1) norm2 = evaluate_moffat1d_normalisation(gamma2, alpha2) - a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) + a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, y_c2=Ny/2+3, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) plt.plot(a) plt.grid() plt.xlabel("y") @@ -413,22 +413,22 @@ def evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alph yc = y - y_c rr = yc * yc rr_gg1 = rr / (gamma1 * gamma1) - rr_gg2 = rr / (gamma2 * gamma2) + rr_gg2 = (y - y_c2) * (y - y_c2) / (gamma2 * gamma2) norm = 1. / norm1 + eta / norm2 a = (1 + rr_gg1) ** -alpha1 + eta * (1 + rr_gg2) ** -alpha2 a *= (amplitude / norm) return a -@njit(["float32[:,:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) -def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, norm_moffat, dnormda, fixed): # pragma: no cover +@njit(["float32[:,:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) +def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, y_g, sigma, norm_moffat, dnormda, fixed): # pragma: no cover r"""Compute a 1D Moffat-Gaussian Jacobian, whose integral is normalised to unity. .. math :: f(y) \propto A \left\lbrace \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha} \times \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} - - \eta e^{-(y-y_c)^2/(2\sigma^2)}\right\rbrace + - \eta e^{-(y-y_g)^2/(2\sigma^2)}\right\rbrace \quad\text{ and } \quad \eta < 0, \alpha > 1/2 Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor for the Moffat @@ -442,13 +442,15 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, amplitude: float Integral :math:`A` of the function. y_c: float - Center :math:`y_c` of the function. + Center :math:`y_c` of the Moffat function. gamma: float Width :math:`\gamma` of the Moffat function. alpha: float Exponent :math:`\alpha` of the Moffat function. eta_gauss: float Relative negative amplitude of the Gaussian function. + y_g: float + Center :math:`y_g` of the Gaussian function. sigma: float Width :math:`\sigma` of the Gaussian function. norm_moffat: float @@ -476,9 +478,9 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, >>> norm = evaluate_moffat1d_normalisation(gamma, alpha) >>> dnormda = evaluate_moffat1d_normalisation_dalpha(norm, alpha) >>> a = evaluate_moffatgauss1d(y, amplitude=amplitude, y_c=Ny/2, gamma=gamma, alpha=alpha, - ... eta_gauss=eta_gauss, sigma=sigma, norm_moffat=norm) + ... eta_gauss=eta_gauss, y_g=Ny/2+3, sigma=sigma, norm_moffat=norm) >>> J = evaluate_moffatgauss1d_jacobian(y, amplitude=amplitude, y_c=Ny/2, gamma=gamma, alpha=alpha, - ... eta_gauss=eta_gauss, sigma=sigma, norm_moffat=norm, dnormda=dnormda, fixed=np.array([False, False, True, False, False])) + ... eta_gauss=eta_gauss, y_g=Ny/2+3, sigma=sigma, norm_moffat=norm, dnormda=dnormda, fixed=np.array([False, False, True, False, False, False])) >>> J.shape (7, 50) >>> J.dtype @@ -493,9 +495,10 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, """ yc = y - y_c + ys = y - y_g rr = yc * yc rr_gg = rr / (gamma * gamma) - rr_ss = rr / (sigma * sigma) + rr_ss = ys * ys / (sigma * sigma) inv_moffat = 1 / (1 + rr_gg) psf_moffat = inv_moffat ** alpha dpsf_moffat = alpha * inv_moffat * psf_moffat @@ -507,7 +510,7 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, J[0] = (norm / amplitude) * psf # amplitude # fixed x_c so J[1] = 0 if not fixed[2]: - J[2] = (norm / (sigma * sigma)) * yc * eta_gauss * psf_gauss + (2 * norm / (gamma * gamma)) * yc * dpsf_moffat # y_c + J[2] = (2 * norm / (gamma * gamma)) * yc * dpsf_moffat # y_c if not fixed[3]: J[3] = (2 * norm / gamma) * rr_gg * dpsf_moffat - (norm * norm / (amplitude * norm_moffat * gamma)) * psf # gamma if not fixed[4]: @@ -515,19 +518,20 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, if not fixed[5]: J[5] = norm * psf_gauss - (np.sqrt(2 * np.pi) * sigma * norm * norm / amplitude) * psf # eta if not fixed[6]: - J[6] = (-eta_gauss * np.sqrt(2*np.pi) * norm * norm / amplitude) * psf + (eta_gauss * norm / sigma) * rr_ss * psf_gauss # sigma + J[6] = (norm / (sigma * sigma)) * ys * eta_gauss * psf_gauss # y_g + if not fixed[7]: + J[7] = (-eta_gauss * np.sqrt(2*np.pi) * norm * norm / amplitude) * psf + (eta_gauss * norm / sigma) * rr_ss * psf_gauss # sigma return J -@njit(["float32[:,:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) -def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alpha2, norm1, dnormda1, norm2, dnormda2, fixed): # pragma: no cover +@njit(["float32[:,:](int64[:], float32, float32, float32, float32, float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) +def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, y_c2, gamma2, alpha2, norm1, dnormda1, norm2, dnormda2, fixed): # pragma: no cover r"""Compute a 1D DoubleMoffat Jacobian, whose integral is normalised to unity. .. math :: f(y) \propto A \left\lbrace - \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha} \times \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} - - \eta e^{-(y-y_c)^2/(2\sigma^2)}\right\rbrace + \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}+ \eta \frac{1}{\left[ 1 +\left(\frac{y-y_{c,2}}{\gamma_2}\right)^2 \right]^{\alpha_2}}\right\rbrace \quad\text{ and } \quad \eta < 0, \alpha > 1/2 Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor for the Moffat @@ -548,6 +552,8 @@ def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, gam Exponent :math:`\alpha_1` of the Moffat function. eta: float Relative amplitude of the second Moffat function. + y_c2: float + Center :math:`y_{c,2}` of the second Moffat function. gamma2: float Width :math:`\gamma_2` of the second Moffat function. alpha2: float @@ -583,10 +589,10 @@ def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, gam >>> dnormda1 = evaluate_moffat1d_normalisation_dalpha(norm1, alpha1) >>> norm2 = eta * evaluate_moffat1d_normalisation(gamma2, alpha2) >>> dnormda2 = evaluate_moffat1d_normalisation_dalpha(norm2, alpha2) - >>> a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) + >>> a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, y_c2=Ny/2+3, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) >>> J = evaluate_doublemoffat1d_jacobian(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, - ... eta=eta, gamma2=gamma2, alpha2=alpha2, norm1=norm1, dnormda1=dnormda1, norm2=norm2, dnormda2=dnormda2, - ... fixed=np.array([False, False, True, False, False, False, False])) + ... eta=eta, y_c2=Ny/2+3, gamma2=gamma2, alpha2=alpha2, norm1=norm1, dnormda1=dnormda1, norm2=norm2, dnormda2=dnormda2, + ... fixed=np.array([False, False, True, False, False, False, False, False])) >>> J.shape (8, 50) >>> J.dtype @@ -603,7 +609,7 @@ def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, gam yc = y - y_c rr = yc * yc rr_gg1 = rr / (gamma1 * gamma1) - rr_gg2 = rr / (gamma2 * gamma2) + rr_gg2 = (y - y_c2) * (y - y_c2) / (gamma2 * gamma2) inv_moffat1 = 1 / (1 + rr_gg1) psf_moffat1 = inv_moffat1 ** alpha1 dpsf_moffat1 = alpha1 * inv_moffat1 * psf_moffat1 @@ -751,7 +757,7 @@ def evaluate_moffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, fixed): >>> yy, xx = np.mgrid[:Ny, :Nx] >>> amplitude = 10 >>> a = evaluate_moffat2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2) - >>> J = evaluate_moffat2d_jacobian(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2, fixed=np.array([False, False, True, False, False])) + >>> J = evaluate_moffat2d_jacobian(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2, fixed=np.array([False, False, True, False, False, False])) >>> J.shape (5, 2500) >>> J.dtype @@ -786,15 +792,15 @@ def evaluate_moffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, fixed): return J -@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) -def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, sigma): # pragma: no cover +@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) +def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, x_g, y_g, sigma): # pragma: no cover r"""Compute a 2D Moffat+Gauss function, whose integral is normalised to unity. .. math :: f(x, y) = \frac{A}{\frac{\pi \gamma^2}{\alpha-1} + 2 \pi \eta \sigma^2}\left\lbrace \frac{1}{ \left[ 1 +\frac{\left(x-x_c\right)^2+\left(y-y_c\right)^2}{\gamma^2} \right]^\alpha} - + \eta e^{-\left[ \left(x-x_c\right)^2+\left(y-y_c\right)^2\right]/(2 \sigma^2)} + + \eta e^{-\left[ \left(x-x_g\right)^2+\left(y-y_g\right)^2\right]/(2 \sigma^2)} \right\rbrace .. math :: @@ -813,15 +819,19 @@ def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, s amplitude: float Integral :math:`A` of the function. x_c: float - X axis center :math:`x_c` of the function. + X axis center :math:`x_c` of the Moffat function. y_c: float - Y axis center :math:`y_c` of the function. + Y axis center :math:`y_c` of the Moffat function. gamma: float Width :math:`\gamma` of the Moffat function. alpha: float Exponent :math:`\alpha` of the Moffat function. eta_gauss: float Relative negative amplitude of the Gaussian function. + x_g: float + X axis center :math:`x_g` of the Gaussian function. + y_g: float + Y axis center :math:`y_g` of the Gaussian function. sigma: float Width :math:`\sigma` of the Gaussian function. @@ -838,9 +848,9 @@ def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, s >>> yy, xx = np.mgrid[:Ny, :Nx] >>> amplitude = 10 >>> a = evaluate_moffatgauss2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2, - ... eta_gauss=-0.1, sigma=1) + ... eta_gauss=0.1, x_g=Nx/2+3, y_g=Ny/2+3, sigma=2) >>> print(f"{np.sum(a):.6f}") - 9.680574 + 9.692955 >>> a.dtype dtype('float32') @@ -858,12 +868,12 @@ def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, s Ny = 50 yy, xx = np.mgrid[:Nx, :Ny] amplitude = 10 - a = evaluate_moffatgauss2d(xx, yy, amplitude, Nx/2, Ny/2, gamma=5, alpha=2, eta_gauss=-0.1, sigma=1) + a = evaluate_moffatgauss2d(xx, yy, amplitude, Nx/2, Ny/2, gamma=5, alpha=2, eta_gauss=0.1, Nx/2+3, Ny/2+3, sigma=2) im = plt.pcolor(xx, yy, a) plt.grid() plt.xlabel("x") plt.ylabel("y") - plt.colorbar(im, label="Moffat 2D") + plt.colorbar(im, label="Moffat+Gauss 2D") plt.show() """ @@ -871,15 +881,15 @@ def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, s yc = y - y_c rr = xc * xc + yc * yc rr_gg = rr / (gamma * gamma) - rr_ss = rr / (sigma * sigma) + rr_ss = ((x - x_g) * (x - x_g) + (y - y_g) * (y - y_g)) / (sigma * sigma) a = (1 + rr_gg) ** -alpha + eta_gauss * np.exp(-(rr_ss / 2)) norm = (np.pi * gamma * gamma) / (alpha - 1) + eta_gauss * 2 * np.pi * sigma * sigma a *= amplitude / norm return a -@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) -def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamma2, alpha2): # pragma: no cover +@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, float32, float32, float32)"], fastmath=True, cache=True) +def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, x_c2, y_c2, gamma2, alpha2): # pragma: no cover r"""Compute a 2D Double Moffat function, whose integral is normalised to unity. .. math :: @@ -887,7 +897,7 @@ def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamm f(x, y) = \frac{A}{\frac{\pi \gamma_1^2}{\alpha_1-1} + \eta \frac{\pi \gamma_2^2}{\alpha_2-1}}\left\lbrace \frac{1}{ \left[ 1 +\frac{\left(x-x_c\right)^2+\left(y-y_c\right)^2}{\gamma_1^2} \right]^\alpha_1} + \eta \frac{1}{ - \left[ 1 +\frac{\left(x-x_c\right)^2+\left(y-y_c\right)^2}{\gamma_2^2} \right]^\alpha_2} + \left[ 1 +\frac{\left(x-x_{c,2}\right)^2+\left(y-y_{c,2}\right)^2}{\gamma_2^2} \right]^\alpha_2} \right\rbrace .. math :: @@ -915,6 +925,10 @@ def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamm Exponent :math:`\alpha_1` of the Moffat function. eta: float Relative amplitude of the second Moffat function. + x_c2: float + X axis center :math:`x_{c,2}` of the second Moffat function. + y_c2: float + Y axis center :math:`y_{c,2}` of the second Moffat function. gamma2: float Width :math:`\gamma_2` of the second Moffat function. alpha2: float @@ -933,9 +947,9 @@ def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamm >>> yy, xx = np.mgrid[:Ny, :Nx] >>> amplitude = 10 >>> a = evaluate_doublemoffat2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma1=5, alpha1=2, - ... eta=2, gamma2=4, alpha2=3) + ... eta=2, x_c2=Nx/2+3, y_c2=Ny/2+3, gamma2=4, alpha2=3) >>> print(f"{np.sum(a):.6f}") - 9.805085 + 9.804738 >>> a.dtype dtype('float32') @@ -953,7 +967,7 @@ def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamm Ny = 50 yy, xx = np.mgrid[:Nx, :Ny] amplitude = 10 - a = evaluate_doublemoffat2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma1=5, alpha1=2, eta=2, gamma2=4, alpha2=1.5) + a = evaluate_doublemoffat2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma1=5, alpha1=2, eta=2, x_c2=Nx/2+3, y_c2=Ny/2+3, gamma2=4, alpha2=1.5) im = plt.pcolor(xx, yy, a) plt.grid() plt.xlabel("x") @@ -966,7 +980,7 @@ def evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamm yc = y - y_c rr = xc * xc + yc * yc rr_gg1 = rr / (gamma1 * gamma1) - rr_gg2 = rr / (gamma2 * gamma2) + rr_gg2 = ((x - x_c2) * (x - x_c2) + (y - y_c2) * (y - y_c2)) / (gamma2 * gamma2) a = (1 + rr_gg1) ** -alpha1 + eta * (1 + rr_gg2) ** -alpha2 norm = (np.pi * gamma1 * gamma1) / (alpha1 - 1) + eta * (np.pi * gamma2 * gamma2) / (alpha2 - 1) a *= amplitude / norm @@ -1255,8 +1269,8 @@ def evaluate_gauss2d_jacobian(x, y, amplitude, x_c, y_c, sigma, fixed): # pragm -@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) -def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, sigma, fixed): # pragma: no cover +@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) +def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, x_g, y_g, sigma, fixed): # pragma: no cover r"""Compute a 2D Moffat+Gauss Jacobian, whose integral is normalised to unity. Parameters @@ -1277,6 +1291,10 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta Exponent :math:`\alpha` of the Moffat function. eta_gauss: float Relative negative amplitude of the Gaussian function. + x_g: float + X axis center :math:`x_g` of the Gaussian function. + y_g: float + Y axis center :math:`y_g` of the Gaussian function. sigma: float Width :math:`\sigma` of the Gaussian function. fixed: array_like @@ -1296,9 +1314,9 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta >>> yy, xx = np.mgrid[:Ny, :Nx] >>> amplitude = 10 >>> a = evaluate_moffatgauss2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2, - ... eta_gauss=-0.1, sigma=1) + ... eta_gauss=0.1, x_g=Nx/2+3, y_g=Ny/2+3, sigma=1) >>> J = evaluate_moffatgauss2d_jacobian(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2, - ... eta_gauss=-0.1, sigma=1, fixed=np.array([False, False, True, False, False, False, False])) + ... eta_gauss=0.1, x_g=Nx/2+3, y_g=Ny/2+3, sigma=1, fixed=np.array([False, False, True, False, False, False, False, False, False])) >>> J.shape (7, 2500) >>> J.dtype @@ -1315,8 +1333,11 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta xc = (x - x_c).ravel() yc = (y - y_c).ravel() rr = xc * xc + yc * yc + xg = (x - x_g).ravel() + yg = (y - y_g).ravel() + rrg = xg * xg + yg * yg rr_gg = rr / (gamma * gamma) - rr_ss = rr / (sigma * sigma) + rr_ss = rrg / (sigma * sigma) inv_moffat = 1 / (1 + rr_gg) psf_moffat = inv_moffat ** alpha dpsf_moffat = alpha * inv_moffat * psf_moffat @@ -1327,22 +1348,26 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta if not fixed[0]: J[0] = (norm / amplitude) * psf # amplitude if not fixed[1]: - J[1] = (norm / (sigma * sigma)) * xc * eta_gauss * psf_gauss + (2 * norm / (gamma * gamma)) * xc * dpsf_moffat # x_c + J[1] = (2 * norm / (gamma * gamma)) * xc * dpsf_moffat # x_c if not fixed[2]: - J[2] = (norm / (sigma * sigma)) * yc * eta_gauss * psf_gauss + (2 * norm / (gamma * gamma)) * yc * dpsf_moffat # y_c + J[2] = (2 * norm / (gamma * gamma)) * yc * dpsf_moffat # y_c if not fixed[3]: J[3] = (-2 * np.pi * gamma * norm * norm / (alpha-1) / amplitude) * psf + (2 * norm / gamma) * rr_gg * dpsf_moffat # gamma if not fixed[4]: J[4] = (np.pi * gamma * gamma) * norm * norm / (amplitude * (alpha-1) * (alpha-1)) * psf - norm * psf_moffat * np.log(1 + rr_gg) # alpha if not fixed[5]: - J[5] = norm * psf_gauss - (2 * np.pi * sigma * sigma * norm * norm / amplitude) * psf + J[5] = norm * psf_gauss - (2 * np.pi * sigma * sigma * norm * norm / amplitude) * psf # eta if not fixed[6]: - J[6] = (-4 * eta_gauss * np.pi * sigma * norm * norm / amplitude) * psf + (eta_gauss * norm / sigma) * rr_ss * psf_gauss + J[6] = (norm / (sigma * sigma)) * xg * eta_gauss * psf_gauss # x_g + if not fixed[7]: + J[7] = (norm / (sigma * sigma)) * yg * eta_gauss * psf_gauss # y_g + if not fixed[8]: + J[8] = (-4 * eta_gauss * np.pi * sigma * norm * norm / amplitude) * psf + (eta_gauss * norm / sigma) * rr_ss * psf_gauss # sigma return J -@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) -def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamma2, alpha2, fixed): # pragma: no cover +@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32, float32, float32, float32, float32, float32, boolean[:])"], fastmath=True, cache=True) +def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, x_c2, y_c2, gamma2, alpha2, fixed): # pragma: no cover r"""Compute a 2D Double Moffat Jacobian, whose integral is normalised to unity. Parameters @@ -1354,15 +1379,19 @@ def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, amplitude: float Integral :math:`A` of the function. x_c: float - X axis center :math:`x_c` of the function. + X axis center :math:`x_c` of the Moffat function. y_c: float - Y axis center :math:`y_c` of the function. + Y axis center :math:`y_c` of the Moffat function. gamma1: float Width :math:`\gamma_1` of the Moffat function. alpha1: float Exponent :math:`\alpha_1` of the Moffat function. eta: float Relative amplitude of the second Moffat function. + x_c2: float + X axis center :math:`x_{c,2}` of the second Moffat function. + y_c2: float + Y axis center :math:`y_{c,2}` of the second Moffat function. gamma2: float Width :math:`\gamma_2` of the second Moffat function. alpha2: float @@ -1384,9 +1413,9 @@ def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, >>> yy, xx = np.mgrid[:Ny, :Nx] >>> amplitude = 10 >>> a = evaluate_doublemoffat2d(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma1=5, alpha1=2, - ... eta=2, gamma2=1, alpha2=3) + ... eta=2, x_c2=Nx/2+3, y_c2=Ny/2+3, gamma2=1, alpha2=3) >>> J = evaluate_doublemoffat2d_jacobian(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma1=5, alpha1=2, - ... eta=2, gamma2=1, alpha2=3, fixed=np.array([False, False, True, False, False, False, False, False])) + ... eta=2, x_c2=Nx/2+3, y_c2=Ny/2+3, gamma2=1, alpha2=3, fixed=np.array([False, False, True, False, False, False, False, False, False, False])) >>> J.shape (8, 2500) >>> J.dtype @@ -1402,9 +1431,12 @@ def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, """ xc = (x - x_c).ravel() yc = (y - y_c).ravel() + xc2 = (x - x_c2).ravel() + yc2 = (y - y_c2).ravel() rr = xc * xc + yc * yc + rr2 = xc2 * xc2 + yc2 * yc2 rr_gg1 = rr / (gamma1 * gamma1) - rr_gg2 = rr / (gamma2 * gamma2) + rr_gg2 = rr2 / (gamma2 * gamma2) inv_moffat1 = 1 / (1 + rr_gg1) psf_moffat1 = inv_moffat1 ** alpha1 inv_moffat2 = 1 / (1 + rr_gg2) @@ -1417,19 +1449,23 @@ def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, if not fixed[0]: J[0] = (norm / amplitude) * psf # amplitude if not fixed[1]: - J[1] = (2 * norm / (gamma1 * gamma1)) * xc * dpsf_moffat1 + eta * (2 * norm / (gamma2 * gamma2)) * xc * dpsf_moffat2 # x_c + J[1] = (2 * norm / (gamma1 * gamma1)) * xc * dpsf_moffat1 # x_c if not fixed[2]: - J[2] = (2 * norm / (gamma1 * gamma1)) * yc * dpsf_moffat1 + eta * (2 * norm / (gamma2 * gamma2)) * yc * dpsf_moffat2 # y_c + J[2] = (2 * norm / (gamma1 * gamma1)) * yc * dpsf_moffat1 # y_c if not fixed[3]: J[3] = (-2 * np.pi * gamma1 * norm * norm / (alpha1-1) / amplitude) * psf + (2 * norm / gamma1) * rr_gg1 * dpsf_moffat1 # gamma1 if not fixed[4]: J[4] = (np.pi * gamma1 * gamma1) * norm * norm / (amplitude * (alpha1-1) * (alpha1-1)) * psf - norm * psf_moffat1 * np.log(1 + rr_gg1) # alpha1 if not fixed[5]: - J[5] = norm * psf_moffat2 - ((np.pi * gamma2 * gamma2) / (alpha2 - 1) * norm * norm / amplitude) * psf + J[5] = norm * psf_moffat2 - ((np.pi * gamma2 * gamma2) / (alpha2 - 1) * norm * norm / amplitude) * psf # eta if not fixed[6]: - J[6] = eta * (-2 * np.pi * gamma2 * norm * norm / (alpha2-1) / amplitude) * psf + eta * (2 * norm / gamma2) * rr_gg2 * dpsf_moffat2 # gamma2 + J[6] = eta * (2 * norm / (gamma2 * gamma2)) * xc2 * dpsf_moffat2 # x_c2 if not fixed[7]: - J[7] = eta * (np.pi * gamma2 * gamma2) * norm * norm / (amplitude * (alpha2-1) * (alpha2-1)) * psf - eta * norm * psf_moffat2 * np.log(1 + rr_gg2) # alpha2 + J[7] = eta * (2 * norm / (gamma2 * gamma2)) * yc2 * dpsf_moffat2 # y_c2 + if not fixed[8]: + J[8] = eta * (-2 * np.pi * gamma2 * norm * norm / (alpha2-1) / amplitude) * psf + eta * (2 * norm / gamma2) * rr_gg2 * dpsf_moffat2 # gamma2 + if not fixed[9]: + J[9] = eta * (np.pi * gamma2 * gamma2) * norm * norm / (amplitude * (alpha2-1) * (alpha2-1)) * psf - eta * norm * psf_moffat2 * np.log(1 + rr_gg2) # alpha2 return J From 17b2a91f9e57657ab179e4b98336749f1056065a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Tue, 26 Nov 2024 11:43:06 +0100 Subject: [PATCH 02/15] for PSF functions with two components, separate the x_c and y_c parameters (functions and jacobians) --- spectractor/extractor/psf.py | 105 ++++++++++++++++++----------------- 1 file changed, 55 insertions(+), 50 deletions(-) diff --git a/spectractor/extractor/psf.py b/spectractor/extractor/psf.py index d8b289e68..36ed57125 100644 --- a/spectractor/extractor/psf.py +++ b/spectractor/extractor/psf.py @@ -480,9 +480,9 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, >>> a = evaluate_moffatgauss1d(y, amplitude=amplitude, y_c=Ny/2, gamma=gamma, alpha=alpha, ... eta_gauss=eta_gauss, y_g=Ny/2+3, sigma=sigma, norm_moffat=norm) >>> J = evaluate_moffatgauss1d_jacobian(y, amplitude=amplitude, y_c=Ny/2, gamma=gamma, alpha=alpha, - ... eta_gauss=eta_gauss, y_g=Ny/2+3, sigma=sigma, norm_moffat=norm, dnormda=dnormda, fixed=np.array([False, False, True, False, False, False])) + ... eta_gauss=eta_gauss, y_g=Ny/2+3, sigma=sigma, norm_moffat=norm, dnormda=dnormda, fixed=np.array([False, False, True, False, False, False, False, False])) >>> J.shape - (7, 50) + (9, 50) >>> J.dtype dtype('float32') >>> np.allclose(J[2], 0) @@ -505,7 +505,7 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, psf_gauss = np.exp(-(rr_ss / 2)) psf = psf_moffat + eta_gauss * psf_gauss norm = amplitude / ((1. / norm_moffat) + eta_gauss * np.sqrt(2 * np.pi) * sigma) - J = np.zeros((7, y.size), dtype=np.float32) + J = np.zeros((9, y.size), dtype=np.float32) if not fixed[0]: J[0] = (norm / amplitude) * psf # amplitude # fixed x_c so J[1] = 0 @@ -517,10 +517,11 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, J[4] = -norm * psf_moffat * np.log(1 + rr_gg) + psf * (norm * norm / amplitude) * (dnormda / (norm_moffat * norm_moffat)) # alpha if not fixed[5]: J[5] = norm * psf_gauss - (np.sqrt(2 * np.pi) * sigma * norm * norm / amplitude) * psf # eta - if not fixed[6]: - J[6] = (norm / (sigma * sigma)) * ys * eta_gauss * psf_gauss # y_g + # fixed x_c so J[6] = 0 if not fixed[7]: - J[7] = (-eta_gauss * np.sqrt(2*np.pi) * norm * norm / amplitude) * psf + (eta_gauss * norm / sigma) * rr_ss * psf_gauss # sigma + J[7] = (norm / (sigma * sigma)) * ys * eta_gauss * psf_gauss # y_g + if not fixed[8]: + J[8] = (-eta_gauss * np.sqrt(2*np.pi) * norm * norm / amplitude) * psf + (eta_gauss * norm / sigma) * rr_ss * psf_gauss # sigma return J @@ -592,9 +593,9 @@ def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, y_c >>> a = evaluate_doublemoffat1d(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, eta=eta, y_c2=Ny/2+3, gamma2=gamma2, alpha2=alpha2, norm1=norm1, norm2=norm2) >>> J = evaluate_doublemoffat1d_jacobian(y, amplitude=amplitude, y_c=Ny/2, gamma1=gamma1, alpha1=alpha1, ... eta=eta, y_c2=Ny/2+3, gamma2=gamma2, alpha2=alpha2, norm1=norm1, dnormda1=dnormda1, norm2=norm2, dnormda2=dnormda2, - ... fixed=np.array([False, False, True, False, False, False, False, False])) + ... fixed=np.array([False, False, True, False, False, False, False, False, False])) >>> J.shape - (8, 50) + (9, 50) >>> J.dtype dtype('float32') >>> np.allclose(J[2], 0) @@ -607,9 +608,10 @@ def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, y_c """ yc = y - y_c + yc2 = y - y_c2 rr = yc * yc rr_gg1 = rr / (gamma1 * gamma1) - rr_gg2 = (y - y_c2) * (y - y_c2) / (gamma2 * gamma2) + rr_gg2 = yc2 * yc2 / (gamma2 * gamma2) inv_moffat1 = 1 / (1 + rr_gg1) psf_moffat1 = inv_moffat1 ** alpha1 dpsf_moffat1 = alpha1 * inv_moffat1 * psf_moffat1 @@ -618,22 +620,25 @@ def evaluate_doublemoffat1d_jacobian(y, amplitude, y_c, gamma1, alpha1, eta, y_c dpsf_moffat2 = alpha2 * inv_moffat2 * psf_moffat2 psf = psf_moffat1 + eta * psf_moffat2 norm = amplitude / (1/norm1 + eta/norm2) - J = np.zeros((8, y.size), dtype=np.float32) + J = np.zeros((10, y.size), dtype=np.float32) if not fixed[0]: J[0] = (norm / amplitude) * psf # amplitude # fixed x_c so J[1] = 0 if not fixed[2]: - J[2] = (2 * norm / (gamma1 * gamma1)) * yc * dpsf_moffat1 + eta * (2 * norm / (gamma2 * gamma2)) * yc * dpsf_moffat2 # y_c + J[2] = (2 * norm / (gamma1 * gamma1)) * yc * dpsf_moffat1 # y_c if not fixed[3]: J[3] = (2 * norm / gamma1) * rr_gg1 * dpsf_moffat1 - (norm * norm / (amplitude * norm1 * gamma1)) * psf # gamma1 if not fixed[4]: J[4] = -norm * psf_moffat1 * np.log(1 + rr_gg1) + psf * (norm * norm / amplitude) * (dnormda1 / (norm1 * norm1)) # alpha1 if not fixed[5]: J[5] = norm * psf_moffat2 - (1/norm2 * norm * norm / amplitude) * psf # eta - if not fixed[6]: - J[6] = eta * (2 * norm / gamma2) * rr_gg2 * dpsf_moffat2 - (eta * norm * norm / (amplitude * norm2 * gamma2)) * psf # gamma2 + # fixed x_c so J[6] = 0 if not fixed[7]: - J[7] = -eta * norm * psf_moffat2 * np.log(1 + rr_gg2) + psf * (norm * norm / amplitude) * (eta * dnormda2 / (norm2 * norm2)) # alpha2 + J[7] = eta * (2 * norm / (gamma2 * gamma2)) * yc2 * dpsf_moffat2 # y_c2 + if not fixed[8]: + J[8] = eta * (2 * norm / gamma2) * rr_gg2 * dpsf_moffat2 - (eta * norm * norm / (amplitude * norm2 * gamma2)) * psf # gamma2 + if not fixed[9]: + J[9] = -eta * norm * psf_moffat2 * np.log(1 + rr_gg2) + psf * (norm * norm / amplitude) * (eta * dnormda2 / (norm2 * norm2)) # alpha2 return J @@ -1300,7 +1305,6 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta fixed: array_like Array of booleans, with True values for fixed parameters. - Returns ------- J: array_like @@ -1318,7 +1322,7 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta >>> J = evaluate_moffatgauss2d_jacobian(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma=5, alpha=2, ... eta_gauss=0.1, x_g=Nx/2+3, y_g=Ny/2+3, sigma=1, fixed=np.array([False, False, True, False, False, False, False, False, False])) >>> J.shape - (7, 2500) + (9, 2500) >>> J.dtype dtype('float32') >>> np.allclose(J[2], 0) @@ -1344,7 +1348,7 @@ def evaluate_moffatgauss2d_jacobian(x, y, amplitude, x_c, y_c, gamma, alpha, eta psf_gauss = np.exp(-(rr_ss / 2)) psf = psf_moffat + eta_gauss * psf_gauss norm = amplitude / ((np.pi * gamma * gamma) / (alpha - 1) + eta_gauss * 2 * np.pi * sigma * sigma) - J = np.zeros((7, x.size), dtype=np.float32) + J = np.zeros((9, x.size), dtype=np.float32) if not fixed[0]: J[0] = (norm / amplitude) * psf # amplitude if not fixed[1]: @@ -1417,7 +1421,7 @@ def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, >>> J = evaluate_doublemoffat2d_jacobian(xx, yy, amplitude=amplitude, x_c=Nx/2, y_c=Ny/2, gamma1=5, alpha1=2, ... eta=2, x_c2=Nx/2+3, y_c2=Ny/2+3, gamma2=1, alpha2=3, fixed=np.array([False, False, True, False, False, False, False, False, False, False])) >>> J.shape - (8, 2500) + (10, 2500) >>> J.dtype dtype('float32') >>> np.allclose(J[2], 0) @@ -1445,7 +1449,7 @@ def evaluate_doublemoffat2d_jacobian(x, y, amplitude, x_c, y_c, gamma1, alpha1, dpsf_moffat2 = alpha2 * inv_moffat2 * psf_moffat2 psf = psf_moffat1 + eta * psf_moffat2 norm = amplitude / ((np.pi * gamma1 * gamma1) / (alpha1 - 1) + eta * (np.pi * gamma2 * gamma2) / (alpha2 - 1)) - J = np.zeros((8, x.size), dtype=np.float32) + J = np.zeros((10, x.size), dtype=np.float32) if not fixed[0]: J[0] = (norm / amplitude) * psf # amplitude if not fixed[1]: @@ -1987,13 +1991,13 @@ class MoffatGauss(PSF): def __init__(self, values=None, clip=False): PSF.__init__(self, clip=clip) - self.values_default = np.array([1, 0, 0, 3, 2, -0.5, 1, 1]).astype(float) + self.values_default = np.array([1, 0, 0, 3, 2, 0.5, 0, 2, 1, 1]).astype(float) if values is None: values = np.copy(self.values_default) - labels = ["amplitude", "x_c", "y_c", "gamma", "alpha", "eta_gauss", "stddev", "saturation"] - axis_names = ["$A$", r"$x_c$", r"$y_c$", r"$\gamma$", r"$\alpha$", r"$\eta$", r"$\sigma$", "saturation"] + labels = ["amplitude", "x_c", "y_c", "gamma", "alpha", "eta_gauss", "x_g", "y_g", "stddev", "saturation"] + axis_names = ["$A$", r"$x_c$", r"$y_c$", r"$\gamma$", r"$\alpha$", r"$\eta$", r"$x_g$", r"$y_g$", r"$\sigma$", "saturation"] bounds = [(0, np.inf), (-np.inf, np.inf), (-np.inf, np.inf), (1, np.inf), (1.1, 10), - (-1, -5e-3), (1, np.inf), (0, np.inf)] + (-1, np.inf), (-np.inf, np.inf), (-np.inf, np.inf), (1, np.inf), (0, np.inf)] self.params = FitParameters(values=values, labels=labels, axis_names=axis_names, bounds=bounds) def apply_max_width_to_bounds(self, max_half_width=None): @@ -2001,7 +2005,7 @@ def apply_max_width_to_bounds(self, max_half_width=None): self.max_half_width = max_half_width self.params.bounds[2] = (-2 * self.max_half_width, 2 * self.max_half_width) self.params.bounds[3] = (1, self.max_half_width) - self.params.bounds[6] = (1, self.max_half_width) + self.params.bounds[8] = (1, self.max_half_width) def evaluate(self, pixels, values=None): r"""Evaluate the MoffatGauss function. @@ -2028,7 +2032,7 @@ def evaluate(self, pixels, values=None): Examples -------- - >>> p = [2,20,30,4,2,-0.5,1,10] + >>> p = [2,20,30,4,2,0.5,23,33,1,10] >>> psf = MoffatGauss(p) >>> yy, xx = np.mgrid[:50, :60] >>> out = psf.evaluate(pixels=np.array([xx, yy]), values=p) @@ -2038,7 +2042,7 @@ def evaluate(self, pixels, values=None): import matplotlib.pyplot as plt import numpy as np from spectractor.extractor.psf import MoffatGauss - p = [2,20,30,4,2,-0.5,1,10] + p = [2,20,30,4,2,0.5,23,33,1,10] psf = MoffatGauss(p) yy, xx = np.mgrid[:50, :60] out = psf.evaluate(pixels=np.array([xx, yy]), values=p) @@ -2051,10 +2055,10 @@ def evaluate(self, pixels, values=None): """ if values is not None: self.params.values = np.asarray(values).astype(float) - amplitude, x_c, y_c, gamma, alpha, eta_gauss, stddev, saturation = self.params.values + amplitude, x_c, y_c, gamma, alpha, eta_gauss, x_g, y_g, stddev, saturation = self.params.values if pixels.ndim == 3 and pixels.shape[0] == 2: x, y = pixels # .astype(np.float32) # float32 to increase rapidity - out = evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, stddev) + out = evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, x_g, y_g, stddev) if self.clip: out = np.clip(out, 0, saturation) return out @@ -2062,7 +2066,7 @@ def evaluate(self, pixels, values=None): y = np.array(pixels) if alpha > 0.5: norm = evaluate_moffat1d_normalisation(gamma, alpha) - out = evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, stddev, norm_moffat=norm) + out = evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, y_g, stddev, norm_moffat=norm) if self.clip: out = np.clip(out, 0, saturation) return out @@ -2098,10 +2102,10 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr Examples -------- - >>> p = [2,20,30,4,2,-0.5,1,10] + >>> p = [2,20,30,4,2,0.5,23,33,1,10] >>> epsilon = [0.001] * len(p) >>> psf = MoffatGauss(p) - >>> psf.params.fixed = [True, True, False, False, False, False, False, True] # fix amplitude, x_c, saturation + >>> psf.params.fixed = [True, True, False, False, False, False, False, False, False, True] # fix amplitude, x_c, saturation 2D case @@ -2129,7 +2133,7 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr >>> np.allclose(J_ana[0:2], 0) True - .. doctest:: + .. doctest:: :hide: >>> assert J_ana.shape == (len(p), y.size) @@ -2137,7 +2141,7 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr """ if epsilon is None and not analytical: raise ValueError(f"If analytical=False, must give epsilon values for numerical differentiation.") - amplitude, x_c, y_c, gamma, alpha, eta_gauss, sigma, saturation = self.params.values.astype(float) + amplitude, x_c, y_c, gamma, alpha, eta_gauss, x_g, y_g, sigma, saturation = self.params.values.astype(float) J = super().jacobian(pixels, params, epsilon=epsilon, model_input=model_input, analytical=analytical) if not analytical: if model_input is None: @@ -2158,10 +2162,10 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr if pixels.ndim == 1: norm = evaluate_moffat1d_normalisation(gamma, alpha) dnormda = evaluate_moffat1d_normalisation_dalpha(norm, alpha) - J[:-1] = evaluate_moffatgauss1d_jacobian(pixels, amplitude, y_c, gamma, alpha, eta_gauss, sigma, norm, dnormda, fixed=fixed) # [:-1] assume saturation is fixed + J[:-1] = evaluate_moffatgauss1d_jacobian(pixels, amplitude, y_c, gamma, alpha, eta_gauss, y_g, sigma, norm, dnormda, fixed=fixed) # [:-1] assume saturation is fixed else: xx, yy = pixels - J[:-1] = evaluate_moffatgauss2d_jacobian(xx, yy, amplitude, x_c, y_c, gamma, alpha, eta_gauss, sigma, fixed=fixed) # [:-1] assume saturation is fixed + J[:-1] = evaluate_moffatgauss2d_jacobian(xx, yy, amplitude, x_c, y_c, gamma, alpha, eta_gauss, x_g, y_g, sigma, fixed=fixed) # [:-1] assume saturation is fixed return J @@ -2169,13 +2173,14 @@ class DoubleMoffat(PSF): def __init__(self, values=None, clip=False): PSF.__init__(self, clip=clip) - self.values_default = np.array([1, 0, 0, 3, 2, 0.005, 3, 1.5, 10]).astype(float) + self.values_default = np.array([1, 0, 0, 3, 2, 0.005, 0, 2, 3, 1.5, 10]).astype(float) if values is None: values = np.copy(self.values_default) - labels = ["amplitude", "x_c", "y_c", "gamma1", "alpha1", "eta", "gamma2", "alpha2", "saturation"] - axis_names = ["$A$", r"$x_c$", r"$y_c$", r"$\gamma_1$", r"$\alpha_1$", r"$\eta$", r"$\gamma_2$", r"$\alpha_2$", "saturation"] + labels = ["amplitude", "x_c", "y_c", "gamma1", "alpha1", "eta", "x_c2", "y_c2", "gamma2", "alpha2", "saturation"] + axis_names = ["$A$", r"$x_c$", r"$y_c$", r"$\gamma_1$", r"$\alpha_1$", r"$\eta$", r"$x_{c,2}$", r"$y_{c,2}$", + r"$\gamma_2$", r"$\alpha_2$", "saturation"] bounds = [(0, np.inf), (-np.inf, np.inf), (-np.inf, np.inf), (1, np.inf), (1.1, 10), - (0, 1), (1, np.inf), (1.1, 10), (0, np.inf)] + (-1, 1), (-np.inf, np.inf), (-np.inf, np.inf), (1, np.inf), (1.1, 10), (0, np.inf)] self.params = FitParameters(values=values, labels=labels, axis_names=axis_names, bounds=bounds) def apply_max_width_to_bounds(self, max_half_width=None): @@ -2183,7 +2188,7 @@ def apply_max_width_to_bounds(self, max_half_width=None): self.max_half_width = max_half_width self.params.bounds[2] = (-2 * self.max_half_width, 2 * self.max_half_width) self.params.bounds[3] = (1, self.max_half_width) - self.params.bounds[5] = (1, self.max_half_width) + self.params.bounds[7] = (1, self.max_half_width) def evaluate(self, pixels, values=None): r"""Evaluate the DoubleMoffat function. @@ -2210,7 +2215,7 @@ def evaluate(self, pixels, values=None): Examples -------- - >>> p = [2,20,30,4,2,2,5,1.5,20] + >>> p = [2,20,30,4,2,2,23,33,5,1.5,20] >>> psf = DoubleMoffat(p) >>> yy, xx = np.mgrid[:50, :60] >>> out = psf.evaluate(pixels=np.array([xx, yy]), values=p) @@ -2220,7 +2225,7 @@ def evaluate(self, pixels, values=None): import matplotlib.pyplot as plt import numpy as np from spectractor.extractor.psf import DoubleMoffat - p = [2,20,30,4,2,2,5,1.5,20] + p = [2,20,30,4,2,2,23,33,5,1.5,20] psf = DoubleMoffat(p) yy, xx = np.mgrid[:50, :60] out = psf.evaluate(pixels=np.array([xx, yy]), values=p) @@ -2233,10 +2238,10 @@ def evaluate(self, pixels, values=None): """ if values is not None: self.params.values = np.asarray(values).astype(float) - amplitude, x_c, y_c, gamma1, alpha1, eta, gamma2, alpha2, saturation = self.params.values + amplitude, x_c, y_c, gamma1, alpha1, eta, x_c2, y_c2, gamma2, alpha2, saturation = self.params.values if pixels.ndim == 3 and pixels.shape[0] == 2: x, y = pixels # .astype(np.float32) # float32 to increase rapidity - out = evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, gamma2, alpha2) + out = evaluate_doublemoffat2d(x, y, amplitude, x_c, y_c, gamma1, alpha1, eta, x_c2, y_c2, gamma2, alpha2) if self.clip: out = np.clip(out, 0, saturation) return out @@ -2245,7 +2250,7 @@ def evaluate(self, pixels, values=None): if alpha1 > 0.5 and alpha2 > 0.5: norm1 = evaluate_moffat1d_normalisation(gamma1, alpha1) norm2 = evaluate_moffat1d_normalisation(gamma2, alpha2) - out = evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, gamma2, alpha2, norm1, norm2) + out = evaluate_doublemoffat1d(y, amplitude, y_c, gamma1, alpha1, eta, y_c2, gamma2, alpha2, norm1, norm2) if self.clip: out = np.clip(out, 0, saturation) return out @@ -2281,10 +2286,10 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr Examples -------- - >>> p = [2,20,30,4,2,0.5,3,3,20] + >>> p = [2,20,30,4,2,2,23,33,5,1.5,20] >>> epsilon = [0.001] * len(p) >>> psf = DoubleMoffat(p) - >>> psf.params.fixed = [True, True, False, False, False, False, False, False, True] # fix amplitude, x_c, saturation + >>> psf.params.fixed = [True, True, False, False, False, False, False, False, False, False, True] # fix amplitude, x_c, saturation 2D case @@ -2320,7 +2325,7 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr """ if epsilon is None and not analytical: raise ValueError(f"If analytical=False, must give epsilon values for numerical differentiation.") - amplitude, x_c, y_c, gamma1, alpha1, eta, gamma2, alpha2, saturation = self.params.values.astype(float) + amplitude, x_c, y_c, gamma1, alpha1, eta, x_c2, y_c2, gamma2, alpha2, saturation = self.params.values.astype(float) J = super().jacobian(pixels, params, epsilon=epsilon, model_input=model_input, analytical=analytical) if not analytical: if model_input is None: @@ -2343,10 +2348,10 @@ def jacobian(self, pixels, params, epsilon=None, model_input=None, analytical=Tr norm2 = evaluate_moffat1d_normalisation(gamma2, alpha2) dnormda1 = evaluate_moffat1d_normalisation_dalpha(norm1, alpha1) dnormda2 = evaluate_moffat1d_normalisation_dalpha(norm2, alpha2) - J[:-1] = evaluate_doublemoffat1d_jacobian(pixels, amplitude, y_c, gamma1, alpha1, eta, gamma2, alpha2, norm1, dnormda1, norm2, dnormda2, fixed=fixed) # [:-1] assume saturation is fixed + J[:-1] = evaluate_doublemoffat1d_jacobian(pixels, amplitude, y_c, gamma1, alpha1, eta, y_c2, gamma2, alpha2, norm1, dnormda1, norm2, dnormda2, fixed=fixed) # [:-1] assume saturation is fixed else: xx, yy = pixels - J[:-1] = evaluate_doublemoffat2d_jacobian(xx, yy, amplitude, x_c, y_c, gamma1, alpha1, eta, gamma2, alpha2, fixed=fixed) # [:-1] assume saturation is fixed + J[:-1] = evaluate_doublemoffat2d_jacobian(xx, yy, amplitude, x_c, y_c, gamma1, alpha1, eta, x_c2, y_c2, gamma2, alpha2, fixed=fixed) # [:-1] assume saturation is fixed return J From b38935c3df43e8d2c0aba597e7afd84691f2e80c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Wed, 4 Dec 2024 10:01:23 +0100 Subject: [PATCH 03/15] remove negative of min values from check_statistical_fit --- spectractor/extractor/images.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/spectractor/extractor/images.py b/spectractor/extractor/images.py index 7e8af86ff..a0527d1b4 100644 --- a/spectractor/extractor/images.py +++ b/spectractor/extractor/images.py @@ -421,10 +421,8 @@ def check_statistical_error(self): raise AttributeError(f"Noise map must be in ADU units to be plotted and analyzed. " f"Currently self.units={self.units}.") data = np.copy(self.data) - min_noz = np.min(data[data > 0]) - data[data <= 0] = min_noz - y = self.err.flatten() ** 2 - x = data.flatten() + y = self.err[data > 0].flatten() ** 2 + x = data[data > 0].flatten() fit, cov, model = fit_poly1d(x, y, order=1) gain = 1 / fit[0] read_out = np.sqrt(fit[1]) * gain From 2ea5e03b3babdebbc6e028990fc3a6804406fed0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Wed, 4 Dec 2024 11:51:18 +0100 Subject: [PATCH 04/15] put zeros if no second order when saving --- spectractor/extractor/spectrum.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/spectractor/extractor/spectrum.py b/spectractor/extractor/spectrum.py index 88c3b87e8..07262fedd 100644 --- a/spectractor/extractor/spectrum.py +++ b/spectractor/extractor/spectrum.py @@ -691,7 +691,13 @@ def save_spectrum(self, output_file_name, overwrite=False): if extname == "SPEC_COV": hdus[extname].data = self.cov_matrix elif extname == "ORDER2": - hdus[extname].data = [self.lambdas, self.data_next_order, self.err_next_order] + if self.data_next_order is None: + data_next_order = np.zeros_like(self.data) + err_next_order = np.zeros_like(self.err) + else: + data_next_order = self.data_next_order + err_next_order = self.err_next_order + hdus[extname].data = [self.lambdas, data_next_order, err_next_order] elif extname == "ORDER0": hdus[extname].data = self.target.image hdus[extname].header["IM_X0"] = self.target.image_x0 From d90b56383a4ba9693c7e7535c6d2881615be8116 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Sun, 8 Dec 2024 16:53:09 +0100 Subject: [PATCH 05/15] debug SPECTRACTOR_FIT_TARGET_CENTROID = guess --- spectractor/extractor/images.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spectractor/extractor/images.py b/spectractor/extractor/images.py index a0527d1b4..a0a0f9f0f 100644 --- a/spectractor/extractor/images.py +++ b/spectractor/extractor/images.py @@ -959,7 +959,7 @@ def find_target(image, guess=None, rotated=False, widths=[parameters.XWINDOW, pa sub_image_y0 = target_pixcoords[1] - y0 + Dy elif parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "guess": Dx, Dy = widths - sub_image_subtracted, x0, y0, Dx, Dy, sub_errors = find_target_init(image=image, guess=target_pixcoords, + sub_image_subtracted, x0, y0, Dx, Dy, sub_errors = find_target_init(image=image, guess=guess, rotated=rotated, widths=(Dx, Dy)) theX, theY = guess sub_image_x0 = theX - x0 + Dx From fe3ae6db6aba270c1b82d777146c27cc7df50c2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:44:04 +0100 Subject: [PATCH 06/15] remove photutils exception warning --- spectractor/extractor/background.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/spectractor/extractor/background.py b/spectractor/extractor/background.py index 4fa382998..be9a6a6b1 100644 --- a/spectractor/extractor/background.py +++ b/spectractor/extractor/background.py @@ -5,6 +5,7 @@ import os import copy +import photutils.utils.exceptions from spectractor import parameters from spectractor.tools import fit_poly1d_outlier_removal, fit_poly2d_outlier_removal, plot_image_simple @@ -15,7 +16,8 @@ from scipy.signal import medfilt2d from scipy import ndimage from scipy.interpolate import RegularGridInterpolator - +import warnings +warnings.filterwarnings("ignore", category=photutils.utils.exceptions.NoDetectionsWarning) def _from_bkgd_interp_to_func(bgd_model_func_interp): def bgd_model_func(x, y): From e2be774cba3d74ea5b63af98f4679cd4095a4c23 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:45:45 +0100 Subject: [PATCH 07/15] better initialisation for order0 centroid fit --- spectractor/extractor/images.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spectractor/extractor/images.py b/spectractor/extractor/images.py index a0a0f9f0f..9f50a67e9 100644 --- a/spectractor/extractor/images.py +++ b/spectractor/extractor/images.py @@ -1184,7 +1184,7 @@ def find_target_Moffat2D(image, sub_image_subtracted, sub_errors=None): # fit a 2D star profile close to this position psf = Moffat(clip=True) total_flux = np.sum(sub_image_subtracted) - psf.params.values[:3] = [total_flux, avX, avY] + psf.params.values[:5] = [total_flux, avX, avY, np.mean([sigX, sigY]), 3] psf.params.values[-1] = image.saturation if image.target_star2D is not None: psf.params.values = image.target_star2D.params.values From 23bf6ccf20a6357ec0aa9d520ae03fad5c1ed709 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:46:23 +0100 Subject: [PATCH 08/15] if order0 method is guess, use the rotated value for the posiotion in the rotated image --- spectractor/extractor/images.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spectractor/extractor/images.py b/spectractor/extractor/images.py index 9f50a67e9..f5112b2cf 100644 --- a/spectractor/extractor/images.py +++ b/spectractor/extractor/images.py @@ -909,7 +909,7 @@ def find_target(image, guess=None, rotated=False, widths=[parameters.XWINDOW, pa if parameters.PdfPages: parameters.PdfPages.savefig() - if parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "fit" or rotated: + if parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "fit": if target_pixcoords[0] == -1 and target_pixcoords[1] == -1: if guess is None: raise ValueError(f"Guess parameter must be set if WCS solution is not found.") From a296d9a653202e5b132661d0b45bc46c57c5a9b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:46:26 +0100 Subject: [PATCH 09/15] if order0 method is guess, use the rotated value for the posiotion in the rotated image --- spectractor/extractor/images.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/spectractor/extractor/images.py b/spectractor/extractor/images.py index f5112b2cf..fe683a824 100644 --- a/spectractor/extractor/images.py +++ b/spectractor/extractor/images.py @@ -959,6 +959,11 @@ def find_target(image, guess=None, rotated=False, widths=[parameters.XWINDOW, pa sub_image_y0 = target_pixcoords[1] - y0 + Dy elif parameters.SPECTRACTOR_FIT_TARGET_CENTROID == "guess": Dx, Dy = widths + if rotated: + guess2 = find_target_after_rotation(image) + x0 = int(guess2[0]) + y0 = int(guess2[1]) + guess = [x0, y0] sub_image_subtracted, x0, y0, Dx, Dy, sub_errors = find_target_init(image=image, guess=guess, rotated=rotated, widths=(Dx, Dy)) theX, theY = guess From 6c230862076bb8be384c750effb587b3687203cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:46:37 +0100 Subject: [PATCH 10/15] new plot order0 --- spectractor/extractor/images.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/spectractor/extractor/images.py b/spectractor/extractor/images.py index fe683a824..81fa4ea4b 100644 --- a/spectractor/extractor/images.py +++ b/spectractor/extractor/images.py @@ -1220,6 +1220,8 @@ def find_target_Moffat2D(image, sub_image_subtracted, sub_errors=None): star2D = psf.evaluate(pixels=np.array([X, Y])) plot_image_simple(ax1, data=sub_image_subtracted, scale="lin", title="", units=image.units, target_pixcoords=[new_avX, new_avY], vmin=vmin, vmax=vmax) + plt.axvline(avX, color='r', linestyle='-', lw=2, label="profile**4 guess") + plt.axhline(avY, color='r', linestyle='-', lw=2) ax1.legend(loc=1) ax1.text(0.05, 0.05, f'Data', color="white", From fcb9e2c709da104568bf1f632c041d3b330633e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:47:54 +0100 Subject: [PATCH 11/15] ensure spectrum.lambdas_binwidths is computed in calibrate_spectrum --- spectractor/extractor/spectrum.py | 1 + 1 file changed, 1 insertion(+) diff --git a/spectractor/extractor/spectrum.py b/spectractor/extractor/spectrum.py index 07262fedd..23fd2ae70 100644 --- a/spectractor/extractor/spectrum.py +++ b/spectractor/extractor/spectrum.py @@ -1893,6 +1893,7 @@ def calibrate_spectrum(spectrum, with_adr=False, niter=5, grid_search=False): """ with_adr = int(with_adr) + spectrum.lambdas_binwidths = np.gradient(spectrum.lambdas) if spectrum.units != "ADU/s": # go back in ADU/s to remove previous lambda*dlambda normalisation spectrum.convert_from_flam_to_ADUrate() distance = spectrum.chromatic_psf.get_algebraic_distance_along_dispersion_axis() From 0c1864b180f06017a6117158d689f133d637117d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 13:48:15 +0100 Subject: [PATCH 12/15] adapt to new PSF MoffatGauss with free gaussian centroid position --- spectractor/extractor/chromaticpsf.py | 42 +++++++++++++++++---------- spectractor/extractor/psf.py | 10 +++++-- 2 files changed, 33 insertions(+), 19 deletions(-) diff --git a/spectractor/extractor/chromaticpsf.py b/spectractor/extractor/chromaticpsf.py index befe351c4..7d141fa54 100644 --- a/spectractor/extractor/chromaticpsf.py +++ b/spectractor/extractor/chromaticpsf.py @@ -222,9 +222,10 @@ def generate_test_poly_params(self): >>> s = ChromaticPSF(psf, Nx=5, Ny=4, deg=1, saturation=20000) >>> params = s.generate_test_poly_params() - .. doctest:: + .. doctest:: :hide: - >>> assert(np.all(np.isclose(params,[10, 50, 100, 150, 200, 0, 0, 0, 0, 5, 0, 2, 0, -0.4, 0, 1, 0, 20000]))) + + >>> assert(np.all(np.isclose(params,[10, 50, 100, 150, 200, 0, 0, 0, 0, 5, 0, 2, 0, 0.4, 0, 0, 0, 4, 0, 1, 0, 20000]))) """ if not isinstance(self.psf, MoffatGauss) and not isinstance(self.psf, Moffat): @@ -241,7 +242,9 @@ def generate_test_poly_params(self): params += [0.] * (self.degrees['gamma'] - 1) + [0, 5] # gamma params += [0.] * (self.degrees['alpha'] - 1) + [0, 2] # alpha if isinstance(self.psf, MoffatGauss): - params += [0.] * (self.degrees['eta_gauss'] - 1) + [0, -0.4] # eta_gauss + params += [0.] * (self.degrees['eta_gauss'] - 1) + [0, 0.4] # eta_gauss + params += [0.] * (self.degrees['x_g'] - 1) + [0, 0] # x mean + params += [0.] * (self.degrees['y_g'] - 1) + [0, 4] # y mean params += [0.] * (self.degrees['stddev'] - 1) + [0, 1] # stddev params += [self.saturation] # saturation poly_params = np.zeros_like(params) @@ -1116,7 +1119,7 @@ def from_profile_params_to_poly_params(self, profile_params, indices=None): .. doctest:: :hide: - >>> assert(np.all(np.isclose(profile_params[0], [10, 0, 50, 5, 2, -0.4, 1, 8e3]))) + >>> assert(np.all(np.isclose(profile_params[0], [10, 0, 50, 5, 2, 0.4, 0, 53, 1, 8e3]))) From the profile parameters to the polynomial parameters: @@ -1143,9 +1146,9 @@ def from_profile_params_to_poly_params(self, profile_params, indices=None): delta = 0 if name != 'amplitude': weights = np.copy(amplitude)[indices] - if name == 'x_c': + if 'x_' in name: delta = self.x0 - if name == 'y_c': + if 'y_' in name: delta = self.y0 if parameters.PSF_POLY_TYPE == "legendre": fit = np.polynomial.legendre.legfit(poly_x, profile_params[indices, k] - delta, @@ -1255,16 +1258,16 @@ def from_poly_params_to_profile_params(self, poly_params, apply_bounds=False): .. doctest:: :hide: - >>> assert np.allclose(profile_params[0], [10, 0, 50, 5, 2, -0.4, 1, 8e3], rtol=1e-3, atol=1e-3) + >>> assert np.allclose(profile_params[0], [10, 0, 50, 5, 2, 0.4, 0, 54, 1, 8e3], rtol=1e-3, atol=1e-3) From the profile parameters to the polynomial parameters: - >>> profile_params = s.from_profile_params_to_poly_params(profile_params) + >>> poly_params = s.from_profile_params_to_poly_params(profile_params) .. doctest:: :hide: - >>> assert np.allclose(profile_params, poly_params_test) + >>> assert np.allclose(poly_params, poly_params_test) From the polynomial parameters to the profile parameters without Moffat amplitudes: @@ -1272,7 +1275,7 @@ def from_poly_params_to_profile_params(self, poly_params, apply_bounds=False): .. doctest:: :hide: - >>> assert np.allclose(profile_params[0], [1, 0, 50, 5, 2, -0.4, 1, 8e3]) + >>> assert np.allclose(profile_params[0], [1, 0, 50, 5, 2, 0.4, 0, 54, 1, 8e3]) """ length = len(self.table) @@ -1299,9 +1302,9 @@ def from_poly_params_to_profile_params(self, poly_params, apply_bounds=False): elif parameters.PSF_POLY_TYPE == "polynomial": profile_params[:, k] = np.polynomial.polynomial.polyval(poly_x, p) shift += self.degrees[name] + 1 - if name == 'x_c': + if 'x_' in name: profile_params[:, k] += self.x0 - if name == 'y_c': + if 'y_' in name: profile_params[:, k] += self.y0 if apply_bounds: for k, name in enumerate(self.psf.params.labels): @@ -1546,6 +1549,8 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model gamma True alpha True eta_gauss True + x_g True + y_g True stddev True saturation True >>> assert(not np.any(np.isclose(s.table['flux_sum'][3:6], np.zeros(s.Nx)[3:6], rtol=1e-3))) @@ -1573,8 +1578,6 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model # guess = [2 * np.nanmax(signal), middle, 0.5 * fwhm, 2, 0, 0.1 * fwhm, saturation] signal_sum = np.nanmax(signal) guess[0] = signal_sum - guess[1] = xmax_index - guess[2] = middle guess[-1] = saturation # bounds = [(0.1 * maxi, 10 * maxi), (middle - w, middle + w), (0.1, min(fwhm, Ny // 2)), (0.1, self.alpha_max), # (-1, 0), @@ -1584,8 +1587,13 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model initial_bounds = np.copy(psf.params.bounds) bounds = np.copy(psf.params.bounds) bounds[0] = (0.1 * signal_sum, 2 * signal_sum) - bounds[2] = (middle - w, middle + w) bounds[-1] = (0, 2 * saturation) + for k, label in enumerate(psf.params.labels): + if "x_" in label: + guess[k] = xmax_index + if "y_" in label: + guess[k] = middle + bounds[k] = (middle - w, middle + w) # moffat_guess = [2 * np.nanmax(signal), middle, 0.5 * fwhm, 2] # moffat_bounds = [(0.1 * maxi, 10 * maxi), (middle - w, middle + w), (0.1, min(fwhm, Ny // 2)), (0.1, 10)] # fit = fit_moffat1d_outlier_removal(index, signal, sigma=sigma, niter=2, @@ -1630,7 +1638,9 @@ def fit_transverse_PSF1D_profile(self, data, err, w, ws, pixel_step=1, bgd_model # bounds[0] = (0.1 * np.nanstd(bgd), 2 * np.nanmax(y[middle - ws[0]:middle + ws[0]])) bounds[0] = (0.1 * signal_sum, 1.5 * signal_sum) guess[0] = signal_sum - guess[1] = x + for k, label in enumerate(psf.params.labels): + if "x_" in label: + guess[k] = x # if guess[4] > -1: # guess[0] = np.max(signal) / (1 + guess[4]) # std = np.sqrt(np.nansum(pdf * (index - mean) ** 2)) diff --git a/spectractor/extractor/psf.py b/spectractor/extractor/psf.py index 36ed57125..b07936050 100644 --- a/spectractor/extractor/psf.py +++ b/spectractor/extractor/psf.py @@ -1991,7 +1991,7 @@ class MoffatGauss(PSF): def __init__(self, values=None, clip=False): PSF.__init__(self, clip=clip) - self.values_default = np.array([1, 0, 0, 3, 2, 0.5, 0, 2, 1, 1]).astype(float) + self.values_default = np.array([1, 0, 0, 3, 2, 0.5, 0, 3, 2, 1]).astype(float) if values is None: values = np.copy(self.values_default) labels = ["amplitude", "x_c", "y_c", "gamma", "alpha", "eta_gauss", "x_g", "y_g", "stddev", "saturation"] @@ -2005,6 +2005,7 @@ def apply_max_width_to_bounds(self, max_half_width=None): self.max_half_width = max_half_width self.params.bounds[2] = (-2 * self.max_half_width, 2 * self.max_half_width) self.params.bounds[3] = (1, self.max_half_width) + self.params.bounds[7] = (-2 * self.max_half_width, 2 * self.max_half_width) self.params.bounds[8] = (1, self.max_half_width) def evaluate(self, pixels, values=None): @@ -2188,7 +2189,8 @@ def apply_max_width_to_bounds(self, max_half_width=None): self.max_half_width = max_half_width self.params.bounds[2] = (-2 * self.max_half_width, 2 * self.max_half_width) self.params.bounds[3] = (1, self.max_half_width) - self.params.bounds[7] = (1, self.max_half_width) + self.params.bounds[7] = (-2 * self.max_half_width, 2 * self.max_half_width) + self.params.bounds[8] = (1, self.max_half_width) def evaluate(self, pixels, values=None): r"""Evaluate the DoubleMoffat function. @@ -2627,7 +2629,9 @@ def __init__(self, psf, data, data_errors, bgd_model_func=None, jacobian_analyti self.Nx = 1 self.psf.apply_max_width_to_bounds(self.Ny) self.pixels = np.arange(self.Ny, dtype=int) - self.params.fixed[1] = True + for k, label in enumerate(self.params.labels): + if "x_" in label: + self.params.fixed[k] = True else: raise ValueError(f"Data array must have dimension 1 or 2. Here pixels.ndim={data.ndim}.") From 1a3b727f2cef15ce56a83dfea699992d30485776 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 16:09:20 +0100 Subject: [PATCH 13/15] adapt to multiplt x_ and y_ parameters --- spectractor/extractor/chromaticpsf.py | 29 ++++++++++++++++----------- spectractor/extractor/extractor.py | 7 +++++-- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/spectractor/extractor/chromaticpsf.py b/spectractor/extractor/chromaticpsf.py index 7d141fa54..a617db475 100644 --- a/spectractor/extractor/chromaticpsf.py +++ b/spectractor/extractor/chromaticpsf.py @@ -202,7 +202,9 @@ def crop_table(self, new_Nx): def set_polynomial_degrees(self, deg): self.deg = deg self.degrees = {key: deg for key in self.psf.params.labels} - self.degrees["x_c"] = max(self.degrees["x_c"], 1) + for key in self.psf.params.labels: + if "x_" in key: + self.degrees[key] = max(self.degrees[key], 1) self.degrees['saturation'] = 0 def generate_test_poly_params(self): @@ -1453,9 +1455,9 @@ def plot_summary(self, truth=None): for i, name in enumerate(self.psf.params.labels): coeffs = [self.params.values[k] for k in range(self.params.ndim) if name in self.params.labels[k]] delta = 0 - if name == 'x_c': + if "x_" in name: delta = self.x0 - if name == 'y_c': + if "y_" in name: delta = self.y0 if parameters.PSF_POLY_TYPE == "legendre": PSF_models.append(np.polynomial.polynomial.legval(rescale_x_to_legendre(all_pixels), coeffs) + delta) @@ -1867,7 +1869,8 @@ def fit_chromatic_psf(self, data, mask=None, bgd_model_func=None, data_errors=No self.cov_matrix = np.copy(w.amplitude_cov_matrix) # add background crop to y_c - self.params.values[w.Nx + w.y_c_0_index] += w.bgd_width + for y0_index in self.y0_params: + self.params.values[self.Nx + y0_index] += w.bgd_width # fill results self.psf.apply_max_width_to_bounds(max_half_width=w.Ny + 2 * w.bgd_width) @@ -1895,7 +1898,7 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None, axis_names=list(np.copy(chromatic_psf.params.axis_names[length:])), fixed=None, truth=truth, filename=file_name) for k, par in enumerate(params.labels): - if "x_c" in par or "saturation" in par: + if "x_" in par or "saturation" in par: params.fixed[k] = True FitWorkspace.__init__(self, params, file_name=file_name, verbose=verbose, plot=plot, live_fit=live_fit) self.my_logger = set_logger(self.__class__.__name__) @@ -1905,10 +1908,10 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None, self.bgd_model_func = bgd_model_func self.analytical = analytical self.poly_params = np.copy(self.chromatic_psf.params.values) - self.y_c_0_index = -1 + self.y0_params = [] for k, par in enumerate(params.labels): - if par == "y_c_0": - self.y_c_0_index = k + if "y_" in par and "_0" in par: + self.y0_params.append(k) break # prepare the fit @@ -1946,8 +1949,8 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None, yy, xx = np.mgrid[:self.Ny, :self.Nx] self.pixels = np.asarray([xx, yy]) - - self.poly_params[self.Nx + self.y_c_0_index] -= self.bgd_width + for y0_index in self.y0_params: + self.poly_params[self.Nx + y0_index] -= self.bgd_width self.profile_params = self.chromatic_psf.from_poly_params_to_profile_params(self.poly_params) self.data_before_mask = np.copy(self.data) self.mask_before_mask = list(np.copy(self.mask)) @@ -2284,7 +2287,8 @@ def simulate(self, *shape_params): # linear regression for the amplitude parameters # prepare the vectors poly_params = np.concatenate([np.ones(self.Nx), shape_params]) - poly_params[self.Nx + self.y_c_0_index] -= self.bgd_width + for y0_index in self.y0_params: + self.poly_params[self.Nx + y0_index] -= self.bgd_width profile_params = self.chromatic_psf.from_poly_params_to_profile_params(poly_params, apply_bounds=True) profile_params[:self.Nx, 0] = 1 profile_params[:self.Nx, 1] = np.arange(self.Nx) @@ -2358,7 +2362,8 @@ def simulate(self, *shape_params): if self.amplitude_priors_method == "fixed": self.model = self.chromatic_psf.evaluate(self.pixels, poly_params) # self.chromatic_psf.params.values is updated in evaluate(): reset to original values - self.chromatic_psf.params.values[self.Nx + self.y_c_0_index] += self.bgd_width + for y0_index in self.y0_params: + self.chromatic_psf.params.values[self.Nx + y0_index] += self.bgd_width self.poly_params = np.copy(poly_params) self.profile_params = np.copy(profile_params) self.model_err = np.zeros_like(self.model) diff --git a/spectractor/extractor/extractor.py b/spectractor/extractor/extractor.py index 847bec1a3..408b49ede 100644 --- a/spectractor/extractor/extractor.py +++ b/spectractor/extractor/extractor.py @@ -477,8 +477,11 @@ def simulate(self, *params): # Fill spectrogram trace as a function of the pixel column x self.psf_profile_params[order][:, 0] = self.params[f"A{order}"] * self.tr[k](self.lambdas) - self.psf_profile_params[order][:, 1] = dispersion_law.real + self.spectrum.spectrogram_x0 - self.psf_profile_params[order][:, 2] += dispersion_law.imag - self.bgd_width + for p, label in enumerate(self.spectrum.chromatic_psf.psf.params.labels): + if "x_" in label: + self.psf_profile_params[order][:, p] = dispersion_law.real + self.spectrum.spectrogram_x0 + if "y_" in label: + self.psf_profile_params[order][:, p] += dispersion_law.imag - self.bgd_width # Matrix filling M_order = self.spectrum.chromatic_psf.build_sparse_M(self.pixels, self.psf_profile_params[order], From 5dfb14788a4ec4f031f5116fa78d208240cae7b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Fri, 20 Dec 2024 16:21:20 +0100 Subject: [PATCH 14/15] adapt to multiplt x_ and y_ parameters --- spectractor/extractor/chromaticpsf.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/spectractor/extractor/chromaticpsf.py b/spectractor/extractor/chromaticpsf.py index a617db475..a841462ed 100644 --- a/spectractor/extractor/chromaticpsf.py +++ b/spectractor/extractor/chromaticpsf.py @@ -1736,6 +1736,7 @@ def fit_chromatic_psf(self, data, mask=None, bgd_model_func=None, data_errors=No >>> parameters.PIXWIDTH_BACKGROUND = 10 >>> parameters.PIXWIDTH_SIGNAL = 30 >>> parameters.DEBUG = True + >>> parameters.VERBOSE = True Build a mock spectrogram with random Poisson noise using the full 2D PSF model: @@ -1869,7 +1870,7 @@ def fit_chromatic_psf(self, data, mask=None, bgd_model_func=None, data_errors=No self.cov_matrix = np.copy(w.amplitude_cov_matrix) # add background crop to y_c - for y0_index in self.y0_params: + for y0_index in w.y0_params: self.params.values[self.Nx + y0_index] += w.bgd_width # fill results @@ -1912,7 +1913,6 @@ def __init__(self, chromatic_psf, data, data_errors, mode, bgd_model_func=None, for k, par in enumerate(params.labels): if "y_" in par and "_0" in par: self.y0_params.append(k) - break # prepare the fit self.Ny, self.Nx = self.data.shape @@ -2288,7 +2288,7 @@ def simulate(self, *shape_params): # prepare the vectors poly_params = np.concatenate([np.ones(self.Nx), shape_params]) for y0_index in self.y0_params: - self.poly_params[self.Nx + y0_index] -= self.bgd_width + poly_params[self.Nx + y0_index] -= self.bgd_width profile_params = self.chromatic_psf.from_poly_params_to_profile_params(poly_params, apply_bounds=True) profile_params[:self.Nx, 0] = 1 profile_params[:self.Nx, 1] = np.arange(self.Nx) From 8373a8ab52023409d4eb0c184ef45eeb1160703f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9my=20Neveu?= Date: Tue, 28 Jan 2025 07:47:30 +0100 Subject: [PATCH 15/15] add hg ar lines for calib --- spectractor/extractor/spectroscopy.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/spectractor/extractor/spectroscopy.py b/spectractor/extractor/spectroscopy.py index c741b0eb1..5f1e37dfe 100644 --- a/spectractor/extractor/spectroscopy.py +++ b/spectractor/extractor/spectroscopy.py @@ -584,7 +584,7 @@ def build_detected_line_table(self, amplitude_units="", calibration_only=False): HG6 = Line(365.015, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) HG7 = Line(404.656, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) HG8 = Line(407.783, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) -HG9 = Line(435.833, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02]) +HG9 = Line(435.833, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) HG10 = Line(546.074, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) HG11 = Line(576.960, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) HG12 = Line(579.066, atmospheric=False, label=r'$Hg$', label_pos=[0.007, 0.02], use_for_calibration=True) @@ -592,18 +592,18 @@ def build_detected_line_table(self, amplitude_units="", calibration_only=False): AR2 = Line(706.722, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR3 = Line(714.704, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR4 = Line(727.294, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR5 = Line(738.393, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR6 = Line(750.387, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR7 = Line(763.511, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) +AR5 = Line(738.393, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) +AR6 = Line(750.387, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) +AR7 = Line(763.511, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) AR8 = Line(772.376, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR9 = Line(794.818, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR10 = Line(800.616, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR11 = Line(811.531, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) +AR9 = Line(794.818, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) +AR10 = Line(800.616, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) +AR11 = Line(811.531, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) AR12 = Line(826.452, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR13 = Line(842.465, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR14 = Line(852.144, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR15 = Line(866.794, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) -AR16 = Line(912.297, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) +AR16 = Line(912.297, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02], use_for_calibration=True) AR17 = Line(922.450, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR18 = Line(965.7786, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02]) AR19 = Line(978.4503, atmospheric=False, label=r'$Ar$', label_pos=[0.007, 0.02])