-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfunctions_photometry.py
224 lines (189 loc) · 6.86 KB
/
functions_photometry.py
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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
import numpy as np
def fAJHK(Hmag,Hmag_err,G2mag,G2mag_err):
'''
Computes the J-, H-, and K-band photometric extinctions using the RJCE techinque.
See equation 1 in Majewski et al. (2011) and Table 1 in Indebetouw et al. (2005)
Input
Hmag : H-band magnitude
Hmag_err : error in H-band magnitude
G2mag : GLIMPSE channel-2 magnitude
G2mag_err : error in GLIMPSE channel-2 magnitude
Output
AJ : J-band extinction
AJerr : error in J-band extinction
AH : H-band extinction
AHerr : error in H-band extinction
AK : K-band extinction
AKerr : error in K-band extinction
'''
# replace None's with infty's so we don't have to deal with missing numbers
if (isinstance(Hmag,float)==True) or (isinstance(Hmag,int)==True):
if Hmag is None:
Hmag = np.infty
if Hmag_err is None:
Hmag_err = np.infty
if G2mag is None:
G2mag = np.infty
if G2mag_err is None:
G2mag_err = np.infty
elif (isinstance(Hmag,np.ndarray)==True) or (isinstance(Hmag,list)==True):
if Hmag is None:
Hmag = np.ones(shape=Hmag) * np.infty
if Hmag_err is None:
Hmag_err = np.ones(shape=Hmag) * np.infty
if G2mag is None:
G2mag = np.ones(shape=Hmag) * np.infty
if G2mag_err is None:
G2mag_err = np.ones(shape=Hmag) * np.infty
# RJCE techinque
AK = 0.918 * (Hmag - G2mag - 0.08)
AKerr = 0.918 * np.sqrt((Hmag_err)**2. + (G2mag_err)**2.)
# scaling relations from Majewski
AJ = 2.5 * AK
AJerr = 2.5 * AKerr
# scaling relations from Majewski
AH = 1.55 * AK
AHerr = 1.55 * AKerr
return AJ, AJerr, AH, AHerr, AK, AKerr
def fAJHK_0(Jmag,Jmag_err,Hmag,Hmag_err,Kmag,Kmag_err,G2mag,G2mag_err):
'''
Computes the J-, H-, and K-band intrinsic photometric magnitudes using the RJCE techinque.
See equation 1 in Majewski et al. (2011) and Table 1 in Indebetouw et al. (2005)
Input
Jmag : J-band magnitude (optional)
Jmag_err : error in J-band magnitude (optional)
Hmag : H-band magnitude (required by the RJCE techinque)
Hmag_err : error in H-band magnitude (required by the RJCE techinque)
Kmag : K-band magnitude (optional)
Kmag_err : error in K-band magnitude (optional)
G2mag : GLIMPSE channel-2 magnitude (required by the RJCE techinque)
G2mag_err : error in GLIMPSE channel-2 magnitude (required by the RJCE techinque)
Output
J_0 : J-band instinsic magnnintude
J_0err : error in J-band instinsic magnnintude
H_0 : H-band instinsic magnnintude
H_0err : error in H-band instinsic magnnintude
K_0 : K-band instinsic magnnintude
K_0err : error in K-band instinsic magnnintude
'''
# replace None's with infty's so we don't have to deal with missing numbers
if (isinstance(Hmag,float)==True) or (isinstance(Hmag,int)==True):
if Jmag is None:
Jmag = np.infty
if Jmag_err is None:
Jmag_err = np.infty
if Hmag is None:
Hmag = np.infty
if Hmag_err is None:
Hmag_err = np.infty
if Kmag is None:
Kmag = np.infty
if Kmag_err is None:
Kmag_err = np.infty
if G2mag is None:
G2mag = np.infty
if G2mag_err is None:
G2mag_err = np.infty
elif (isinstance(Hmag,np.ndarray)==True) or (isinstance(Hmag,list)==True):
if Jmag is None:
Jmag = np.ones(shape=Hmag) * np.infty
if Jmag_err is None:
Jmag_err = np.ones(shape=Hmag) * np.infty
if Hmag is None:
Hmag = np.ones(shape=Hmag) * np.infty
if Hmag_err is None:
Hmag_err = np.ones(shape=Hmag) * np.infty
if Kmag is None:
Kmag = np.ones(shape=Hmag) * np.infty
if Kmag_err is None:
Kmag_err = np.ones(shape=Hmag) * np.infty
# RJCE techinque (required)
AK = 0.918 * (Hmag - G2mag - 0.08)
AH = 1.55 * AK
AJ = 2.5 * AK
# uncertainty in AK, AJ, AH
AKerr = 0.918 * np.sqrt((Hmag_err)**2. + (G2mag_err)**2.)
AJerr = 2.5 * AKerr
AHerr = 1.55*0.918 * np.sqrt((Hmag_err)**2. + (G2mag_err)**2.)
# intrinsic K-, J-, and H-band magnitudes
H0 = Hmag - AH
J0 = Jmag - AJ
K0 = Kmag - AK
# uncertainty in intrinsic K-, J-, and H-band magnitudes
H0_err = np.sqrt((1.-1.55*0.918)**2.*Hmag_err**2. + (1.55*0.918)**2.*G2mag_err**2.)
J0_err = np.sqrt(Jmag_err**2. + AJerr**2.)
K0_err = np.sqrt(Kmag_err**2. + AKerr**2.)
return J0,J0_err,H0,H0_err,K0,K0_err
def fLJHK(J0,H0,K0,p_mas):
'''
Computes the distance-corrected intrinsic J-, H-, and K-band luminosities usinng Gaia parallaxes,
Input:
J0 : extinction-corrected J-band magintude
H0 : extinction-corrected H-band magintude
K0 : extinction-corrected K-band magintude
p_mas : Gaia parallaxes in milli-arcsec
'''
# replace None's with infty's so we don't have to deal with missing numbers
if (isinstance(p_mas,float)==True) or (isinstance(p_mas,int)==True):
if J0 is None:
J0 = np.infty
if H0 is None:
H0 = np.infty
if K0 is None:
K0 = np.infty
elif (isinstance(p_mas,np.ndarray)==True) or (isinstance(p_mas,list)==True):
if J0 is None:
J0 = np.ones(shape=p_mas.shape) * np.infty
if H0 is None:
H0 = np.ones(shape=p_mas.shape) * np.infty
if K0 is None:
K0 = np.ones(shape=p_mas.shape) * np.infty
# Gaia parallaxes --> distances
p_arcsec,d_pc,d_10pc,d_Mpc = fGaiaParallax(p_mas)
# correct intrinsic magnintudes for Gaia distances
L_J = J0 - 5.*np.log10(d_10pc)
L_H = H0 - 5.*np.log10(d_10pc)
L_K = K0 - 5.*np.log10(d_10pc)
return L_J,L_H,L_K
def fGaiaParallax(p_mas):
'''
Computes distance using Gaia parallax.
Input
p_mas : Gaia parallax in milli-arcseconds
Output
p_as : parallax
d_pc :
d_10pc :
'''
# convert parallax from milli-arcsec to arcsec
p_arcsec = p_mas*1E-3
# compute distances
d_pc = 1./p_arcsec
d_10pc = d_pc/10.
d_Mpc = d_pc/1E6
return p_arcsec,d_pc,d_10pc,d_Mpc
def fmagerrtosnr(magerr):
'''
Converts photometric magnitude uncertainty to signal-to-noise ratio.
'''
magsnr = 1./magerr
return magsnr
def fmagsnrtoerr(magsnr):
'''
Converts photometric magnitude signal-to-noise ratio to uncertainty.
'''
magerr = 1./abs(magsnr)
return magerr
def fconvert_AB_vega(mag_AB,zeropoint):
'''
Converts AB magnitudes to the Vega magnitude scale.
'''
mag_Vega = mag_AB - zeropoint
return mag_Vega
def fcolour_err(mag1,mag1err,mag2,mag2err):
'''
Computes photometric colour and its uncertainty.
'''
colour = mag1 - mag2
colour_err = np.sqrt(mag1err**2. + mag2err**2.)
return colour,colour_err