python/c10teff.py at master · astroChasqui/python · 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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
#!/usr/bin/env python
import irtools
import numpy as np
from astropy.io import ascii
import logging
import os


logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger(__name__)

def many(infile, outfile):
    """Calculates Teff using Casagrande et al. (2010) calibrations.

    One column of input file must have 'id' as header
    """
    c10_coef = get_c10coef()

    obs = ascii.read(infile)
    k_colors = list(set(obs.dtype.names) & set(c10_coef['color']))
    kc = ['teff_'+x for x in k_colors]
    more_keys = [prefix+elt for elt in kc for prefix in ('','err_') ]
    res_keys = ['id', 'teff_color', 'err_teff_color']
    res_keys.extend(more_keys)
    res = dict.fromkeys(res_keys, '')
    res['id'] = obs['id'][:]
    n_stars = len(obs['id'])

    for color in k_colors:
        tt, ett = [], []
        for i in range(n_stars):
            if obs[color][i] != '':
                value = float(obs[color][i])

                ebv = float(obs['ebv'][i])

                ebv = 0
                if 'ebv' in obs.dtype.names:
                    if obs['ebv'][i] == '':
                        ebv = 0
                    else:
                        ebv = float(obs['ebv'][i])

                err_value = 0
                if 'err_'+color in obs.dtype.names:
                    if obs['err_'+color][i] == '':
                        err_value = 0
                    else:
                        err_value = float(obs['err_'+color][i])
                if obs['feh'][i] == '':
                    obs['feh'][i] = 0
                feh = float(obs['feh'][i])
                x = one(color, value, feh, c10_coef, ebv=ebv, err_value=err_value)
            else:
                x = None
            if x:
                tt.append(x[0])
                ett.append(x[1])
            else:
                tt.append(None)
                ett.append(None)
        res['teff_'+color] = tt
        res['err_teff_'+color] = ett

    teff = []
    eteff = []
    for i in range(n_stars):
        tt = []
        ett = []
        for color in k_colors:
            if res['teff_'+color][i] > 0:
                tt.append(res['teff_'+color][i])
                ett.append(res['err_teff_'+color][i])
        if len(tt) > 0:
            tv = irtools.wmean(tt, ett)
            teff.append(int(tv[0]))
            if tv[2] != None:
                eteff.append(int(tv[1])+int(tv[2]))
            else:
                eteff.append(int(tv[1]))
        else:
            teff.append(None)
            eteff.append(None)
    res['teff_color'] = teff
    res['err_teff_color'] = eteff

    #get rid of None s
    for k, v in res.iteritems():
        for i in range(len(v)):
          if v[i] is None:
              res[k][i] = ''

    ascii.write(res, outfile, delimiter=',', names=res_keys)


def one(color, value, feh, c10_coef, ebv=0, err_value=0, err_feh=0):
    """Calculates Teff using Casagrande et al. (2010) calibrations.

    Example:
    >>>c10_coef = c10teff.get_c10coef()
    >>>c10teff.one('bv', 0.35, -0.15, c10_coef, ebv=0, err_value=0.01, err_feh=0.05)
    """

    if(color not in c10_coef['color']):
        print('Calibration not available for color given.')
        return None, None

    err_msg = ''
    min_value = c10_coef['min_value'][c10_coef['color'] == color]
    max_value = c10_coef['max_value'][c10_coef['color'] == color]
    min_feh = c10_coef['min_feh'][c10_coef['color'] == color]
    max_feh = c10_coef['max_feh'][c10_coef['color'] == color]

    if ebv > 0:
        k_ebv = c10_coef['k'][c10_coef['color'] == color]
        value -= k_ebv * ebv

    if value < min_value or value > max_value:
        err_msg += 'color value is outside of limits of applicability'
    if feh < min_feh or feh > max_feh:
        err_msg +=  '[Fe/H] is outside of limits of applicability'
    if err_msg != '':
        logger.warning(err_msg)
        return None, None

    a = []
    for i in range(6):
        a.append(c10_coef['a'+str(i)][c10_coef['color'] == color])
    teff = theta(a, value, feh)

    etv = etf = 0
    if err_value > 0:
        tmv = theta(a, value-err_value, feh)
        tpv = theta(a, value+err_value, feh)
        etv = (tmv-tpv)/2
    if err_feh > 0:
        tmf = theta(a, value, feh-err_feh)
        tpf = theta(a, value, feh+err_feh)
        etf = (tmf-tpf)/2
    err_clbr = c10_coef['err_clbr'][c10_coef['color'] == color]
    err_teff = np.sqrt(etv**2+etf**2+err_clbr**2)

    return int(teff), int(err_teff)


def theta(a, value, feh):
    return 5040/(a[0]+a[1]*value+a[2]*value**2+a[3]*value*feh+
                 a[4]*feh+a[5]*feh**2)


def get_c10coef():
    """Silly function to read the C10 coefficients file.

    Instead of running:
    >>>c10_coef = ascii.read('c10teff.csv')
    You do:
    >>>c10_coef = c10teff.get_c10coef()

    In this way you don't have to worry about where the file actually is,
    as long as it is in the same directory where the c10teff.py code lives.
    """
    path = os.path.dirname(os.path.realpath(__file__))
    return ascii.read(os.path.join(path, 'c10teff.csv'))


if __name__ == "__main__":
    import sys
    many(sys.argv[1], sys.argv[2])