activity-spirou/tools.py at main · piacortes/activity-spirou · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
"""
@author: Pia Cortes-Zuleta <pia.cortes@lam.fr>
Laboratoire d'Astropysique de Marseille

Description: Tools needed to compute stellar activity indicators.
This includes getRV, bisector, gauss and bigauss profiles.

Last updated: 30/08/2021
    """
import numpy as np
from scipy.optimize import curve_fit
from scipy.interpolate import InterpolatedUnivariateSpline
from uncertainties import ufloat
from uncertainties.umath import *

def getRV(rv,ccf):
    """ Computes RV using a gauss profile fitting on a given CCF

    Parameters:
    rv (array): array of radial velocities
    ccf (array): cross-correlation function (CCF)

    Return:
    float: RV value

     """

    # Min of CCF
    imin = int(np.argmin(ccf))

    # Good RV guess
    rv_p0 = rv[imin]

    # Define blue and red wings of the CCF where the derivative changes sign
    width_blue = imin - np.max(np.where(np.gradient(ccf[0:imin])>0))
    width_red = np.min(np.where(np.gradient(ccf[imin:])<0))


    width = int(np.min([width_blue, width_red]))

    if width < 20: # Some cases the width is very short
        width = int(np.max([width_blue, width_red]))
    #ccf = ccf[imin-width_blue:imin+width_red]

    # Normalize CCF between 0 and 1 (continium in 1)
    ccf -= np.min(ccf)
    ccf /= np.min( ccf[ [imin - width,imin + width] ])

    ccf = ccf[imin-width_blue : imin+width_red]
    rv = rv[imin-width_blue : imin+width_red]

    ## Get RV from gaussian profile
    pInit , cov = curve_fit(gauss, rv, ccf,
                            p0=[-(np.max(ccf)-np.min(ccf)),rv_p0,1.0, np.max(ccf)])
    rv_target = pInit[1]
    return rv_target

def gauss(x,amp, x0,sigma, cont):
    """ Obtain inverted Gaussian profile of a given set of parameters """
    y = - amp * np.exp(-0.5*(x-x0)**2/sigma**2) +cont
    return y

def gauss_with_err(x, amp, x0, sigma, cont):
    amp = ufloat(amp[0], amp[1])
    x0 = ufloat(x0[0], x0[1])
    sigma = ufloat(sigma[0], sigma[1])
    cont = ufloat(cont[0], cont[1])
    y = - amp * exp(-0.5*(x-x0)**2/sigma**2) +cont
    return y

def biGauss(x,amp,x0,sigma,cont,A):
    """ Computes BiGaussian profile from Figueira et al. 2013

    Parameters:
    x (array): array of values to compute the function
    amp (float): amplitude of the Gaussian
    x0 (float): center of the Gaussian
    sigma (float): sigma of the Gaussian
    continium (float): continium
    A (float): asymmetry parameter applied in the Gaussian's sigma

    Return:
    array: BiGaussian profile

    """
    return np.piecewise(x,
                        [x<x0, x>=x0],
                         [lambda x: amp * np.exp(-0.5 * ((x-x0)/(sigma*(1-A)))**2) + cont,
                          lambda x: amp * np.exp(-0.5 * ((x-x0)/(sigma*(1+A)))**2) + cont])


def bisector(rv, ccf, low_high_cut = 0.1):
    """ Computes the bisector line of a function of the depth of a given CCF.
        Based on Etienne Artigau's (UdM) version.

    Parameters:
    rv (array): array of radial velocities
    ccf (array): cross-correlation function (CCF)
    low_high_cut (float): limit value in the depth of the CCF

    Returns:
    array: depth
    array: bisector
    """
    # get minima
    imin = int(np.argmin(ccf))

    # get point where the derivative changes sign at the edge of the line
    # blue and red wings of the CCF
    width_blue =  imin - np.max(np.where(np.gradient(ccf[0:imin])>0))
    width_red = np.min(np.where(np.gradient(ccf[imin:])<0))

    # get the width from the side of the center that reaches
    # that point first
    width = int(np.min([width_blue, width_red]))

    # Some cases the width is very short because of low-quality CCFs
    # Let's keep the max instead
    if width < 20:
        width = int(np.max([width_blue, width_red]))

    # set depth to zero
    ccf -= np.min(ccf)

    # set continuum to one
    ccf /= np.min( ccf[ [imin - width,imin + width] ])

    # interpolate each side of the ccf slope at a range of depths
    depth = np.arange(low_high_cut,1-low_high_cut,0.001)

    # blue and red side of line
    g1 = (ccf[imin:imin - width:-1]>low_high_cut) & (ccf[imin:imin - width:-1]<(1-low_high_cut))
    spline1 = InterpolatedUnivariateSpline(ccf[imin:imin - width:-1][g1],rv[imin:imin - width:-1 ][g1], k=2)

    g2 = (ccf[imin : imin + width]>low_high_cut) & (ccf[imin : imin + width]<(1-low_high_cut))
    spline2 = InterpolatedUnivariateSpline(ccf[imin : imin + width][g2],rv[imin : imin + width][g2], k=2)

    # get midpoint
    bisector = (spline2(depth)+spline1(depth))/2

    # get bisector widht
    width_ccf = (spline2(depth)-spline1(depth))

    # define depth in the same way as Perryman, 0 is top, 1 is bottom
    depth = 1-depth

    return depth, bisector