1
2 """
3 Extinction models of Arenou, Drimmel and Marschall
4
5 Uniformly rewritten by K. Smolders
6 """
7
8 from ivs.catalogs import vizier
9 from ivs.sed.reddening import get_law
10 from ivs.aux.decorators import memoized
11 from ivs.aux import loggers
12 from ivs import config
13
14 import numpy as np
15 from numpy import (abs, arange, array, ceil, cos, dot, floor, int, logical_and,
16 logical_or, max, min, ones, pi, sin, sqrt, where, zeros, exp)
17 import scipy as sc
18 import astropy.io.fits as pf
19 import logging
20
21 logger = logging.getLogger("SED.EXT")
22 logger.addHandler(loggers.NullHandler)
23
24
25
26
27 -def findext(lng, lat, model='drimmel', distance=None, **kwargs):
28 """
29 Get the "model" extinction at a certain longitude and latitude.
30
31 Find the predicted V-band extinction (Av) based on three
32 dimensional models of the galactic interstellar extinction.
33 The user can choose between different models by setting the
34 model keyword:
35
36 1) "arenou": model from Arenou et al. (1992).
37 2) "schlegel": model from Schlegel et al. (1998)
38 3) "drimmel": model from Drimmel et al. (2003)
39 4) "marshall": model from Marshall et al. (2006)
40
41 example useage:
42
43 1. Find the total galactic extinction for a star at galactic lattitude 10.2
44 and longitude 59.0 along the line of sight, as given by the model of
45 Arenou et al. (1992)
46
47 >>> lng = 10.2
48 >>> lat = 59.0
49 >>> av = findext(lng, lat, model='arenou')
50 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av))
51 Av at lng = 10.20, lat = 59.00 is 0.05 magnitude
52
53 2. Find the extinction for a star at galactic lattitude 107.05 and
54 longitude -34.93 and a distance of 144.65 parsecs, as given by the
55 model of Arenou et al. (1992)
56
57 >>> lng = 107.05
58 >>> lat = -34.93
59 >>> dd = 144.65
60 >>> av = findext(lng, lat, distance = dd, model='arenou')
61 >>> print("Av at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, av))
62 Av at lng = 107.05, lat = -34.93 and distance = 144.65 parsecs is 0.15 magnitude
63
64 3. Find the Marschall extinction for a star at galactic longitude 10.2 and
65 latitude 9.0. If the distance is not given, we return the
66 complete extinction along the line of sight (i.e. put the star somewhere
67 out of the galaxy).
68
69 >>> lng = 10.2
70 >>> lat = 9.0
71 >>> av = findext(lng, lat, model='marshall')
72 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av))
73 Av at lng = 10.20, lat = 9.00 is 10.67 magnitude
74
75 4. Find the Marshall extinction for a star at galactic lattitude 271.05 and
76 longitude -4.93 and a distance of 144.65 parsecs, but convert to Av using
77 Rv = 2.5 instead of Rv = 3.1
78
79 >>> lng = 271.05
80 >>> lat = -4.93
81 >>> dd = 144.65
82 >>> av = findext(lng, lat, distance = dd, model='marshall', Rv=2.5)
83 >>> print("Av at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, av))
84 Av at lng = 271.05, lat = -4.93 and distance = 144.65 parsecs is 13.95 magnitude
85
86 5. Find the extinction for a star at galactic longitude 10.2 and
87 latitude 9.0, using the schlegel model, using Rv=2.9 instead of
88 Rv=3.1
89
90 >>> lng = 58.2
91 >>> lat = 24.0
92 >>> distance = 848.
93 >>> av = findext(lng, lat, distance=distance)
94 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av))
95 Av at lng = 58.20, lat = 24.00 is 0.12 magnitude
96
97 REMARKS:
98 a) Schlegel actually returns E(B-V), this value is then converted to Av (the desired value for Rv can be set as a keyword; standard sets Rv=3.1)
99 b) Schlegel is very dubious for latitudes between -5 and 5 degrees
100 c) Marschall actually returns Ak, this value is then converted to Av (the reddening law and Rv can be set as keyword; standard sets Rv=3.1, redlaw='cardelli1989')
101 d) Marschall is only available for certain longitudes and latitudes:
102 0 < lng < 100 or 260 < lng < 360 and -10 < lat < 10
103
104 @param lng: Galactic Longitude (in degrees)
105 @type lng: float
106 @param lat: Galactic Lattitude (in degrees)
107 @type lat: float
108 @param model: the name of the extinction model: ("arenou", "schlegel", "drimmel" or "marshall"; if none given, the program uses "drimmel")
109 @type model: str
110 @param distance: Distance to the source (in parsecs), if the distance is not given, the total galactic extinction along the line of sight is calculated
111 @type distance: float
112 @return: The extinction in Johnson V-band
113 @rtype: float
114 """
115
116 if model.lower() == 'drimmel':
117 av = findext_drimmel(lng, lat, distance=distance, **kwargs)
118 elif model.lower() == 'marshall' or model.lower() == 'marschall':
119 av = findext_marshall(lng, lat, distance=distance, **kwargs)
120 elif model.lower() == 'arenou':
121 av = findext_arenou(lng, lat, distance=distance, **kwargs)
122 elif model.lower() == 'schlegel':
123 av = findext_schlegel(lng, lat, distance=distance, **kwargs)
124 return(av)
125
126
127
128
129
130
131
132
133
134 -def findext_arenou(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='Av'):
135 """
136 Find the predicted V-band extinction (Av) according to the
137 3D model for galactic extinction of Arenou et al, "Modelling
138 the Galactic interstellar extinction distribution in three
139 dimensions", Arenou et al, "A tridimensional model of the
140 galactic interstellar extinction" published in Astronomy and
141 Astrophysics (ISSN 0004-6361), vol. 258, no. 1, p. 104-111, 1992
142
143 Example usage:
144
145 1. Find the extinction for a star at galactic lattitude 10.2 and
146 longitude 59.0. If the distance is not given, we use the
147 maximal distance r0 for that line of sight, as given in the
148 Appendix of Arenou et al. (1992)
149
150 >>> lng = 10.2
151 >>> lat = 59.0
152 >>> av = findext_arenou(lng, lat)
153 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av))
154 Av at lng = 10.20, lat = 59.00 is 0.05 magnitude
155
156 2. Find the extinction for a star at galactic lattitude 107.05 and
157 longitude -34.93 and a distance of 144.65 parsecs
158
159 >>> lng = 107.05
160 >>> lat = -34.93
161 >>> dd = 144.65
162 >>> av = findext_arenou(lng, lat, distance = dd)
163 >>> print("Av at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, av))
164 Av at lng = 107.05, lat = -34.93 and distance = 144.65 parsecs is 0.15 magnitude
165
166 @param ll: Galactic Longitude (in degrees)
167 @type ll: float
168 @param bb: Galactic Lattitude (in degrees)
169 @type bb: float
170 @param distance: Distance to the source (in parsecs)
171 @type distance: float
172 @return: The extinction in Johnson V-band
173 @rtype: float
174 """
175
176 if bb < -90. or bb > 90:
177 logger.error("galactic lattitude outside [-90,90] degrees")
178 elif ll < 0. or ll > 360:
179 logger.error("galactic longitude outside [0,360] degrees")
180 elif distance < 0 and distance is not None:
181 logger.error("distance is negative")
182
183
184 alpha, beta, gamma, rr0, saa = _getarenouparams(ll, bb)
185 logger.info("Arenou params: alpha = %.2f, beta = %.2f, gamma = %.2f, r0 = %.2f and saa = %.2f" %(alpha, beta, gamma, rr0, saa))
186
187
188
189 if distance is None:
190 av = alpha*rr0 + beta*rr0**2.
191 else:
192 distance = distance/1e3
193 if distance <= rr0:
194 av = alpha*distance + beta*distance**2.
195 else:
196 av = alpha*rr0 + beta*rr0**2. + (distance-rr0)*gamma
197
198
199 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'])
200
201 return av/redflux[0]
202
204 """
205 Input: galactic coordinates
206 Output: Arenou 1992 alpha, beta, gamma, rr0, saa
207
208 @param ll: Galactic Longitude (in degrees)
209 @type ll: float
210 @param bb: Galactic Lattitude (in degrees)
211 @type bb: float
212 """
213 if -90 < bb < -60:
214 gamma = 0
215 if 0 <= ll < 29:
216 alpha = 2.22538 ; beta = -6.00212 ; rr0 = 0.052 ; saa = 13.
217 elif 29 <= ll < 57:
218 alpha = 3.35436 ; beta = -14.74567 ; rr0 = 0.035 ; saa = 40
219 elif 57 <= ll < 85:
220 alpha = 2.77637 ; beta = -9.62706 ; rr0 = 0.042 ; saa = 15
221 elif 85 <= ll < 110:
222 alpha = 4.44190 ; beta = -19.92097 ; rr0 = 0.025 ; saa = 36
223 elif 110 <= ll < 150:
224 alpha = 4.46685 ; beta = -26.07305 ; rr0 = 0.026 ; saa = 28
225 elif 150 <= ll < 180:
226 alpha = 7.63699 ; beta = -46.10856 ; rr0 = 0.014 ; saa = 38
227 elif 180 <= ll < 210:
228 alpha = 2.43412 ; beta = -8.69913 ; rr0 = 0.050 ; saa = 36
229 elif 210 <= ll < 240:
230 alpha = 3.34481 ; beta = -13.93228 ; rr0 = 0.035 ; saa = 38
231 elif 240 <= ll < 270:
232 alpha = 1.40733 ; beta = -13.43418 ; rr0 = 0.091 ; saa = 30
233 elif 270 <= ll < 300:
234 alpha = 1.64466 ; beta = -3.97380 ; rr0 = 0.074 ; saa = 28
235 elif 300 <= ll < 330:
236 alpha = 2.12696 ; beta = -6.05682 ; rr0 = 0.056 ; saa = 14
237 elif 330 <= ll < 360:
238 alpha = 2.34636 ; beta = -8.17784 ; rr0 = 0.052 ; saa = 16
239
240 if -60 < bb < -45:
241 gamma = 0
242 if 0 <= ll < 30:
243 alpha = 2.77060 ; beta = -9.52310 ; rr0 = 0.145 ; saa = 16
244 elif 30 <= ll < 60:
245 alpha = 1.96533 ; beta = -9.52310 ; rr0 = 0.174 ; saa = 06
246 elif 60 <= ll < 110:
247 alpha = 1.93622 ; beta = -13.31757 ; rr0 = 0.073 ; saa = 26
248 elif 110 <= ll < 180:
249 alpha = 1.05414 ; beta = -2.236540 ; rr0 = 0.223 ; saa = 74
250 elif 180 <= ll < 210:
251 alpha = 1.39990 ; beta = -1.35325 ; rr0 = 0.252 ; saa = 10
252 elif 210 <= ll < 240:
253 alpha = 2,73481 ; beta = -11.70266 ; rr0 = 0.117 ; saa = 8
254 elif 240 <= ll < 270:
255 alpha = 2.99784 ; beta = -11.64272 ; rr0 = 0.129 ; saa = 3
256 elif 270 <= ll < 300:
257 alpha = 3.23802 ; beta = -11.63810 ; rr0 = 0.139 ; saa = 7
258 elif 300 <= ll < 330:
259 alpha = 1.72679 ; beta = -6.05085 ; rr0 = 0.143 ; saa = 7
260 elif 330 <= ll < 360:
261 alpha = 1.88890 ; beta = -5.51861 ; rr0 = 0.171 ; saa = 14
262
263 if -45 < bb < -30:
264 gamma = 0
265 if 0 <= ll < 30:
266 alpha = 1.98973 ; beta = -4.86206 ; rr0 = 0.205 ; saa = 6
267 elif 30 <= ll < 60:
268 alpha = 1.49901 ; beta = -3.75837 ; rr0 = 0.199 ; saa = 16
269 elif 60 <= ll < 90:
270 alpha = 0.90091 ; beta = -1.30459 ; rr0 = 0.329 ; saa = 73
271 elif 90 <= ll < 120:
272 alpha = 1.94200 ; beta = -6.26833 ; rr0 = 0.155 ; saa = 18
273 elif 120 <= ll < 160:
274 alpha = -0.37804 ; beta = 10.77372 ; rr0 = 0.210 ; saa = 100
275 elif 160 <= ll < 200:
276 alpha = -0.15710 ; beta = 5.03190 ; rr0 = 0.294 ; saa = 24
277 elif 200 <= ll < 235:
278 alpha = 3.20162 ; beta = -10.59297 ; rr0 = 0.151 ; saa = 9
279 elif 235 <= ll < 265:
280 alpha = 1.95079 ; beta = 4.73280 ; rr0 = 0.206 ; saa = 21
281 elif 265 <= ll < 300:
282 alpha = 1.91200 ; beta = 4.97640 ; rr0 = 0.192 ; saa = 17
283 elif 300 <= ll < 330:
284 alpha = 2.50487 ; beta = -9.63106 ; rr0 = 0.145 ; saa = 28
285 elif 330 <= ll < 360:
286 alpha = 2.44394 ; beta = -9.17612 ; rr0 = 0.133 ; saa = 7
287
288 if -30 < bb < -15:
289 gamma = 0
290 if 0 <= ll < 20:
291 alpha = 2.82440 ; beta = -4.78217 ; rr0 = 0.295 ; saa = 32
292 elif 20 <= ll < 40:
293 alpha = 3.84362 ; beta = -8.04690 ; rr0 = 0.239 ; saa = 46
294 elif 40 <= ll < 80:
295 alpha = 0.60365 ; beta = 0.07893 ; rr0 = 0.523 ; saa = 22
296 elif 80 <= ll < 100:
297 alpha = 0.58307 ; beta = -0.21053 ; rr0 = 0.523 ; saa = 53
298 elif 100 <= ll < 120:
299 alpha = 2.03861 ; beta = -4.40843 ; rr0 = 0.231 ; saa = 60
300 elif 120 <= ll < 140:
301 alpha = 1.14271 ; beta = -1.35635 ; rr0 = 0.421 ; saa = 34
302 elif 140 <= ll < 160:
303 alpha = 0.79908 ; beta = 1.48074 ; rr0 = 0.513 ; saa = 61
304 elif 160 <= ll < 180:
305 alpha = 0.94260 ; beta = 8.16346 ; rr0 = 0.441 ; saa = 42
306 elif 180 <= ll < 200:
307 alpha = 1.66398 ; beta = 0.26775 ; rr0 = 0.523 ; saa = 42
308 elif 200 <= ll < 220:
309 alpha = 1.08760 ; beta = -1.02443 ; rr0 = 0.523 ; saa = 45
310 elif 220 <= ll < 240:
311 alpha = 1.20087 ; beta = -2.45407 ; rr0 = 0.245 ; saa = 6
312 elif 240 <= ll < 260:
313 alpha = 1.13147 ; beta = -1.87916 ; rr0 = 0.301 ; saa = 16
314 elif 260 <= ll < 280:
315 alpha = 0.97804 ; beta = -2.92838 ; rr0 = 0.338 ; saa = 21
316 elif 290 <= ll < 300:
317 alpha = 1.40086 ; beta = -1.12403 ; rr0 = 0.523 ; saa = 19
318 elif 300 <= ll < 320:
319 alpha = 2.06355 ; beta = -3.68278 ; rr0 = 0.280 ; saa = 42
320 elif 320 <= ll < 340:
321 alpha = 1.59260 ; beta = -2.18754 ; rr0 = 0.364 ; saa = 15
322 elif 310 <= ll < 360:
323 alpha = 1.45589 ; beta = -1.90598 ; rr0 = 0.382 ; saa = 21
324
325 if -15 < bb < -5:
326 gamma = 0
327 if 0 <= ll < 10:
328 alpha = 2.56438 ; beta = -2.31586 ; rr0 = 0.554 ; saa = 37
329 elif 10 <= ll < 20:
330 alpha = 3.24095 ; beta = -2.78217 ; rr0 = 0.582 ; saa = 38
331 elif 20 <= ll < 30:
332 alpha = 2.95627 ; beta = -2.57422 ; rr0 = 0.574 ; saa = 41 ; gamma = 0.08336
333 elif 30 <= ll < 40:
334 alpha = 1.85158 ; beta = -0.67716 ; rr0 = 1.152 ; saa = 4
335 elif 40 <= ll < 50:
336 alpha = 1.60770 ; beta = 0.35279 ; rr0 = 0.661 ; saa = 24
337 elif 50 <= ll < 60:
338 alpha = 0.69920 ; beta = -0.09146 ; rr0 = 0.952 ; saa = 2 ; gamma = 0.12839
339 elif 60 <= ll < 80:
340 alpha = 1.36189 ; beta = -1.05290 ; rr0 = 0.647 ; saa = 45 ; gamma = 0.16258
341 elif 80 <= ll < 90:
342 alpha = 0.33179 ; beta = 0.37629 ; rr0 = 1.152 ; saa = 62
343 elif 90 <= ll < 100:
344 alpha = 1.70303 ; beta = -0.75246 ; rr0 = 1.132 ; saa = 31
345 elif 100 <= ll < 110:
346 alpha = 1.97414 ; beta = -1.59784 ; rr0 = 0.618 ; saa = 35 ; gamma = 0.12847
347 elif 110 <= ll < 120:
348 alpha = 1.07407 ; beta = -0.40066 ; rr0 = 1.152 ; saa = 14 ; gamma = 0.17698
349 elif 120 <= ll < 130:
350 alpha = 1.69495 ; beta = -1.00071 ; rr0 = 0.847 ; saa = 28 ; gamma = 0.08567
351 elif 130 <= ll < 140:
352 alpha = 1.51449 ; beta = -0.08441 ; rr0 = 0.897 ; saa = 12
353 elif 140 <= ll < 150:
354 alpha = 1.87894 ; beta = -0.73314 ; rr0 = 1.152 ; saa = 23
355 elif 150 <= ll < 160:
356 alpha = 1.43670 ; beta = 0.67706 ; rr0 = 0.778 ; saa = 3 ; gamma = 0.42624
357 elif 160 <= ll < 180:
358 alpha = 6.84802 ; beta = -5.06864 ; rr0 = 0.676 ; saa = 50
359 elif 180 <= ll < 190:
360 alpha = 4.16321 ; beta = -5.80016 ; rr0 = 0.359 ; saa = 51 ; gamma = 0.60387
361 elif 190 <= ll < 200:
362 alpha = 0.78135 ; beta = -0.27826 ; rr0 = 1.152 ; saa = 4
363 elif 200 <= ll < 210:
364 alpha = 0.85535 ; beta = 0.20848 ; rr0 = 1.152 ; saa = 17
365 elif 210 <= ll < 220:
366 alpha = 0.52521 ; beta = 0.65726 ; rr0 = 1.152 ; saa = 7
367 elif 220 <= ll < 230:
368 alpha = 0.88376 ; beta = -0.44519 ; rr0 = 0.993 ; saa = 40 ; gamma = 0.06013
369 elif 230 <= ll < 240:
370 alpha = 0.42228 ; beta = 0.26304 ; rr0 = 0.803 ; saa = 26
371 elif 240 <= ll < 250:
372 alpha = 0.71318 ; beta = -0.67229 ; rr0 = 0.530 ; saa = 55 ; gamma = 0.20994
373 elif 250 <= ll < 260:
374 alpha = 0.99606 ; beta = -0.70103 ; rr0 = 0.710 ; saa = 48 ; gamma = 0.01323
375 elif 260 <= ll < 270:
376 alpha = 0.91519 ; beta = -0.39690 ; rr0 = 1.152 ; saa = 48 ; gamma = 0.01961
377 elif 270 <= ll < 280:
378 alpha = 0.85791 ; beta = -0.29115 ; rr0 = 1.152 ; saa = 19
379 elif 280 <= ll < 290:
380 alpha = 1.44226 ; beta = -1.09775 ; rr0 = 0.657 ; saa = 39
381 elif 290 <= ll < 300:
382 alpha = 2.55486 ; beta = -1.68293 ; rr0 = 0.759 ; saa = 31
383 elif 300 <= ll < 310:
384 alpha = 3.18047 ; beta = -2.69796 ; rr0 = 0.589 ; saa = 40
385 elif 210 <= ll < 320:
386 alpha = 2.11235 ; beta = -1.77506 ; rr0 = 0.595 ; saa = 29
387 elif 320 <= ll < 330:
388 alpha = 1.75884 ; beta = -1.38574 ; rr0 = 0.635 ; saa = 25
389 elif 330 <= ll < 340:
390 alpha = 1.97257 ; beta = -1.55545 ; rr0 = 0.634 ; saa = 34 ; gamma = 0.00043
391 elif 340 <= ll < 350:
392 alpha = 1.41497 ; beta = -1.05722 ; rr0 = 0.669 ; saa = 46 ; gamma = 0.03264
393 elif 350 <= ll < 360:
394 alpha = 1.17795 ; beta = -0.95012 ; rr0 = 0.620 ; saa = 46 ; gamma = 0.03339
395
396 if -5 < bb < 5:
397 gamma = 0
398 if 0 <= ll < 10:
399 alpha = 2.62556 ; beta = -1.11097 ; rr0 = 1.182 ; saa = 57 ; gamma = 0.00340
400 elif 10 <= ll < 20:
401 alpha = 3.14461 ; beta = -1.01140 ; rr0 = 1.555 ; saa = 42
402 elif 20 <= ll < 30:
403 alpha = 4.26624 ; beta = -1.61242 ; rr0 = 1.323 ; saa = 34
404 elif 30 <= ll < 40:
405 alpha = 2.54447 ; beta = -0.12771 ; rr0 = 1.300 ; saa = 30
406 elif 40 <= ll < 50:
407 alpha = 2.27030 ; beta = -0.68720 ; rr0 = 1.652 ; saa = 52 ; gamma = 0.04928
408 elif 50 <= ll < 60:
409 alpha = 1.34359 ; beta = -0.05416 ; rr0 = 2.000 ; saa = 32
410 elif 60 <= ll < 70:
411 alpha = 1.76327 ; beta = -0.26407 ; rr0 = 2.000 ; saa = 37
412 elif 70 <= ll < 80:
413 alpha = 2.20666 ; beta = -0.41651 ; rr0 = 2.000 ; saa = 36
414 elif 80 <= ll < 90:
415 alpha = 1.50130 ; beta = -0.09855 ; rr0 = 1.475 ; saa = 45
416 elif 90 <= ll < 100:
417 alpha = 2.43965 ; beta = -0.82128 ; rr0 = 1.485 ; saa = 36 ; gamma = 0.01959
418 elif 100 <= ll < 110:
419 alpha = 3.35775 ; beta = -1.16400 ; rr0 = 0.841 ; saa = 35 ; gamma = 0.00298
420 elif 110 <= ll < 120:
421 alpha = 2.60621 ; beta = -0.68687 ; rr0 = 1.897 ; saa = 36
422 elif 120 <= ll < 130:
423 alpha = 2.90112 ; beta = -0.97988 ; rr0 = 1.480 ; saa = 32
424 elif 130 <= ll < 140:
425 alpha = 2.55377 ; beta = -0.71214 ; rr0 = 1.793 ; saa = 38
426 elif 140 <= ll < 150:
427 alpha = 3.12598 ; beta = -1.21437 ; rr0 = 1.287 ; saa = 23 ; gamma = 0.15298
428 elif 150 <= ll < 160:
429 alpha = 3.66930 ; beta = -2.29731 ; rr0 = 0.799 ; saa = 40 ; gamma = 0.33473
430 elif 160 <= ll < 170:
431 alpha = 2.15465 ; beta = -0.70690 ; rr0 = 1.524 ; saa = 37 ; gamma = 0.14017
432 elif 170 <= ll < 180:
433 alpha = 1.82465 ; beta = -0.60223 ; rr0 = 1.515 ; saa = 29 ; gamma = 0.20730
434 elif 180 <= ll < 190:
435 alpha = 1.76269 ; beta = -0.35945 ; rr0 = 2.000 ; saa = 28 ; gamma = 0.08052
436 elif 190 <= ll < 200:
437 alpha = 1.06085 ; beta = -0.14211 ; rr0 = 2.000 ; saa = 48
438 elif 200 <= ll < 210:
439 alpha = 1.21333 ; beta = -0.23225 ; rr0 = 2.000 ; saa = 57
440 elif 210 <= ll < 220:
441 alpha = 0.58326 ; beta = -0.06097 ; rr0 = 2.000 ; saa = 30 ; gamma = 0.36962
442 elif 220 <= ll < 230:
443 alpha = 0.74200 ; beta = -0.19293 ; rr0 = 1.923 ; saa = 48 ; gamma = 0.07459
444 elif 230 <= ll < 240:
445 alpha = 0.67520 ; beta = -0.21041 ; rr0 = 1.604 ; saa = 49 ; gamma = 0.16602
446 elif 240 <= ll < 250:
447 alpha = 0.62609 ; beta = -0.25312 ; rr0 = 1.237 ; saa = 73 ; gamma = 0.14437
448 elif 250 <= ll < 260:
449 alpha = 0.61415 ; beta = -0.13788 ; rr0 = 2.000 ; saa = 42 ; gamma = 0.26859
450 elif 260 <= ll < 270:
451 alpha = 0.58108 ; beta = 0.01195 ; rr0 = 2.000 ; saa = 40 ; gamma = 0.07661
452 elif 270 <= ll < 280:
453 alpha = 0.68352 ; beta = -0.10743 ; rr0 = 2.000 ; saa = 50 ; gamma = 0.00849
454 elif 280 <= ll < 290:
455 alpha = 0.61747 ; beta = 0.02675 ; rr0 = 2,000 ; saa = 49
456 elif 290 <= ll < 300:
457 alpha = 0.06827 ; beta = -0.26290 ; rr0 = 2.000 ; saa = 44
458 elif 300 <= ll < 310:
459 alpha = 1.53631 ; beta = -0.36833 ; rr0 = 2.000 ; saa = 37 ; gamma = 0.02960
460 elif 210 <= ll < 320:
461 alpha = 1.94300 ; beta = -0.71445 ; rr0 = 1.360 ; saa = 36 ; gamma = 0.15643
462 elif 320 <= ll < 330:
463 alpha = 1.22185 ; beta = -0.40185 ; rr0 = 1.520 ; saa = 48 ; gamma = 0.07354
464 elif 330 <= ll < 340:
465 alpha = 1.79429 ; beta = -0.48657 ; rr0 = 1.844 ; saa = 43
466 elif 340 <= ll < 350:
467 alpha = 2.29545 ; beta = -0.84096 ; rr0 = 1.365 ; saa = 32
468 elif 350 <= ll < 360:
469 alpha = 2.07408 ; beta = -0.64745 ; rr0 = 1.602 ; saa = 36 ; gamma = 0.12750
470
471
472 if 5 < bb < 15:
473 gamma = 0
474 if 0 <= ll < 10:
475 alpha = 2.94213 ; beta = -2.09158 ; rr0 = 0.703 ; saa = 41 ; gamma = 0.05490
476 elif 10 <= ll < 30:
477 alpha = 3.04627 ; beta = 7.71159 ; rr0 = 0.355 ; saa = 45
478 elif 30 <= ll < 40:
479 alpha = 3.78033 ; beta = -3.91956 ; rr0 = 0.482 ; saa = 42
480 elif 40 <= ll < 50:
481 alpha = 2.18119 ; beta = -2.4050 ; rr0 = 0.453 ; saa = 27
482 elif 50 <= ll < 60:
483 alpha = 1.45372 ; beta = -0.49522 ; rr0 = 1.152 ; saa = 31
484 elif 60 <= ll < 70:
485 alpha = 1.05051 ; beta = -1.01704 ; rr0 = 0.516 ; saa = 2
486 elif 70 <= ll < 80:
487 alpha = 0.48416 ; beta = -0.27182 ; rr0 = 0.891 ; saa = 94 ; gamma = 0.08639
488 elif 80 <= ll < 90:
489 alpha = 0.61963 ; beta = 0.41697 ; rr0 = 1.152 ; saa = 35 ; gamma = 0.47171
490 elif 90 <= ll < 100:
491 alpha = 4.40348 ; beta = -2.95011 ; rr0 = 0.745 ; saa = 52
492 elif 100 <= ll < 120:
493 alpha = 2.50938 ; beta = -0.56541 ; rr0 = 1.152 ; saa = 27
494 elif 120 <= ll < 130:
495 alpha = 0.44180 ; beta = 1.58923 ; rr0 = 0.949 ; saa = 4
496 elif 130 <= ll < 140:
497 alpha = 3.96081 ; beta = -3.37374 ; rr0 = 0.587 ; saa = 40 ; gamma = 0.34109
498 elif 140 <= ll < 160:
499 alpha = 2.53335 ; beta = -0.40541 ; rr0 = 1.152 ; saa = 38
500 elif 160 <= ll < 170:
501 alpha = 2.03760 ; beta = -0.66136 ; rr0 = 1.152 ; saa = 23
502 elif 170 <= ll < 200:
503 alpha = 1.06946 ; beta = -0.87395 ; rr0 = 0.612 ; saa = 29 ; gamma = 0.29230
504 elif 200 <= ll < 210:
505 alpha = 0.86348 ; beta = -0.65870 ; rr0 = 0.655 ; saa = 79 ; gamma = 0.09089
506 elif 210 <= ll < 230:
507 alpha = 0.30117 ; beta = -0.16136 ; rr0 = 0.933 ; saa = 17 ; gamma = 0.07495
508 elif 230 <= ll < 240:
509 alpha = 0.75171 ; beta = -0.57143 ; rr0 = 0.658 ; saa = 12 ; gamma = 0.00534
510 elif 240 <= ll < 250:
511 alpha = 1.97427 ; beta = -2.02654 ; rr0 = 0.487 ; saa = 67
512 elif 250 <= ll < 260:
513 alpha = 1.25208 ; beta = -1.47763 ; rr0 = 0.424 ; saa = 19 ; gamma = 0.09089
514 elif 260 <= ll < 270:
515 alpha = 0.89448 ; beta = -0.43870 ; rr0 = 1.019 ; saa = 5
516 elif 270 <= ll < 280:
517 alpha = 0.81141 ; beta = -0.51001 ; rr0 = 0.795 ; saa = 27 ; gamma = 0.03505
518 elif 280 <= ll < 290:
519 alpha = 0.83781 ; beta = -0.44138 ; rr0 = 0.949 ; saa = 50 ; gamma = 0.02820
520 elif 290 <= ll < 300:
521 alpha = 1.10600 ; beta = -0.86263 ; rr0 = 0.641 ; saa = 28 ; gamma = 0.03402
522 elif 300 <= ll < 310:
523 alpha = 1.37040 ; beta = -1.02779 ; rr0 = 0.667 ; saa = 28 ; gamma = 0.05608
524 elif 310 <= ll < 320:
525 alpha = 1.77590 ; beta = -1.26951 ; rr0 = 0.699 ; saa = 37 ; gamma = 0.06972
526 elif 320 <= ll < 330:
527 alpha = 1.20865 ; beta = -0.70679 ; rr0 = 0.855 ; saa = 35 ; gamma = 0.02902
528 elif 330 <= ll < 340:
529 alpha = 2.28830 ; beta = -1.71890 ; rr0 = 0.666 ; saa = 42 ; gamma = 0.22887
530 elif 340 <= ll < 350:
531 alpha = 3.26278 ; beta = -0.94181 ; rr0 = 1.152 ; saa = 38
532 elif 350 <= ll < 360:
533 alpha = 2.58100 ; beta = -1.69237 ; rr0 = 0.763 ; saa = 53
534
535 if 15 < bb < 30:
536 gamma = 0
537 if 0 <= ll < 20:
538 alpha = 6.23279 ; beta = -10.30384 ; rr0 = 0.302 ; saa = 42
539 elif 20 <= ll < 40:
540 alpha = -4.47693 ; beta = -7.28366 ; rr0 = 0.307 ; saa = 29
541 elif 40 <= ll < 60 :
542 alpha = 1.22938 ; beta = -1.19030 ; rr0 = 0.516 ; saa = 5
543 elif 60 <= ll < 80 :
544 alpha = 0.84291 ; beta = -1.59338 ; rr0 = 0.265 ; saa = 4
545 elif 80 <= ll < 100 :
546 alpha = 0.23996 ; beta = 0.06304 ; rr0 = 0.523 ; saa = 32
547 elif 100 <= ll < 140 :
548 alpha = 0.40062 ; beta = -1.75628 ; rr0 = 0.114 ; saa = 16
549 elif 140 <= ll < 180 :
550 alpha = 0.56898 ; beta = -0.53331 ; rr0 = 0.523 ; saa = 41
551 elif 180 <= ll < 200 :
552 alpha = -0.95721 ; beta = 11.6917 ; rr0 = 0.240 ; saa = 2
553 elif 200 <= ll < 220 :
554 alpha = -0.19051 ; beta = 1.45670 ; rr0 = 0.376 ; saa = 1
555 elif 220 <= ll < 240 :
556 alpha = 2.31305 ; beta = -7.82531 ; rr0 = 0.148 ; saa = 95
557 elif 240 <= ll < 260:
558 alpha = 1.39169 ; beta = -1.72984 ; rr0 = 0.402 ; saa = 6
559 elif 260 <= ll < 260:
560 alpha = 1.59418 ; beta = -1.28296 ; rr0 = 0.523 ; saa = 36
561 elif 280 <= ll < 300 :
562 alpha = 1.57082 ; beta = -197295 ; rr0 = 0.398 ; saa = 10
563 elif 300 <= ll < 320 :
564 alpha = 1.95998 ; beta = -3.26159 ; rr0 = 0.300 ; saa = 11
565 elif 320 <= ll < 340:
566 alpha = 2.59567 ; beta = -4.84133 ; rr0 = 0.168 ; saa = 37
567 elif 340 <= ll < 360:
568 alpha = 5.30273 ; beta = -7.43033 ; rr0 = 0.357 ; saa = 37
569
570 if 30 < bb < 45:
571 gamma = 0
572 if 0 <= ll < 20:
573 alpha = 2.93960 ; beta = -6.48019 ; rr0 = 0.227 ; saa = 77
574 elif 20 <= ll < 50:
575 alpha = 1.65864 ; beta = -9.99317 ; rr0 = 0.083 ; saa = 99
576 elif 50 <= ll < 80:
577 alpha = 1.71831 ; beta = -7.25286 ; rr0 = 0.118 ; saa = 28
578 elif 80 <= ll < 110:
579 alpha = 1.33617 ; beta = -10.39799 ; rr0 = 0.064 ; saa = 99
580 elif 110 <= ll < 160:
581 alpha = -0.31330 ; beta = 1.35622 ; rr0 = 0.329 ; saa = 24
582 elif 160 <= ll < 190:
583 alpha = 1.51984 ; beta = -8.69502 ; rr0 = 0.087 ; saa = 99
584 elif 190 <= ll < 220:
585 alpha = -0.50758 ; beta = 4.73320 ; rr0 = 0.250 ; saa = 78
586 elif 220 <= ll < 250:
587 alpha = 1.25864 ; beta = -12.59627 ; rr0 = 0.050 ; saa = 70
588 elif 250 <= ll < 280:
589 alpha = 1.54243 ; beta = -3.75065 ; rr0 = 0.205 ; saa = 10
590 elif 280 <= ll < 320:
591 alpha = 2.72258 ; beta = -7.47806 ; rr0 = 0.182 ; saa = 5
592 elif 320 <= ll < 340:
593 alpha = 2.81435 ; beta = -5.52139 ; rr0 = 0.255 ; saa = 10
594 elif 340 <= ll < 360:
595 alpha = 2.23818 ; beta = 0.81772 ; rr0 = 0.329 ; saa = 19
596
597 if 45 < bb < 60:
598 gamma = 0
599 if 0 <= ll < 60:
600 alpha = 1.38587 ; beta = -9.06536 ; rr0 = 0.076 ; saa = 3
601 elif 60 <= ll < 90:
602 alpha = 2.28570 ; beta = -9.88812 ; rr0 = 0.116 ; saa = 3
603 elif 90 <= ll < 110:
604 alpha = 1.36385 ; beta = -8.10127 ; rr0 = 0.084 ; saa = 4
605 elif 110 <= ll < 170:
606 alpha = 0.05943 ; beta = -1.08126 ; rr0 = 0.027 ; saa = 50
607 elif 170 <= ll < 200:
608 alpha = 1.40171 ; beta = -3.21783 ; rr0 = 0.218 ; saa = 99
609 elif 200 <= ll < 230:
610 alpha = 0.14718 ; beta = 3.92670 ; rr0 = 0.252 ; saa = 14
611 elif 230 <= ll < 290:
612 alpha = 0.57124 ; beta = -4.30242 ; rr0 = 0.066 ; saa = 10
613 elif 290 <= ll < 330:
614 alpha = 3.69891 ; beta = -19.62204 ; rr0 = 0.094 ; saa = 5
615 elif 330 <= ll < 360:
616 alpha = 1.19563 ; beta = -0.45043 ; rr0 = 0.252 ; saa = 9
617
618 if 60 < bb < 90:
619 gamma = 0
620 if 0 <= ll < 30:
621 alpha = 0.69443 ; beta = -0.27600 ; rr0 = 0.153 ; saa = 99
622 elif 30 <= ll < 60:
623 alpha = 1.11811 ; beta = 0.71179 ; rr0 = 0.085 ; saa = 73
624 elif 60 <= ll < 90:
625 alpha = 1.10427 ; beta = -2.37654 ; rr0 = 0.123 ; saa = 99
626 elif 90 <= ll < 120:
627 alpha = -0.42211 ; beta = 5.24037 ; rr0 = 0.184 ; saa = 12
628 elif 120 <= ll < 150:
629 alpha = 0.87576 ; beta = -4.38033 ; rr0 = 0.100 ; saa = 35
630 elif 150 <= ll < 180:
631 alpha = 1.27477 ; beta = -4.98307 ; rr0 = 0.128 ; saa = 72
632 elif 180 <= ll < 210:
633 alpha = 1.19512 ; beta = -6.58464 ; rr0 = 0.091 ; saa = 49
634 elif 210 <= ll < 240:
635 alpha = 0.97581 ; beta = -4.89869 ; rr0 = 0.100 ; saa = 95
636 elif 240 <= ll < 270:
637 alpha = 0.54379 ; beta = -0.84403 ; rr0 = 0.207 ; saa = 35
638 elif 270 <= ll < 300:
639 alpha = 0.85054 ; beta = 13.01249 ; rr0 = 0.126 ; saa = 39
640 elif 300 <= ll < 330:
641 alpha = 0.74347 ; beta = 1.39825 ; rr0 = 0.207 ; saa = 10
642 elif 330 <= ll < 360:
643 alpha = 0.77310 ; beta = -4.45005 ; rr0 = 0.087 ; saa = 16
644
645 return alpha, beta, gamma, rr0, saa
646
656 """
657 Read in the Marshall data
658 """
659
660
661 filen = config.get_datafile('catalogs','extinction_marshall.tsv')
662 data_ma, units_ma, comments_ma = vizier.tsv2recarray(filen)
663 return data_ma, units_ma, comments_ma
664
665 -def findext_marshall(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):
666 """
667 Find the V-band extinction according to the reddening model of
668 Marshall et al. (2006) published in Astronomy and Astrophysics,
669 Volume 453, Issue 2, July II 2006, pp.635-651
670
671 The band in which the extinction is calculated is actually optional, and given
672 with the keyword C{norm}. If you set C{norm='Av'}, you get visual extinction
673 (in JOHNSON.V) band, if you set C{norm='Ak'}, you get infrared extinction
674 (in JOHNSON.K). If you want the extinction in different passbands, you can
675 set them via C{norm='GENEVA.V'}. So, C{norm='Av'} is equivalent to setting
676 C{norm='JOHNSON.V'}.
677
678 Example usage:
679
680 1. Find the extinction for a star at galactic longitude 10.2 and
681 latitude 9.0. If the distance is not given, we return the
682 complete extinction along the line of sight (i.e. put the star somewhere
683 out of the galaxy).
684
685 >>> lng = 10.2
686 >>> lat = 9.0
687 >>> ak = findext_marshall(lng, lat, norm='Ak')
688 >>> print("Ak at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, ak))
689 Ak at lng = 10.20, lat = 9.00 is 0.15 magnitude
690
691 2. Find the extinction for a star at galactic lattitude 271.05 and
692 longitude -4.93 and a distance of 144.65 parsecs
693
694 >>> lng = 271.05
695 >>> lat = -4.93
696 >>> dd = 144.65
697 >>> ak = findext_marshall(lng, lat, distance = dd, norm='Ak')
698 >>> print("Ak at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, ak))
699 Ak at lng = 271.05, lat = -4.93 and distance = 144.65 parsecs is 0.02 magnitude
700
701 3. The extinction for galactic longitude 10.2 and latitude 59.0 can not
702 be calculated using the Marshall models, because they are only valid at
703 certain lng and lat (0 < lng < 100 or 260 < lng < 360, -10 < lat < 10)
704
705 >>> lng = 10.2
706 >>> lat = 59.0
707 >>> ak = findext_marshall(lng, lat, norm='Ak')
708 >>> print(ak)
709 None
710
711
712 @param ll: Galactic Longitude (in degrees) should be between 0 and 100 or 260 and 360 degrees
713 @type ll: float
714 @param bb: Galactic Lattitude (in degrees) should be between -10 and 10 degrees
715 @type bb: float
716 @param distance: Distance to the source (in parsecs)
717 @type distance: float
718 @param redlaw: the used reddening law (standard: 'cardelli1989')
719 @type redlaw: str
720 @param Rv: Av/E(B-V) (standard: 3.1)
721 @type Rv: float
722 @return: The extinction in K-band
723 @rtype: float
724 """
725
726
727 data_ma, units_ma, comments_ma = get_marshall_data()
728
729
730 if (ll > 100. and ll < 260.) or ll < 0 or ll > 360:
731 Av = None
732 logger.error("Galactic longitude invalid")
733 return Av
734 elif bb > 10. or bb < -10.:
735 Av = None
736 logger.error("Galactic lattitude invalid")
737 return Av
738
739
740 dist = np.sqrt((data_ma.GLAT - bb)**2. + (data_ma.GLON - ll)**2.)
741 kma = np.where(dist == dist.min())
742
743 if min(dist) > .5:
744 logger.error("Could not find a good model value")
745 return Av
746
747
748 nb = data_ma.nb[kma]
749 rr = [] ; e_rr = [] ; ext = [] ; e_ext = []
750 for i in np.arange(nb)+1.:
751 rr.append(data_ma["r%i"%i][kma])
752 e_rr.append(data_ma["e_r%i"%i][kma])
753 ext.append(data_ma["ext%i"%i][kma])
754 e_ext.append(data_ma["e_ext%i"%i][kma])
755 rr = np.array(rr)[:,0] ; e_rr = np.array(e_rr)[:,0] ; ext = np.array(ext)[:,0] ; e_ext = np.array(e_ext)[:,0]
756
757
758 if distance is None:
759 logger.info("No distance given")
760 dd = max(rr)
761 ak = ext[-1]
762 else:
763 dd = distance/1e3
764
765 if dd < min(rr):
766 logger.info("Distance below distance to first bin")
767 ak = (dd/rr[0])*ext[0]
768 elif dd > max(rr):
769 logger.info("Distance more than distance to last bin")
770 ak = ext[-1]
771 dd = max(rr)
772 else:
773 logger.info("Interpolating linearly")
774 ak = sc.interp(dd, rr, ext)
775
776
777
778 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.K'])
779 return ak/redflux[0]
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804 avdisk = pf.getdata(config.get_datafile('drimmel',"avdisk.fits" ))
805 avloc = pf.getdata(config.get_datafile('drimmel',"avloc.fits" ))
806 avspir = pf.getdata(config.get_datafile('drimmel',"avspir.fits" ))
807 avdloc = pf.getdata(config.get_datafile('drimmel',"avdloc.fits" ))
808 avori = pf.getdata(config.get_datafile('drimmel',"avori.fits" ))
809 coordinates = pf.getdata(config.get_datafile('drimmel',"coordinates.fits" ))
810 avgrid = pf.getdata(config.get_datafile('drimmel',"avgrid.fits" ))
811 avori2 = pf.getdata(config.get_datafile('drimmel',"avori2.fits" ))
812 rf_allsky = pf.getdata(config.get_datafile('drimmel',"rf_allsky.fits" ))
813 glat = rf_allsky.glat
814 glng = rf_allsky.glng
815 ncomp = rf_allsky.ncomp
816 pnum = rf_allsky.pnum
817 rfac = rf_allsky.rfac
818
819 g2e = array([[-0.054882486, -0.993821033, -0.096476249], [0.494116468, -0.110993846, 0.862281440], [-0.867661702, -0.000346354, 0.497154957]])
820
821 -def findext_drimmel(lng, lat, distance=None, rescaling=True,
822 redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):
823 """
824 Procedure to retrieve the absorption in V from three-dimensional
825 grids, based on the Galactic dust distribution of Drimmel & Spergel.
826
827 Example usage:
828
829 1. Find the extinction for a star at galactic longitude 10.2 and
830 latitude 9.0. If the distance is not given, we return the
831 complete extinction along the line of sight (i.e. put the star somewhere
832 out of the galaxy).
833
834 >>> lng = 10.2
835 >>> lat = 9.0
836 >>> av = findext_drimmel(lng, lat)
837 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av))
838 Av at lng = 10.20, lat = 9.00 is 1.24 magnitude
839
840 2. Find the extinction for a star at galactic lattitude 107.05 and
841 longitude -34.93 and a distance of 144.65 parsecs
842
843 >>> lng = 271.05
844 >>> lat = -4.93
845 >>> dd = 144.65
846 >>> ak = findext_marshall(lng, lat, distance = dd)
847 >>> print("Ak at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, ak))
848 Ak at lng = 271.05, lat = -4.93 and distance = 144.65 parsecs is 0.02 magnitude
849
850
851 @param lng: Galactic Longitude (in degrees)
852 @type lng: float
853 @param lat: Galactic Lattitude (in degrees)
854 @type lat: float
855 @param distance: Distance to the source (in parsecs)
856 @type distance: float
857 @param rescaling: Rescaling needed or not?
858 @type rescaling: boolean
859 @return: extinction in V band with/without rescaling
860 @rtype: float
861 """
862
863 deg2rad = pi/180.
864 nsky = 393216
865
866
867 xsun = -8.0
868 zsun = 0.015
869
870
871 if distance is None:
872 d = 1e10
873 else:
874 d = distance/1e3
875
876
877 dfac = ones(nsky)
878 sfac = ones(nsky)
879 lfac = ones(nsky)
880 indx = where(ncomp == 1)
881 dfac[indx] = rfac[indx]
882 indx = where(ncomp == 2)
883 sfac[indx] = rfac[indx]
884 indx = where(ncomp == 3)
885 lfac[indx] = rfac[indx]
886
887
888 num = array(d).size
889 out = zeros(num)
890 avloc = zeros(num)
891 abspir = zeros(num)
892 absdisk = zeros(num)
893
894 l = lng*deg2rad
895 b = lat*deg2rad
896
897
898
899
900 vectoarr = arange(768*512)
901 vectoarr.shape = (512,768)
902 vectoarr = vectoarr.transpose()
903 res = 9
904 pxindex = _ll2pix(lng, lat, res)
905 xout, yout = _pix2xy(pxindex, res, sixpack=True)
906 tblindex = vectoarr[xout, yout]
907 dmax = ones(num)*100.
908 if b != 0.:
909 dmax = .49999/abs(sin(b)) - zsun/sin(b)
910 if cos(l) != 0.:
911 dmax = min([dmax,(14.9999/abs(cos(l)) - xsun/cos(l))])
912 if sin(l) != 0.:
913 dmax = min([dmax, 14.9999/abs(sin(l))])
914
915
916 r = array([d])
917 indx = where(array([d]) > dmax)[0]
918 if len(indx) > 0:
919 r[indx] = dmax
920
921
922 x = array([r*cos(b)*cos(l)])
923 y = array([r*cos(b)*sin(l)])
924 z = array([r*sin(b) + zsun])
925
926
927 i = where(logical_and(abs(x) < 1.,abs(y) < 2.))[0]
928 ni = len(i)
929 j = where(logical_or(abs(x) >= 1., abs(y) >= 2.))[0]
930 nj = len(j)
931
932 if ni > 0:
933
934 dx = 0.02
935 dy = 0.02
936 dz = 0.02
937 nx = 101
938 ny = 201
939 nz = 51
940
941
942 xi = x[i]/dx + float(nx - 1)/2.
943 yj = y[i]/dy + float(ny - 1)/2.
944 zk = z[i]/dz + float(nz - 1)/2.
945
946
947 avloc[i] = _trint(avori2, xi, yj, zk, missing=0.)
948
949
950 k = where(logical_and(abs(x) < 0.75, abs(y) < 0.75))[0]
951 nk = len(k)
952 m = where(logical_or(abs(x) >= 0.75, abs(y) >= 0.75))[0]
953 nm = len(m)
954
955 if nk > 0:
956
957
958 dx = 0.05
959 dy = 0.05
960 dz = 0.02
961 nx = 31
962 ny = 31
963 nz = 51
964
965
966 xi = x[k]/dx + float(nx - 1)/2.
967 yj = y[k]/dy + float(ny - 1)/2.
968 zk = z[k]/dz + float(nz - 1)/2.
969
970
971 absdisk[k] = _trint(avdloc,xi,yj,zk,missing=0.)
972
973
974 x = x + xsun
975
976
977 if nj > 0:
978
979 dmax = 100.
980 if b != 0.:
981 dmax = .49999/abs(sin(b)) - zsun/sin(b)
982 if cos(l) > 0.:
983 dmax = min([dmax, 2.374999/abs(cos(l))])
984 if cos(l) < 0.:
985 dmax = min([dmax, 1.374999/abs(cos(l))])
986 if sin(l) != 0.:
987 dmax = min([dmax, 3.749999/abs(sin(l))])
988
989
990 r1 = array([d])
991 indx = where(array([d]) >= dmax)[0]
992 n = len(indx)
993 if n > 0:
994 r1[indx] = dmax
995
996
997 x1 = array([r1*cos(b)*cos(l) + xsun])
998 y1 = array([r1*cos(b)*sin(l)])
999 z1 = array([r1*sin(b) + zsun])
1000
1001
1002 dx = 0.05
1003 dy = 0.05
1004 dz = 0.02
1005 nx = 76
1006 ny = 151
1007 nz = 51
1008
1009
1010 xi = x1[j]/dx + 2.5*float(nx - 1)
1011 yj = y1[j]/dy + float(ny - 1)/2.
1012 zk = z1[j]/dz + float(nz - 1)/2.
1013
1014
1015 avloc[j] = _trint(avori,xi,yj,zk,missing=0.)
1016
1017
1018 dx = 0.2
1019 dy = 0.2
1020 dz = 0.02
1021 nx = 151
1022 ny = 151
1023 nz = 51
1024
1025
1026 xi = x/dx + float(nx - 1)/2.
1027 yj = y/dy + float(ny - 1)/2.
1028 zk = z/dz + float(nz - 1)/2.
1029
1030
1031 abspir = _trint(avspir,xi,yj,zk,missing=0.)
1032 if nm > 0:
1033 absdisk[m] = _trint(avdisk,xi[m],yj[m],zk[m],missing=0.)
1034
1035
1036 if rescaling:
1037 out = (dfac[tblindex]*absdisk + sfac[tblindex]*abspir + lfac[tblindex]*avloc).flatten()
1038 else:
1039 out = (absdisk + abspir + avloc).flatten()
1040
1041
1042 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'])
1043
1044 return(out/redflux[0])
1045
1046 -def _pix2xy(pixel, resolution, sixpack=0):
1047 """
1048 _pix2xy creates a raster image (sky cube or face) given a pixelindex and a
1049 resolution The data array can be either a vector or two-dimensional array.
1050 In the latter case, the data for each raster image can be stored in either
1051 the columns or rows. The procedure also returns the x and y raster
1052 coordinates of each pixel. Only right oriented, ecliptic coordinate images
1053 are built.
1054 SMOLDERS SEAL OF APPROVAL
1055 """
1056
1057 pixel = array(pixel)
1058 resolution = int(resolution)
1059 newshape = []
1060 for s in pixel.shape:
1061 if s > 1:
1062 newshape.append()
1063 pixel.shape = newshape
1064 pix_sz = (pixel)
1065
1066 if max(pixel) > 6*4**(resolution-1):
1067 raise('Maximum pixel number too large for resolution')
1068
1069
1070 data = [-1]
1071 face = -1
1072 bad_pixval = 0.0
1073
1074 xout, yout = _rastr(pixel,resolution)
1075 return(xout, yout)
1076
1077 -def _rastr(pixel,resolution):
1078 """
1079 SMOLDERS SEAL OF APPROVAL
1080 """
1081 pixel = array(pixel)
1082 n = pixel.size
1083 i0 = 3
1084 j0 = 2
1085 offx = [0,0,1,2,2,1]
1086 offy = [1,0,0,0,1,1]
1087 fij = _pix2fij(pixel,resolution)
1088 cube_side = 2**(resolution-1)
1089 lenc = i0*cube_side
1090 if n > 1:
1091 x_out = offx[fij[0,:]] * cube_side + fij[1,:]
1092 x_out = lenc - (x_out+1)
1093 y_out = offy[fij[0,:]] * cube_side + fij[2,:]
1094 else:
1095 x_out = offx[fij[0]] * cube_side + fij[1]
1096 x_out = lenc - (x_out+1)
1097 y_out = offy[fij[0]] * cube_side + fij[2]
1098 return(x_out, y_out)
1099
1101 """
1102 This function takes an n-element pixel array and generates an n by 3 element
1103 array containing the corresponding face, column, and row number (the latter
1104 two within the face).
1105 SMOLDERS SEAL OF APPROVAL
1106 """
1107
1108 pixel = array(pixel)
1109 n = pixel.size
1110 output = array(zeros((3,n)), int)
1111 res1 = resolution - 1
1112 num_pix_face = 4**res1
1113
1114 face = pixel/num_pix_face
1115 fpix = pixel-num_pix_face*face
1116 output[0,:] = face
1117 pow_2 = 2**arange(16)
1118 ii = array(zeros(n), int)
1119 jj = array(zeros(n), int)
1120
1121 if n > 1:
1122 for i in arange(n):
1123 for bit in arange(res1):
1124 ii[i] = ii[i] | (pow_2[bit]*(1 & fpix))
1125 fpix = fpix >> 1
1126 jj[i] = jj[i] | (pow_2[bit]*(1 & fpix))
1127 fpix = fpix >> 1
1128 else:
1129 for bit in arange(res1):
1130 ii = ii | (pow_2[bit]*(1 & fpix))
1131 fpix = fpix >> 1
1132 jj = jj | (pow_2[bit]*(1 & fpix))
1133 fpix = fpix >> 1
1134 output[1,:] = ii
1135 output[2,:] = jj
1136 return output
1137
1139 gstar = 1.37484847732
1140 g = -0.13161671474
1141 m = 0.004869491981
1142 w1 = -0.159596235474
1143 c00 = 0.141189631152
1144 c10 = 0.0809701286525
1145 c01 = -0.281528535557
1146 c11 = 0.15384112876
1147 c20 = -0.178251207466
1148 c02 = 0.106959469314
1149 d0 = 0.0759196200467
1150 d1 = -0.0217762490699
1151 r0 = 0.577350269
1152 aa = alpha**2.
1153 bb = beta**2
1154 a4 = aa**2
1155 b4 = bb**2
1156 onmaa = 1.-aa
1157 onmbb = 1.-bb
1158 x = alpha*(gstar + aa*(1.-gstar) + onmaa*(bb*(g+(m-g)*aa + onmbb*(c00+c10*aa+c01*bb+c11*aa*bb+c20*a4+c02*b4)) + aa*(w1-onmaa*(d0+d1*aa))))
1159 y = beta *(gstar + bb*(1.-gstar) + onmbb*(aa*(g+(m-g)*bb + onmaa*(c00+c10*bb+c01*aa+c11*bb*aa+c20*b4+c02*a4)) + bb*(w1-onmbb*(d0+d1*bb))))
1160 return x, y
1161
1163 """
1164 Convert longitude, latitude to unit vectors
1165 SMOLDERS SEAL OF APPROVAL
1166
1167 @param lng : galactic longitude
1168 @type lng : float
1169 @param lat : galactic lattitude
1170 @type lat : float
1171 @return : unitvector
1172 @rtype : ndarray
1173 """
1174 vector = zeros(3)
1175 d2r = pi/180
1176 lng = lng * d2r
1177 lat = lat * d2r
1178 vector[0] = cos(lat) * cos(lng)
1179 vector[1] = cos(lat) * sin(lng)
1180 vector[2] = sin(lat)
1181 return vector
1182
1184 """
1185 Unitvector for galactic coordinates to unitvector for ecliptic coordinates
1186 SMOLDERS SEAL OF APPROVAL
1187 """
1188 out_uvec = dot(in_uvec, g2e)
1189 return out_uvec
1190
1192 """
1193 converts unitvector into nface number (0-5) and X,Y in range 0-1
1194 SMOLDERS SEAL OF APPROVAL
1195 """
1196 vector = array(vector)
1197 vs = vector.shape
1198 if len(vs) > 1:
1199 vec0 = array(vector[:,0])
1200 vec1 = array(vector[:,1])
1201 vec2 = array(vector[:,2])
1202 else:
1203 vec0 = array([vector[0]])
1204 vec1 = array([vector[1]])
1205 vec2 = array([vector[2]])
1206 abs_yx = abs(vec1/vec0)
1207 abs_zx = abs(vec2/vec0)
1208 abs_zy = abs(vec2/vec1)
1209
1210 nface = zeros(len(vec0))
1211 for i in arange(len(vec0)):
1212 nface[i] = (0 * (abs_zx[i] >= 1 and abs_zy[i] >= 1 and vec2[i] >= 0) +
1213 5 * (abs_zx[i] >= 1 and abs_zy[i] >= 1 and vec2[i] < 0) +
1214 1 * (abs_zx[i] < 1 and abs_yx[i] < 1 and vec0[i] >= 0) +
1215 3 * (abs_zx[i] < 1 and abs_yx[i] < 1 and vec0[i] < 0) +
1216 2 * (abs_zy[i] < 1 and abs_yx[i] >= 1 and vec1[i] >= 0) +
1217 4 * (abs_zy[i] < 1 and abs_yx[i] >= 1 and vec1[i] < 0))
1218
1219 nface_0 = (nface == 0)*1.
1220 nface_1 = (nface == 1)*1.
1221 nface_2 = (nface == 2)*1.
1222 nface_3 = (nface == 3)*1.
1223 nface_4 = (nface == 4)*1.
1224 nface_5 = (nface == 5)*1.
1225
1226 eta = (vec2/vec0)*(nface_1 - nface_3) + (vec2/vec1)*(nface_2 - nface_4) - (vec0/vec2)*(nface_0 + nface_5)
1227 xi = (vec1/vec0)*(nface_1 + nface_3) - (vec0/vec1)*(nface_2 + nface_4) + (vec1/vec2) * (nface_0 - nface_5)
1228 x, y = _incube(xi,eta)
1229 x = (x+1.)/2.
1230 y = (y+1.)/2.
1231 return x,y,array(nface, dtype=int)
1232
1234 """
1235 This function takes an n by 3 element vector containing the face, column, and
1236 row number (the latter two within the face) of a pixel and converts it into an
1237 n-element pixel array for a given resolution.
1238 """
1239
1240 fij = array(fij)
1241 n = fij.size/3
1242
1243 pixel = array(zeros(n), dtype=int)
1244 pixel_1 = array(zeros(n), dtype=int)
1245
1246 if n > 1:
1247 ff = array(fij[0,:], dtype=int)
1248 ii = array(fij[1,:], dtype=int)
1249 jj = array(fij[2,:], dtype=int)
1250 else:
1251 ff = int(fij[0])
1252 ii = int(fij[1])
1253 jj = int(fij[2])
1254
1255 num_pix_face = 4**(res-1)
1256 pow_2 = 2**arange(16)
1257
1258
1259 if n > 1:
1260 for i in arange(n):
1261 for bit in arange(res-1):
1262 pixel_1[i] = pixel_1[i] | ((pow_2[bit] & ii[i]) << bit)
1263 pixel_1[i] = pixel_1[i] | ((pow_2[bit] & jj[i]) << (bit+1))
1264 else:
1265 for bit in arange(res-1):
1266 pixel_1 = pixel_1 | ((pow_2[bit] & ii) << bit)
1267 pixel_1 = pixel_1 | ((pow_2[bit] & jj) << (bit+1))
1268
1269 pixel = ff*num_pix_face + pixel_1
1270 return pixel
1271
1273 """
1274 Routine returns pixel number given unit vector pointing to center of pixel
1275 resolution of the cube.
1276 SMOLDERS SEAL OF APPROVAL
1277 """
1278 two = 2
1279 vector = array(vector)
1280 n_vec = int(vector.size/3)
1281 res1 = resolution - 1
1282 num_pix_side = int(two**res1)
1283 x, y, face = _axisxy(vector)
1284 ia = array(x*num_pix_side, dtype=int)
1285 ib = array([2.**res1 - 1]*n_vec, dtype=int)
1286 ja = array(y*num_pix_side, dtype=int)
1287 jb = array([2.**res1 - 1]*n_vec, dtype=int)
1288 i = zeros(n_vec)
1289 j = zeros(n_vec)
1290 for k in arange(n_vec):
1291 i[k] = min([ia[k],ib[k]])
1292 j[k] = min([ja[k],jb[k]])
1293 pixel = _fij2pix(array([face,i,j]),resolution)
1294 return pixel
1295
1297 """
1298 _ll2pix is a python function to convert galactic coordinates to DIRBE pixels
1299 (coorconv([lng,lat], infmt='L', outfmt='P', inco='G', outco='R9'))
1300
1301 Input
1302 @param lng : galactic longitude
1303 @type lng : float
1304 @param lat : galactic lattitude
1305 @type lat : float
1306 @return: output coordinate array
1307 @rtype: ndarray
1308
1309 (Note: The default coordinate system for pixels is ecliptic.)
1310 """
1311
1312 uvec = _ll2uv(lng, lat)
1313 uvec = _galvec2eclvec(uvec)
1314 output = _uv2pix(uvec,res)
1315 return output
1316
1317 -def _trint(mm,x,y,z,missing=None):
1318 """
1319 Pixel-based trilinear interpolation, where mm is a 3d data cube.
1320 """
1321 xmi = array(floor(x),int)
1322 xd = abs(xmi-x)
1323 xma = array(ceil(x),int)
1324 ymi = array(floor(y), int)
1325 yd = abs(ymi-y)
1326 yma = array(ceil(y),int)
1327 zmi = array(floor(z), int)
1328 zd = abs(zmi-z)
1329 zma = array(ceil(z), int)
1330 i1 = mm[zmi,ymi,xmi]*(1.-zd) + mm[zma,ymi,xmi]*(zd)
1331 i2 = mm[zmi,yma,xmi]*(1.-zd) + mm[zma,yma,xmi]*(zd)
1332 j1 = mm[zmi,ymi,xma]*(1.-zd) + mm[zma,ymi,xma]*(zd)
1333 j2 = mm[zmi,yma,xma]*(1.-zd) + mm[zma,yma,xma]*(zd)
1334 w1 = i1*(1-yd) + i2*yd
1335 w2 = j1*(1-yd) + j2*yd
1336 iv = w1*(1.-xd) + w2*xd
1337 return iv
1338
1351
1352 dustname = config.get_datafile('schlegel',"SFD_dust_4096_sgp.fits")
1353 maskname = config.get_datafile('schlegel',"SFD_mask_4096_sgp.fits")
1354 data = pf.getdata(dustname)
1355 mask = pf.getdata(maskname)
1356 return data, mask
1357
1359
1360 dustname = config.get_datafile('schlegel',"SFD_dust_4096_ngp.fits")
1361 maskname = config.get_datafile('schlegel',"SFD_mask_4096_ngp.fits")
1362 data = pf.getdata(dustname)
1363 mask = pf.getdata(maskname)
1364 return data, mask
1365
1367 """
1368 Converts coordinates in lattitude and longitude to coordinates
1369 in x and y pixel coordinates.
1370
1371 The x and y coordinates you find with these formulas are not
1372 the coordinates you can read in ds9, but 1 smaller. Hence
1373 (x+1, y+1) are the coordinates you find in DS9.
1374
1375 Input
1376 @param ll : galactic longitude
1377 @type ll : float
1378 @param bb : galactic lattitude
1379 @type bb : float
1380 @return: output coordinate array
1381 @rtype: ndarray
1382 """
1383 deg2rad = pi/180.
1384
1385 if bb <= 0 : hs = -1.
1386 elif bb > 0: hs = +1.
1387
1388 yy = 2048 * sqrt(1. - hs * sin(bb*deg2rad)) * cos(ll*deg2rad) + 2047.5
1389 xx = -2048 * hs * sqrt(1 - hs * sin(bb*deg2rad)) * sin(ll*deg2rad) + 2047.5
1390
1391 return xx, yy
1392
1393 -def findext_schlegel(ll, bb, distance = None, redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):
1394 """
1395 Get the "Schlegel" extinction at certain longitude and latitude
1396
1397 This function returns the E(B-V) maps of Schlegel et al. (1998),
1398 depending on wether the distance is given or not, the E(B-V) value
1399 is corrected. If distance is set we use the distance-corrected
1400 values:
1401
1402 E(B-V) = E(B-V)_maps * (1 - exp(-10 * r * sin(|b|)))
1403 where E(B-V) is the value to be used and E(B-V)_maps the value
1404 as found with the Schlegel dust maps
1405
1406 Then we convert the E(B-V) to Av. Standard we use Av = E(B-V)*Rv with Rv=3.1, but the value of Rv can be given as a keyword.
1407
1408 ! WARNING: the schlegel maps are not usefull when |b| < 5 degrees !
1409 """
1410 deg2rad = pi/180.
1411 dd = distance
1412 if distance is not None:
1413 dd = distance/1.e3
1414
1415
1416
1417 xx, yy = _lb2xy_schlegel(ll,bb)
1418
1419
1420 if bb <= 0:
1421 data, mask = get_schlegel_data_south()
1422 elif bb > 0:
1423 data, mask = get_schlegel_data_north()
1424
1425 if abs(bb) < 10.:
1426 logger.warning("Schlegel is not good for lattitudes > 10 degrees")
1427
1428
1429 xl = floor(xx)
1430 yl = floor(yy)
1431 xh = xl + 1.
1432 yh = yl + 1.
1433
1434
1435 w1 = (xl-xx)**2 + (yl-yy)**2
1436 w2 = (xl-xx)**2 + (yh-yy)**2
1437 w3 = (xh-xx)**2 + (yl-yy)**2
1438 w4 = (xh-xx)**2 + (yh-yy)**2
1439
1440
1441 v1 = data[xl, yl]
1442 v2 = data[xl, yh]
1443 v3 = data[xh, yl]
1444 v4 = data[xh, yh]
1445 f1 = mask[xl, yl]
1446 f2 = mask[xl, yh]
1447 f3 = mask[xh, yl]
1448 f4 = mask[xh, yh]
1449
1450 ebv = (w1*v1 + w2*v2 + w3*v3 + w4*v4) / (w1 + w2 + w3 + w4)
1451
1452
1453 logger.info("flag of pixel 1 is: %i" %f1)
1454 logger.info("flag of pixel 2 is: %i" %f2)
1455 logger.info("flag of pixel 3 is: %i" %f3)
1456 logger.info("flag of pixel 4 is: %i" %f4)
1457
1458 if dd is not None:
1459 ebv = ebv * (1. - exp(-10. * dd * sin(abs(bb*deg2rad))))
1460
1461
1462 av = ebv*Rv
1463
1464
1465 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'])
1466
1467 return av/redflux[0]
1468
1469
1470
1471 if __name__ == "__main__":
1472 print findext_marshall(10.2,9.)
1473 print findext_marshall(10.2,9.,norm='Ak')
1474