From aebf35d17a7c7fa0bd0ca9995aaeb44d841a36d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jo=C3=A3o=20Faria?= Date: Thu, 25 Jul 2024 14:00:25 +0200 Subject: [PATCH] calculate error on contrast, following pipeline; shift CCFs --- iCCF/gaussian.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/iCCF/gaussian.py b/iCCF/gaussian.py index 96d54e0..15ccf43 100644 --- a/iCCF/gaussian.py +++ b/iCCF/gaussian.py @@ -1,8 +1,9 @@ -from typing import List, Union, Optional +from typing import List, Optional import numpy as np from numpy import exp, log, sqrt from scipy import optimize +from scipy.interpolate import CubicSpline from .utils import numerical_gradient @@ -212,12 +213,12 @@ def FWHMerror(rv, ccf, eccf): ccf : array The values of the CCF profile. eccf : array - The errors on each value of the CCF profile. + The uncertainties on each value of the CCF profile. """ return 2.0 * RVerror(rv, ccf, eccf) -def contrast(rv, ccf): +def contrast(rv, ccf, eccf=None, error=False, **kwargs): """ Calculate the contrast (depth, measured in percentage) of the CCF. @@ -227,6 +228,22 @@ def contrast(rv, ccf): The velocity values where the CCF is defined. ccf : array The values of the CCF profile. + eccf : array + The uncertainties on each value of the CCF profile. """ - A, _, _, continuum = gaussfit(rv, ccf) - return abs(100 * A / continuum) + A, _, _, continuum = gaussfit(rv, ccf, yerr=eccf, **kwargs) + if not error: + return abs(100 * A / continuum) + else: + ## error propagation + # (A, _, _, continuum), (sA, _, _, sc) = gaussfit(rv, ccf, return_errors=True) + # return abs(100 * (A / continuum) * np.hypot(sA / A, sc / continuum)) + ## as in the ESPRESSO pipeline + fwhm = FWHM(rv, ccf, eccf, **kwargs) + snr = RVerror(rv, ccf, eccf) # 1.0 / sqrt(ccf_sum) + return abs(snr * 2 / fwhm * A); + + +def rv_shift(rv, ccf, radial_velocity): + """ Shift a CCF profile by `radial_velocity`. """ + return CubicSpline(rv, ccf)(rv - radial_velocity) \ No newline at end of file