ComboCode.cc.ivs.sed.extinctionmodels

1 # -*- coding: utf-8 -*- 2 """ 3 Extinction models of Arenou, Drimmel and Marschall 4 5 Uniformly rewritten by K. Smolders 6 """ 7 8 from cc.ivs.catalogs import vizier 9 from cc.ivs.sed.reddening import get_law 10 from cc.ivs.aux.decorators import memoized 11 from cc.ivs.aux import loggers 12 import cc.path 13 #from ivs import config 14 15 import os 16 import numpy as np 17 from numpy import (abs, arange, array, ceil, cos, dot, floor, int, logical_and, 18 logical_or, max, min, ones, pi, sin, sqrt, where, zeros, exp) 19 import scipy as sc 20 from astropy.io import fits as pf 21 import logging 22 23 logger = logging.getLogger("SED.EXT") 24 logger.addHandler(loggers.NullHandler)

25 26 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 27 #{ Top level wrapper 28 # # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 29 -def findext(lng, lat, model='drimmel', distance=None, **kwargs):

30 """ 31 Get the "model" extinction at a certain longitude and latitude. 32 33 Find the predicted V-band extinction (Av) based on three 34 dimensional models of the galactic interstellar extinction. 35 The user can choose between different models by setting the 36 model keyword: 37 38 1) "arenou": model from Arenou et al. (1992). 39 2) "schlegel": model from Schlegel et al. (1998) 40 3) "drimmel": model from Drimmel et al. (2003) 41 4) "marshall": model from Marshall et al. (2006) 42 43 example useage: 44 45 1. Find the total galactic extinction for a star at galactic lattitude 10.2 46 and longitude 59.0 along the line of sight, as given by the model of 47 Arenou et al. (1992) 48 49 >>> lng = 10.2 50 >>> lat = 59.0 51 >>> av = findext(lng, lat, model='arenou') 52 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av)) 53 Av at lng = 10.20, lat = 59.00 is 0.05 magnitude 54 55 2. Find the extinction for a star at galactic lattitude 107.05 and 56 longitude -34.93 and a distance of 144.65 parsecs, as given by the 57 model of Arenou et al. (1992) 58 59 >>> lng = 107.05 60 >>> lat = -34.93 61 >>> dd = 144.65 62 >>> av = findext(lng, lat, distance = dd, model='arenou') 63 >>> print("Av at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, av)) 64 Av at lng = 107.05, lat = -34.93 and distance = 144.65 parsecs is 0.15 magnitude 65 66 3. Find the Marschall extinction for a star at galactic longitude 10.2 and 67 latitude 9.0. If the distance is not given, we return the 68 complete extinction along the line of sight (i.e. put the star somewhere 69 out of the galaxy). 70 71 >>> lng = 10.2 72 >>> lat = 9.0 73 >>> av = findext(lng, lat, model='marshall') 74 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av)) 75 Av at lng = 10.20, lat = 9.00 is 10.67 magnitude 76 77 4. Find the Marshall extinction for a star at galactic lattitude 271.05 and 78 longitude -4.93 and a distance of 144.65 parsecs, but convert to Av using 79 Rv = 2.5 instead of Rv = 3.1 80 81 >>> lng = 271.05 82 >>> lat = -4.93 83 >>> dd = 144.65 84 >>> av = findext(lng, lat, distance = dd, model='marshall', Rv=2.5) 85 >>> print("Av at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, av)) 86 Av at lng = 271.05, lat = -4.93 and distance = 144.65 parsecs is 13.95 magnitude 87 88 5. Find the extinction for a star at galactic longitude 10.2 and 89 latitude 9.0, using the schlegel model, using Rv=2.9 instead of 90 Rv=3.1 91 92 >>> lng = 58.2 93 >>> lat = 24.0 94 >>> distance = 848. 95 >>> av = findext(lng, lat, distance=distance) 96 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av)) 97 Av at lng = 58.20, lat = 24.00 is 0.12 magnitude 98 99 REMARKS: 100 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) 101 b) Schlegel is very dubious for latitudes between -5 and 5 degrees 102 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') 103 d) Marschall is only available for certain longitudes and latitudes: 104 0 < lng < 100 or 260 < lng < 360 and -10 < lat < 10 105 106 @param lng: Galactic Longitude (in degrees) 107 @type lng: float 108 @param lat: Galactic Lattitude (in degrees) 109 @type lat: float 110 @param model: the name of the extinction model: ("arenou", "schlegel", "drimmel" or "marshall"; if none given, the program uses "drimmel") 111 @type model: str 112 @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 113 @type distance: float 114 @return: The extinction in Johnson V-band 115 @rtype: float 116 """ 117 118 if model.lower() == 'drimmel': 119 av = findext_drimmel(lng, lat, distance=distance, **kwargs) 120 elif model.lower() == 'marshall' or model.lower() == 'marschall': 121 av = findext_marshall(lng, lat, distance=distance, **kwargs) 122 elif model.lower() == 'arenou': 123 av = findext_arenou(lng, lat, distance=distance, **kwargs) 124 elif model.lower() == 'schlegel': 125 av = findext_schlegel(lng, lat, distance=distance, **kwargs) 126 return(av)

127

128 #} 129 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 130 #{ Arenou 3D extinction model 131 # (based on Arenou et al, "A tridimensional model of the 132 # galactic interstellar extinction" published in Astronomy and 133 # Astrophysics (ISSN 0004-6361), vol. 258, no. 1, p. 104-111.) 134 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 135 136 -def findext_arenou(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):

137 """ 138 Find the predicted V-band extinction (Av) according to the 139 3D model for galactic extinction of Arenou et al, "Modelling 140 the Galactic interstellar extinction distribution in three 141 dimensions", Arenou et al, "A tridimensional model of the 142 galactic interstellar extinction" published in Astronomy and 143 Astrophysics (ISSN 0004-6361), vol. 258, no. 1, p. 104-111, 1992 144 145 Example usage: 146 147 1. Find the extinction for a star at galactic lattitude 10.2 and 148 longitude 59.0. If the distance is not given, we use the 149 maximal distance r0 for that line of sight, as given in the 150 Appendix of Arenou et al. (1992) 151 152 >>> lng = 10.2 153 >>> lat = 59.0 154 >>> av = findext_arenou(lng, lat) 155 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av)) 156 Av at lng = 10.20, lat = 59.00 is 0.05 magnitude 157 158 2. Find the extinction for a star at galactic lattitude 107.05 and 159 longitude -34.93 and a distance of 144.65 parsecs 160 161 >>> lng = 107.05 162 >>> lat = -34.93 163 >>> dd = 144.65 164 >>> av = findext_arenou(lng, lat, distance = dd) 165 >>> print("Av at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, av)) 166 Av at lng = 107.05, lat = -34.93 and distance = 144.65 parsecs is 0.15 magnitude 167 168 @param ll: Galactic Longitude (in degrees) 169 @type ll: float 170 @param bb: Galactic Lattitude (in degrees) 171 @type bb: float 172 @param distance: Distance to the source (in parsecs) 173 @type distance: float 174 @return: The extinction in Johnson V-band 175 @rtype: float 176 """ 177 # make sure that the values for b and l are within the correct range 178 if bb < -90. or bb > 90: 179 logger.error("galactic lattitude outside [-90,90] degrees") 180 elif ll < 0. or ll > 360: 181 logger.error("galactic longitude outside [0,360] degrees") 182 elif distance < 0 and distance is not None: 183 logger.error("distance is negative") 184 185 # find the Arenou paramaters in the Appendix of Arenou et al. (1992) 186 alpha, beta, gamma, rr0, saa = _getarenouparams(ll, bb) 187 logger.info("Arenou params: alpha = %.2f, beta = %.2f, gamma = %.2f, r0 = %.2f and saa = %.2f" %(alpha, beta, gamma, rr0, saa)) 188 189 # compute the visual extinction from the Arenou paramaters using Equation 5 190 # and 5bis 191 if distance is None: 192 av = alpha*rr0 + beta*rr0**2. 193 else: 194 distance = distance/1e3 # to kparsec 195 if distance <= rr0: 196 av = alpha*distance + beta*distance**2. 197 else: 198 av = alpha*rr0 + beta*rr0**2. + (distance-rr0)*gamma 199 200 #-- Marshall is standard in Ak, but you can change this: 201 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'],**kwargs) 202 203 return av/redflux[0]

204

205 -def _getarenouparams(ll,bb):

206 """ 207 Input: galactic coordinates 208 Output: Arenou 1992 alpha, beta, gamma, rr0, saa 209 210 @param ll: Galactic Longitude (in degrees) 211 @type ll: float 212 @param bb: Galactic Lattitude (in degrees) 213 @type bb: float 214 """ 215 if -90 < bb < -60: 216 gamma = 0 217 if 0 <= ll < 29: 218 alpha = 2.22538 ; beta = -6.00212 ; rr0 = 0.052 ; saa = 13. 219 elif 29 <= ll < 57: 220 alpha = 3.35436 ; beta = -14.74567 ; rr0 = 0.035 ; saa = 40 221 elif 57 <= ll < 85: 222 alpha = 2.77637 ; beta = -9.62706 ; rr0 = 0.042 ; saa = 15 223 elif 85 <= ll < 110: 224 alpha = 4.44190 ; beta = -19.92097 ; rr0 = 0.025 ; saa = 36 225 elif 110 <= ll < 150: 226 alpha = 4.46685 ; beta = -26.07305 ; rr0 = 0.026 ; saa = 28 227 elif 150 <= ll < 180: 228 alpha = 7.63699 ; beta = -46.10856 ; rr0 = 0.014 ; saa = 38 229 elif 180 <= ll < 210: 230 alpha = 2.43412 ; beta = -8.69913 ; rr0 = 0.050 ; saa = 36 231 elif 210 <= ll < 240: 232 alpha = 3.34481 ; beta = -13.93228 ; rr0 = 0.035 ; saa = 38 233 elif 240 <= ll < 270: 234 alpha = 1.40733 ; beta = -13.43418 ; rr0 = 0.091 ; saa = 30 235 elif 270 <= ll < 300: 236 alpha = 1.64466 ; beta = -3.97380 ; rr0 = 0.074 ; saa = 28 237 elif 300 <= ll < 330: 238 alpha = 2.12696 ; beta = -6.05682 ; rr0 = 0.056 ; saa = 14 239 elif 330 <= ll < 360: 240 alpha = 2.34636 ; beta = -8.17784 ; rr0 = 0.052 ; saa = 16 241 242 if -60 < bb < -45: 243 gamma = 0 244 if 0 <= ll < 30: 245 alpha = 2.77060 ; beta = -9.52310 ; rr0 = 0.145 ; saa = 16 246 elif 30 <= ll < 60: 247 alpha = 1.96533 ; beta = -9.52310 ; rr0 = 0.174 ; saa = 06 248 elif 60 <= ll < 110: 249 alpha = 1.93622 ; beta = -13.31757 ; rr0 = 0.073 ; saa = 26 250 elif 110 <= ll < 180: 251 alpha = 1.05414 ; beta = -2.236540 ; rr0 = 0.223 ; saa = 74 252 elif 180 <= ll < 210: 253 alpha = 1.39990 ; beta = -1.35325 ; rr0 = 0.252 ; saa = 10 254 elif 210 <= ll < 240: 255 alpha = 2,73481 ; beta = -11.70266 ; rr0 = 0.117 ; saa = 8 256 elif 240 <= ll < 270: 257 alpha = 2.99784 ; beta = -11.64272 ; rr0 = 0.129 ; saa = 3 258 elif 270 <= ll < 300: 259 alpha = 3.23802 ; beta = -11.63810 ; rr0 = 0.139 ; saa = 7 260 elif 300 <= ll < 330: 261 alpha = 1.72679 ; beta = -6.05085 ; rr0 = 0.143 ; saa = 7 262 elif 330 <= ll < 360: 263 alpha = 1.88890 ; beta = -5.51861 ; rr0 = 0.171 ; saa = 14 264 265 if -45 < bb < -30: 266 gamma = 0 267 if 0 <= ll < 30: 268 alpha = 1.98973 ; beta = -4.86206 ; rr0 = 0.205 ; saa = 6 269 elif 30 <= ll < 60: 270 alpha = 1.49901 ; beta = -3.75837 ; rr0 = 0.199 ; saa = 16 271 elif 60 <= ll < 90: 272 alpha = 0.90091 ; beta = -1.30459 ; rr0 = 0.329 ; saa = 73 273 elif 90 <= ll < 120: 274 alpha = 1.94200 ; beta = -6.26833 ; rr0 = 0.155 ; saa = 18 275 elif 120 <= ll < 160: 276 alpha = -0.37804 ; beta = 10.77372 ; rr0 = 0.210 ; saa = 100 277 elif 160 <= ll < 200: 278 alpha = -0.15710 ; beta = 5.03190 ; rr0 = 0.294 ; saa = 24 279 elif 200 <= ll < 235: 280 alpha = 3.20162 ; beta = -10.59297 ; rr0 = 0.151 ; saa = 9 281 elif 235 <= ll < 265: 282 alpha = 1.95079 ; beta = 4.73280 ; rr0 = 0.206 ; saa = 21 283 elif 265 <= ll < 300: 284 alpha = 1.91200 ; beta = 4.97640 ; rr0 = 0.192 ; saa = 17 285 elif 300 <= ll < 330: 286 alpha = 2.50487 ; beta = -9.63106 ; rr0 = 0.145 ; saa = 28 287 elif 330 <= ll < 360: 288 alpha = 2.44394 ; beta = -9.17612 ; rr0 = 0.133 ; saa = 7 289 290 if -30 < bb < -15: 291 gamma = 0 292 if 0 <= ll < 20: 293 alpha = 2.82440 ; beta = -4.78217 ; rr0 = 0.295 ; saa = 32 294 elif 20 <= ll < 40: 295 alpha = 3.84362 ; beta = -8.04690 ; rr0 = 0.239 ; saa = 46 296 elif 40 <= ll < 80: 297 alpha = 0.60365 ; beta = 0.07893 ; rr0 = 0.523 ; saa = 22 298 elif 80 <= ll < 100: 299 alpha = 0.58307 ; beta = -0.21053 ; rr0 = 0.523 ; saa = 53 300 elif 100 <= ll < 120: 301 alpha = 2.03861 ; beta = -4.40843 ; rr0 = 0.231 ; saa = 60 302 elif 120 <= ll < 140: 303 alpha = 1.14271 ; beta = -1.35635 ; rr0 = 0.421 ; saa = 34 304 elif 140 <= ll < 160: 305 alpha = 0.79908 ; beta = 1.48074 ; rr0 = 0.513 ; saa = 61 306 elif 160 <= ll < 180: 307 alpha = 0.94260 ; beta = 8.16346 ; rr0 = 0.441 ; saa = 42 308 elif 180 <= ll < 200: 309 alpha = 1.66398 ; beta = 0.26775 ; rr0 = 0.523 ; saa = 42 310 elif 200 <= ll < 220: 311 alpha = 1.08760 ; beta = -1.02443 ; rr0 = 0.523 ; saa = 45 312 elif 220 <= ll < 240: 313 alpha = 1.20087 ; beta = -2.45407 ; rr0 = 0.245 ; saa = 6 314 elif 240 <= ll < 260: 315 alpha = 1.13147 ; beta = -1.87916 ; rr0 = 0.301 ; saa = 16 316 elif 260 <= ll < 280: 317 alpha = 0.97804 ; beta = -2.92838 ; rr0 = 0.338 ; saa = 21 318 elif 290 <= ll < 300: 319 alpha = 1.40086 ; beta = -1.12403 ; rr0 = 0.523 ; saa = 19 320 elif 300 <= ll < 320: 321 alpha = 2.06355 ; beta = -3.68278 ; rr0 = 0.280 ; saa = 42 322 elif 320 <= ll < 340: 323 alpha = 1.59260 ; beta = -2.18754 ; rr0 = 0.364 ; saa = 15 324 elif 310 <= ll < 360: 325 alpha = 1.45589 ; beta = -1.90598 ; rr0 = 0.382 ; saa = 21 326 327 if -15 < bb < -5: 328 gamma = 0 329 if 0 <= ll < 10: 330 alpha = 2.56438 ; beta = -2.31586 ; rr0 = 0.554 ; saa = 37 331 elif 10 <= ll < 20: 332 alpha = 3.24095 ; beta = -2.78217 ; rr0 = 0.582 ; saa = 38 333 elif 20 <= ll < 30: 334 alpha = 2.95627 ; beta = -2.57422 ; rr0 = 0.574 ; saa = 41 ; gamma = 0.08336 335 elif 30 <= ll < 40: 336 alpha = 1.85158 ; beta = -0.67716 ; rr0 = 1.152 ; saa = 4 337 elif 40 <= ll < 50: 338 alpha = 1.60770 ; beta = 0.35279 ; rr0 = 0.661 ; saa = 24 339 elif 50 <= ll < 60: 340 alpha = 0.69920 ; beta = -0.09146 ; rr0 = 0.952 ; saa = 2 ; gamma = 0.12839 341 elif 60 <= ll < 80: 342 alpha = 1.36189 ; beta = -1.05290 ; rr0 = 0.647 ; saa = 45 ; gamma = 0.16258 343 elif 80 <= ll < 90: 344 alpha = 0.33179 ; beta = 0.37629 ; rr0 = 1.152 ; saa = 62 345 elif 90 <= ll < 100: 346 alpha = 1.70303 ; beta = -0.75246 ; rr0 = 1.132 ; saa = 31 347 elif 100 <= ll < 110: 348 alpha = 1.97414 ; beta = -1.59784 ; rr0 = 0.618 ; saa = 35 ; gamma = 0.12847 349 elif 110 <= ll < 120: 350 alpha = 1.07407 ; beta = -0.40066 ; rr0 = 1.152 ; saa = 14 ; gamma = 0.17698 351 elif 120 <= ll < 130: 352 alpha = 1.69495 ; beta = -1.00071 ; rr0 = 0.847 ; saa = 28 ; gamma = 0.08567 353 elif 130 <= ll < 140: 354 alpha = 1.51449 ; beta = -0.08441 ; rr0 = 0.897 ; saa = 12 355 elif 140 <= ll < 150: 356 alpha = 1.87894 ; beta = -0.73314 ; rr0 = 1.152 ; saa = 23 357 elif 150 <= ll < 160: 358 alpha = 1.43670 ; beta = 0.67706 ; rr0 = 0.778 ; saa = 3 ; gamma = 0.42624 359 elif 160 <= ll < 180: 360 alpha = 6.84802 ; beta = -5.06864 ; rr0 = 0.676 ; saa = 50 361 elif 180 <= ll < 190: 362 alpha = 4.16321 ; beta = -5.80016 ; rr0 = 0.359 ; saa = 51 ; gamma = 0.60387 363 elif 190 <= ll < 200: 364 alpha = 0.78135 ; beta = -0.27826 ; rr0 = 1.152 ; saa = 4 365 elif 200 <= ll < 210: 366 alpha = 0.85535 ; beta = 0.20848 ; rr0 = 1.152 ; saa = 17 367 elif 210 <= ll < 220: 368 alpha = 0.52521 ; beta = 0.65726 ; rr0 = 1.152 ; saa = 7 369 elif 220 <= ll < 230: 370 alpha = 0.88376 ; beta = -0.44519 ; rr0 = 0.993 ; saa = 40 ; gamma = 0.06013 371 elif 230 <= ll < 240: 372 alpha = 0.42228 ; beta = 0.26304 ; rr0 = 0.803 ; saa = 26 373 elif 240 <= ll < 250: 374 alpha = 0.71318 ; beta = -0.67229 ; rr0 = 0.530 ; saa = 55 ; gamma = 0.20994 375 elif 250 <= ll < 260: 376 alpha = 0.99606 ; beta = -0.70103 ; rr0 = 0.710 ; saa = 48 ; gamma = 0.01323 377 elif 260 <= ll < 270: 378 alpha = 0.91519 ; beta = -0.39690 ; rr0 = 1.152 ; saa = 48 ; gamma = 0.01961 379 elif 270 <= ll < 280: 380 alpha = 0.85791 ; beta = -0.29115 ; rr0 = 1.152 ; saa = 19 381 elif 280 <= ll < 290: 382 alpha = 1.44226 ; beta = -1.09775 ; rr0 = 0.657 ; saa = 39 383 elif 290 <= ll < 300: 384 alpha = 2.55486 ; beta = -1.68293 ; rr0 = 0.759 ; saa = 31 385 elif 300 <= ll < 310: 386 alpha = 3.18047 ; beta = -2.69796 ; rr0 = 0.589 ; saa = 40 387 elif 210 <= ll < 320: 388 alpha = 2.11235 ; beta = -1.77506 ; rr0 = 0.595 ; saa = 29 389 elif 320 <= ll < 330: 390 alpha = 1.75884 ; beta = -1.38574 ; rr0 = 0.635 ; saa = 25 391 elif 330 <= ll < 340: 392 alpha = 1.97257 ; beta = -1.55545 ; rr0 = 0.634 ; saa = 34 ; gamma = 0.00043 393 elif 340 <= ll < 350: 394 alpha = 1.41497 ; beta = -1.05722 ; rr0 = 0.669 ; saa = 46 ; gamma = 0.03264 395 elif 350 <= ll < 360: 396 alpha = 1.17795 ; beta = -0.95012 ; rr0 = 0.620 ; saa = 46 ; gamma = 0.03339 397 398 if -5 < bb < 5: 399 gamma = 0 400 if 0 <= ll < 10: 401 alpha = 2.62556 ; beta = -1.11097 ; rr0 = 1.182 ; saa = 57 ; gamma = 0.00340 402 elif 10 <= ll < 20: 403 alpha = 3.14461 ; beta = -1.01140 ; rr0 = 1.555 ; saa = 42 404 elif 20 <= ll < 30: 405 alpha = 4.26624 ; beta = -1.61242 ; rr0 = 1.323 ; saa = 34 406 elif 30 <= ll < 40: 407 alpha = 2.54447 ; beta = -0.12771 ; rr0 = 1.300 ; saa = 30 408 elif 40 <= ll < 50: 409 alpha = 2.27030 ; beta = -0.68720 ; rr0 = 1.652 ; saa = 52 ; gamma = 0.04928 410 elif 50 <= ll < 60: 411 alpha = 1.34359 ; beta = -0.05416 ; rr0 = 2.000 ; saa = 32 412 elif 60 <= ll < 70: 413 alpha = 1.76327 ; beta = -0.26407 ; rr0 = 2.000 ; saa = 37 414 elif 70 <= ll < 80: 415 alpha = 2.20666 ; beta = -0.41651 ; rr0 = 2.000 ; saa = 36 416 elif 80 <= ll < 90: 417 alpha = 1.50130 ; beta = -0.09855 ; rr0 = 1.475 ; saa = 45 418 elif 90 <= ll < 100: 419 alpha = 2.43965 ; beta = -0.82128 ; rr0 = 1.485 ; saa = 36 ; gamma = 0.01959 420 elif 100 <= ll < 110: 421 alpha = 3.35775 ; beta = -1.16400 ; rr0 = 0.841 ; saa = 35 ; gamma = 0.00298 422 elif 110 <= ll < 120: 423 alpha = 2.60621 ; beta = -0.68687 ; rr0 = 1.897 ; saa = 36 424 elif 120 <= ll < 130: 425 alpha = 2.90112 ; beta = -0.97988 ; rr0 = 1.480 ; saa = 32 426 elif 130 <= ll < 140: 427 alpha = 2.55377 ; beta = -0.71214 ; rr0 = 1.793 ; saa = 38 428 elif 140 <= ll < 150: 429 alpha = 3.12598 ; beta = -1.21437 ; rr0 = 1.287 ; saa = 23 ; gamma = 0.15298 430 elif 150 <= ll < 160: 431 alpha = 3.66930 ; beta = -2.29731 ; rr0 = 0.799 ; saa = 40 ; gamma = 0.33473 432 elif 160 <= ll < 170: 433 alpha = 2.15465 ; beta = -0.70690 ; rr0 = 1.524 ; saa = 37 ; gamma = 0.14017 434 elif 170 <= ll < 180: 435 alpha = 1.82465 ; beta = -0.60223 ; rr0 = 1.515 ; saa = 29 ; gamma = 0.20730 436 elif 180 <= ll < 190: 437 alpha = 1.76269 ; beta = -0.35945 ; rr0 = 2.000 ; saa = 28 ; gamma = 0.08052 438 elif 190 <= ll < 200: 439 alpha = 1.06085 ; beta = -0.14211 ; rr0 = 2.000 ; saa = 48 440 elif 200 <= ll < 210: 441 alpha = 1.21333 ; beta = -0.23225 ; rr0 = 2.000 ; saa = 57 442 elif 210 <= ll < 220: 443 alpha = 0.58326 ; beta = -0.06097 ; rr0 = 2.000 ; saa = 30 ; gamma = 0.36962 444 elif 220 <= ll < 230: 445 alpha = 0.74200 ; beta = -0.19293 ; rr0 = 1.923 ; saa = 48 ; gamma = 0.07459 446 elif 230 <= ll < 240: 447 alpha = 0.67520 ; beta = -0.21041 ; rr0 = 1.604 ; saa = 49 ; gamma = 0.16602 448 elif 240 <= ll < 250: 449 alpha = 0.62609 ; beta = -0.25312 ; rr0 = 1.237 ; saa = 73 ; gamma = 0.14437 450 elif 250 <= ll < 260: 451 alpha = 0.61415 ; beta = -0.13788 ; rr0 = 2.000 ; saa = 42 ; gamma = 0.26859 452 elif 260 <= ll < 270: 453 alpha = 0.58108 ; beta = 0.01195 ; rr0 = 2.000 ; saa = 40 ; gamma = 0.07661 454 elif 270 <= ll < 280: 455 alpha = 0.68352 ; beta = -0.10743 ; rr0 = 2.000 ; saa = 50 ; gamma = 0.00849 456 elif 280 <= ll < 290: 457 alpha = 0.61747 ; beta = 0.02675 ; rr0 = 2,000 ; saa = 49 458 elif 290 <= ll < 300: 459 alpha = 0.06827 ; beta = -0.26290 ; rr0 = 2.000 ; saa = 44 460 elif 300 <= ll < 310: 461 alpha = 1.53631 ; beta = -0.36833 ; rr0 = 2.000 ; saa = 37 ; gamma = 0.02960 462 elif 210 <= ll < 320: 463 alpha = 1.94300 ; beta = -0.71445 ; rr0 = 1.360 ; saa = 36 ; gamma = 0.15643 464 elif 320 <= ll < 330: 465 alpha = 1.22185 ; beta = -0.40185 ; rr0 = 1.520 ; saa = 48 ; gamma = 0.07354 466 elif 330 <= ll < 340: 467 alpha = 1.79429 ; beta = -0.48657 ; rr0 = 1.844 ; saa = 43 468 elif 340 <= ll < 350: 469 alpha = 2.29545 ; beta = -0.84096 ; rr0 = 1.365 ; saa = 32 470 elif 350 <= ll < 360: 471 alpha = 2.07408 ; beta = -0.64745 ; rr0 = 1.602 ; saa = 36 ; gamma = 0.12750 472 473 474 if 5 < bb < 15: 475 gamma = 0 476 if 0 <= ll < 10: 477 alpha = 2.94213 ; beta = -2.09158 ; rr0 = 0.703 ; saa = 41 ; gamma = 0.05490 478 elif 10 <= ll < 30: 479 alpha = 3.04627 ; beta = 7.71159 ; rr0 = 0.355 ; saa = 45 480 elif 30 <= ll < 40: 481 alpha = 3.78033 ; beta = -3.91956 ; rr0 = 0.482 ; saa = 42 482 elif 40 <= ll < 50: 483 alpha = 2.18119 ; beta = -2.4050 ; rr0 = 0.453 ; saa = 27 484 elif 50 <= ll < 60: 485 alpha = 1.45372 ; beta = -0.49522 ; rr0 = 1.152 ; saa = 31 486 elif 60 <= ll < 70: 487 alpha = 1.05051 ; beta = -1.01704 ; rr0 = 0.516 ; saa = 2 488 elif 70 <= ll < 80: 489 alpha = 0.48416 ; beta = -0.27182 ; rr0 = 0.891 ; saa = 94 ; gamma = 0.08639 490 elif 80 <= ll < 90: 491 alpha = 0.61963 ; beta = 0.41697 ; rr0 = 1.152 ; saa = 35 ; gamma = 0.47171 492 elif 90 <= ll < 100: 493 alpha = 4.40348 ; beta = -2.95011 ; rr0 = 0.745 ; saa = 52 494 elif 100 <= ll < 120: 495 alpha = 2.50938 ; beta = -0.56541 ; rr0 = 1.152 ; saa = 27 496 elif 120 <= ll < 130: 497 alpha = 0.44180 ; beta = 1.58923 ; rr0 = 0.949 ; saa = 4 498 elif 130 <= ll < 140: 499 alpha = 3.96081 ; beta = -3.37374 ; rr0 = 0.587 ; saa = 40 ; gamma = 0.34109 500 elif 140 <= ll < 160: 501 alpha = 2.53335 ; beta = -0.40541 ; rr0 = 1.152 ; saa = 38 502 elif 160 <= ll < 170: 503 alpha = 2.03760 ; beta = -0.66136 ; rr0 = 1.152 ; saa = 23 504 elif 170 <= ll < 200: 505 alpha = 1.06946 ; beta = -0.87395 ; rr0 = 0.612 ; saa = 29 ; gamma = 0.29230 506 elif 200 <= ll < 210: 507 alpha = 0.86348 ; beta = -0.65870 ; rr0 = 0.655 ; saa = 79 ; gamma = 0.09089 508 elif 210 <= ll < 230: 509 alpha = 0.30117 ; beta = -0.16136 ; rr0 = 0.933 ; saa = 17 ; gamma = 0.07495 510 elif 230 <= ll < 240: 511 alpha = 0.75171 ; beta = -0.57143 ; rr0 = 0.658 ; saa = 12 ; gamma = 0.00534 512 elif 240 <= ll < 250: 513 alpha = 1.97427 ; beta = -2.02654 ; rr0 = 0.487 ; saa = 67 514 elif 250 <= ll < 260: 515 alpha = 1.25208 ; beta = -1.47763 ; rr0 = 0.424 ; saa = 19 ; gamma = 0.09089 516 elif 260 <= ll < 270: 517 alpha = 0.89448 ; beta = -0.43870 ; rr0 = 1.019 ; saa = 5 518 elif 270 <= ll < 280: 519 alpha = 0.81141 ; beta = -0.51001 ; rr0 = 0.795 ; saa = 27 ; gamma = 0.03505 520 elif 280 <= ll < 290: 521 alpha = 0.83781 ; beta = -0.44138 ; rr0 = 0.949 ; saa = 50 ; gamma = 0.02820 522 elif 290 <= ll < 300: 523 alpha = 1.10600 ; beta = -0.86263 ; rr0 = 0.641 ; saa = 28 ; gamma = 0.03402 524 elif 300 <= ll < 310: 525 alpha = 1.37040 ; beta = -1.02779 ; rr0 = 0.667 ; saa = 28 ; gamma = 0.05608 526 elif 310 <= ll < 320: 527 alpha = 1.77590 ; beta = -1.26951 ; rr0 = 0.699 ; saa = 37 ; gamma = 0.06972 528 elif 320 <= ll < 330: 529 alpha = 1.20865 ; beta = -0.70679 ; rr0 = 0.855 ; saa = 35 ; gamma = 0.02902 530 elif 330 <= ll < 340: 531 alpha = 2.28830 ; beta = -1.71890 ; rr0 = 0.666 ; saa = 42 ; gamma = 0.22887 532 elif 340 <= ll < 350: 533 alpha = 3.26278 ; beta = -0.94181 ; rr0 = 1.152 ; saa = 38 534 elif 350 <= ll < 360: 535 alpha = 2.58100 ; beta = -1.69237 ; rr0 = 0.763 ; saa = 53 536 537 if 15 < bb < 30: 538 gamma = 0 539 if 0 <= ll < 20: 540 alpha = 6.23279 ; beta = -10.30384 ; rr0 = 0.302 ; saa = 42 541 elif 20 <= ll < 40: 542 alpha = -4.47693 ; beta = -7.28366 ; rr0 = 0.307 ; saa = 29 543 elif 40 <= ll < 60 : 544 alpha = 1.22938 ; beta = -1.19030 ; rr0 = 0.516 ; saa = 5 545 elif 60 <= ll < 80 : 546 alpha = 0.84291 ; beta = -1.59338 ; rr0 = 0.265 ; saa = 4 547 elif 80 <= ll < 100 : 548 alpha = 0.23996 ; beta = 0.06304 ; rr0 = 0.523 ; saa = 32 549 elif 100 <= ll < 140 : 550 alpha = 0.40062 ; beta = -1.75628 ; rr0 = 0.114 ; saa = 16 551 elif 140 <= ll < 180 : 552 alpha = 0.56898 ; beta = -0.53331 ; rr0 = 0.523 ; saa = 41 553 elif 180 <= ll < 200 : 554 alpha = -0.95721 ; beta = 11.6917 ; rr0 = 0.240 ; saa = 2 555 elif 200 <= ll < 220 : 556 alpha = -0.19051 ; beta = 1.45670 ; rr0 = 0.376 ; saa = 1 557 elif 220 <= ll < 240 : 558 alpha = 2.31305 ; beta = -7.82531 ; rr0 = 0.148 ; saa = 95 559 elif 240 <= ll < 260: 560 alpha = 1.39169 ; beta = -1.72984 ; rr0 = 0.402 ; saa = 6 561 elif 260 <= ll < 260: 562 alpha = 1.59418 ; beta = -1.28296 ; rr0 = 0.523 ; saa = 36 563 elif 280 <= ll < 300 : 564 alpha = 1.57082 ; beta = -197295 ; rr0 = 0.398 ; saa = 10 565 elif 300 <= ll < 320 : 566 alpha = 1.95998 ; beta = -3.26159 ; rr0 = 0.300 ; saa = 11 567 elif 320 <= ll < 340: 568 alpha = 2.59567 ; beta = -4.84133 ; rr0 = 0.168 ; saa = 37 569 elif 340 <= ll < 360: 570 alpha = 5.30273 ; beta = -7.43033 ; rr0 = 0.357 ; saa = 37 571 572 if 30 < bb < 45: 573 gamma = 0 574 if 0 <= ll < 20: 575 alpha = 2.93960 ; beta = -6.48019 ; rr0 = 0.227 ; saa = 77 576 elif 20 <= ll < 50: 577 alpha = 1.65864 ; beta = -9.99317 ; rr0 = 0.083 ; saa = 99 578 elif 50 <= ll < 80: 579 alpha = 1.71831 ; beta = -7.25286 ; rr0 = 0.118 ; saa = 28 580 elif 80 <= ll < 110: 581 alpha = 1.33617 ; beta = -10.39799 ; rr0 = 0.064 ; saa = 99 582 elif 110 <= ll < 160: 583 alpha = -0.31330 ; beta = 1.35622 ; rr0 = 0.329 ; saa = 24 584 elif 160 <= ll < 190: 585 alpha = 1.51984 ; beta = -8.69502 ; rr0 = 0.087 ; saa = 99 586 elif 190 <= ll < 220: 587 alpha = -0.50758 ; beta = 4.73320 ; rr0 = 0.250 ; saa = 78 588 elif 220 <= ll < 250: 589 alpha = 1.25864 ; beta = -12.59627 ; rr0 = 0.050 ; saa = 70 590 elif 250 <= ll < 280: 591 alpha = 1.54243 ; beta = -3.75065 ; rr0 = 0.205 ; saa = 10 592 elif 280 <= ll < 320: 593 alpha = 2.72258 ; beta = -7.47806 ; rr0 = 0.182 ; saa = 5 594 elif 320 <= ll < 340: 595 alpha = 2.81435 ; beta = -5.52139 ; rr0 = 0.255 ; saa = 10 596 elif 340 <= ll < 360: 597 alpha = 2.23818 ; beta = 0.81772 ; rr0 = 0.329 ; saa = 19 598 599 if 45 < bb < 60: 600 gamma = 0 601 if 0 <= ll < 60: 602 alpha = 1.38587 ; beta = -9.06536 ; rr0 = 0.076 ; saa = 3 603 elif 60 <= ll < 90: 604 alpha = 2.28570 ; beta = -9.88812 ; rr0 = 0.116 ; saa = 3 605 elif 90 <= ll < 110: 606 alpha = 1.36385 ; beta = -8.10127 ; rr0 = 0.084 ; saa = 4 607 elif 110 <= ll < 170: 608 alpha = 0.05943 ; beta = -1.08126 ; rr0 = 0.027 ; saa = 50 609 elif 170 <= ll < 200: 610 alpha = 1.40171 ; beta = -3.21783 ; rr0 = 0.218 ; saa = 99 611 elif 200 <= ll < 230: 612 alpha = 0.14718 ; beta = 3.92670 ; rr0 = 0.252 ; saa = 14 613 elif 230 <= ll < 290: 614 alpha = 0.57124 ; beta = -4.30242 ; rr0 = 0.066 ; saa = 10 615 elif 290 <= ll < 330: 616 alpha = 3.69891 ; beta = -19.62204 ; rr0 = 0.094 ; saa = 5 617 elif 330 <= ll < 360: 618 alpha = 1.19563 ; beta = -0.45043 ; rr0 = 0.252 ; saa = 9 619 620 if 60 < bb < 90: 621 gamma = 0 622 if 0 <= ll < 30: 623 alpha = 0.69443 ; beta = -0.27600 ; rr0 = 0.153 ; saa = 99 624 elif 30 <= ll < 60: 625 alpha = 1.11811 ; beta = 0.71179 ; rr0 = 0.085 ; saa = 73 626 elif 60 <= ll < 90: 627 alpha = 1.10427 ; beta = -2.37654 ; rr0 = 0.123 ; saa = 99 628 elif 90 <= ll < 120: 629 alpha = -0.42211 ; beta = 5.24037 ; rr0 = 0.184 ; saa = 12 630 elif 120 <= ll < 150: 631 alpha = 0.87576 ; beta = -4.38033 ; rr0 = 0.100 ; saa = 35 632 elif 150 <= ll < 180: 633 alpha = 1.27477 ; beta = -4.98307 ; rr0 = 0.128 ; saa = 72 634 elif 180 <= ll < 210: 635 alpha = 1.19512 ; beta = -6.58464 ; rr0 = 0.091 ; saa = 49 636 elif 210 <= ll < 240: 637 alpha = 0.97581 ; beta = -4.89869 ; rr0 = 0.100 ; saa = 95 638 elif 240 <= ll < 270: 639 alpha = 0.54379 ; beta = -0.84403 ; rr0 = 0.207 ; saa = 35 640 elif 270 <= ll < 300: 641 alpha = 0.85054 ; beta = 13.01249 ; rr0 = 0.126 ; saa = 39 642 elif 300 <= ll < 330: 643 alpha = 0.74347 ; beta = 1.39825 ; rr0 = 0.207 ; saa = 10 644 elif 330 <= ll < 360: 645 alpha = 0.77310 ; beta = -4.45005 ; rr0 = 0.087 ; saa = 16 646 647 return alpha, beta, gamma, rr0, saa

648

649 #} 650 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 651 #{ Marshall 3D extinction model 652 # (Marshall et al. (2006) published in Astronomy and Astrophysics, 653 # Volume 453, Issue 2, July II 2006, pp.635-651 654 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 655 656 @memoized 657 -def get_marshall_data():

658 """ 659 Read in the Marshall data 660 """ 661 662 #data_ma, units_ma, comments_ma = vizier.search("J/A+A/453/635") 663 #filen = config.get_datafile('catalogs','extinction_marshall.tsv') 664 filen = os.path.join(cc.path.ivsdata,'catalogs','extinction_marshall.tsv') 665 data_ma, units_ma, comments_ma = vizier.tsv2recarray(filen) 666 return data_ma, units_ma, comments_ma

667

668 -def findext_marshall(ll, bb, distance=None, redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):

669 """ 670 Find the V-band extinction according to the reddening model of 671 Marshall et al. (2006) published in Astronomy and Astrophysics, 672 Volume 453, Issue 2, July II 2006, pp.635-651 673 674 The band in which the extinction is calculated is actually optional, and given 675 with the keyword C{norm}. If you set C{norm='Av'}, you get visual extinction 676 (in JOHNSON.V) band, if you set C{norm='Ak'}, you get infrared extinction 677 (in JOHNSON.K). If you want the extinction in different passbands, you can 678 set them via C{norm='GENEVA.V'}. So, C{norm='Av'} is equivalent to setting 679 C{norm='JOHNSON.V'}. 680 681 Example usage: 682 683 1. Find the extinction for a star at galactic longitude 10.2 and 684 latitude 9.0. If the distance is not given, we return the 685 complete extinction along the line of sight (i.e. put the star somewhere 686 out of the galaxy). 687 688 >>> lng = 10.2 689 >>> lat = 9.0 690 >>> ak = findext_marshall(lng, lat, norm='Ak') 691 >>> print("Ak at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, ak)) 692 Ak at lng = 10.20, lat = 9.00 is 0.15 magnitude 693 694 2. Find the extinction for a star at galactic lattitude 271.05 and 695 longitude -4.93 and a distance of 144.65 parsecs 696 697 >>> lng = 271.05 698 >>> lat = -4.93 699 >>> dd = 144.65 700 >>> ak = findext_marshall(lng, lat, distance = dd, norm='Ak') 701 >>> print("Ak at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, ak)) 702 Ak at lng = 271.05, lat = -4.93 and distance = 144.65 parsecs is 0.02 magnitude 703 704 3. The extinction for galactic longitude 10.2 and latitude 59.0 can not 705 be calculated using the Marshall models, because they are only valid at 706 certain lng and lat (0 < lng < 100 or 260 < lng < 360, -10 < lat < 10) 707 708 >>> lng = 10.2 709 >>> lat = 59.0 710 >>> ak = findext_marshall(lng, lat, norm='Ak') 711 >>> print(ak) 712 None 713 714 715 @param ll: Galactic Longitude (in degrees) should be between 0 and 100 or 260 and 360 degrees 716 @type ll: float 717 @param bb: Galactic Lattitude (in degrees) should be between -10 and 10 degrees 718 @type bb: float 719 @param distance: Distance to the source (in parsecs) 720 @type distance: float 721 @param redlaw: the used reddening law (standard: 'cardelli1989') 722 @type redlaw: str 723 @param Rv: Av/E(B-V) (standard: 3.1) 724 @type Rv: float 725 @return: The extinction in V-band (or norm) 726 @rtype: float 727 """ 728 729 # get Marshall data 730 data_ma, units_ma, comments_ma = get_marshall_data() 731 732 # Check validity of the coordinates 733 if (ll > 100. and ll < 260.) or ll < 0 or ll > 360: 734 Av = None 735 logger.error("Galactic longitude invalid") 736 return Av 737 elif bb > 10. or bb < -10.: 738 Av = None 739 logger.error("Galactic lattitude invalid") 740 return Av 741 742 # Find the galactic lattitude and longitude of the model, closest to your star 743 dist = np.sqrt((data_ma.GLAT - bb)**2. + (data_ma.GLON - ll)**2.) 744 kma = np.where(dist == dist.min()) 745 746 if min(dist) > .5: 747 logger.error("Could not find a good model value") 748 return Av 749 750 # find the correct index for the distance 751 nb = data_ma.nb[kma] 752 rr = [] ; e_rr = [] ; ext = [] ; e_ext = [] 753 for i in np.arange(nb)+1.: 754 rr.append(data_ma["r%i"%i][kma]) 755 e_rr.append(data_ma["e_r%i"%i][kma]) 756 ext.append(data_ma["ext%i"%i][kma]) 757 e_ext.append(data_ma["e_ext%i"%i][kma]) 758 rr = np.array(rr)[:,0] ; e_rr = np.array(e_rr)[:,0] ; ext = np.array(ext)[:,0] ; e_ext = np.array(e_ext)[:,0] 759 760 # Interpolate linearly in distance. If beyond furthest bin, keep that value. 761 if distance is None: 762 logger.info("No distance given") 763 dd = max(rr) 764 ak = ext[-1] 765 else: 766 dd = distance/1e3 767 768 if dd < min(rr): 769 logger.info("Distance below distance to first bin") 770 ak = (dd/rr[0])*ext[0] 771 elif dd > max(rr): 772 logger.info("Distance more than distance to last bin") 773 ak = ext[-1] 774 dd = max(rr) 775 else: 776 logger.info("Interpolating linearly") 777 ak = sc.interp(dd, rr, ext) 778 779 #-- Marshall is standard in Ak, but you can change this: 780 #redwave, redflux = get_law(redlaw,Rv=Rv,wave_units='micron',norm='Av', wave=array([0.54,2.22])) 781 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.K'],**kwargs) 782 return ak/redflux[0]

783 784 #} 785 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 786 #{ Drimmel 3D extinction model presented in Drimmel et al. in 787 # Astronomy and Astrophysics, v.409, p.205-215 (2003) and 788 # Proceedings of the Gaia Symposium "The Three-Dimensional 789 # Universe with Gaia" 790 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 791 792 # DISCLAIMER: This software is based on the IDL scripts written at the 793 # Cosmology Data Analysis Center in support of the Cosmic Background Explorer 794 # (COBE) Project under NASA contract number NAS5-30750. 795 # This software may be used, copied, modified or redistributed so long as it 796 # is not sold and this disclaimer is distributed along with the software. If 797 # you modify the software please indicate your modifications in a prominent 798 # place in the source code. All routines are provided "as is" without any 799 # express or implied warranties whatsoever. All routines are distributed 800 # without guarantee of support. If errors are found in this code it is 801 # requested that you contact us by sending email to the address below to 802 # report the errors but we make no claims regarding timely fixes. This 803 # software has been used for analysis of COBE data but has not been validated 804 # and has not been used to create validated data sets of any type. 805 # Please send bug reports to CGIS@ZWICKY.GSFC.NASA.GOV. 806 807 fn_base = os.path.join(cc.path.ivsdata,'drimmel') 808 # avdisk = pf.getdata(config.get_datafile('drimmel',"avdisk.fits" )) 809 # avloc = pf.getdata(config.get_datafile('drimmel',"avloc.fits" )) 810 # avspir = pf.getdata(config.get_datafile('drimmel',"avspir.fits" )) 811 # avdloc = pf.getdata(config.get_datafile('drimmel',"avdloc.fits" )) 812 # avori = pf.getdata(config.get_datafile('drimmel',"avori.fits" )) 813 # coordinates = pf.getdata(config.get_datafile('drimmel',"coordinates.fits" )) 814 # avgrid = pf.getdata(config.get_datafile('drimmel',"avgrid.fits" )) 815 # avori2 = pf.getdata(config.get_datafile('drimmel',"avori2.fits" )) 816 # rf_allsky = pf.getdata(config.get_datafile('drimmel',"rf_allsky.fits" )) 817 avdisk = pf.getdata(os.path.join(fn_base,"avdisk.fits" )) 818 avloc = pf.getdata(os.path.join(fn_base,"avloc.fits" )) 819 avspir = pf.getdata(os.path.join(fn_base,"avspir.fits" )) 820 avdloc = pf.getdata(os.path.join(fn_base,"avdloc.fits" )) 821 avori = pf.getdata(os.path.join(fn_base,"avori.fits" )) 822 coordinates = pf.getdata(os.path.join(fn_base,"coordinates.fits" )) 823 avgrid = pf.getdata(os.path.join(fn_base,"avgrid.fits" )) 824 avori2 = pf.getdata(os.path.join(fn_base,"avori2.fits" )) 825 rf_allsky = pf.getdata(os.path.join(fn_base,"rf_allsky.fits" )) 826 glat = rf_allsky.glat 827 glng = rf_allsky.glng 828 ncomp = rf_allsky.ncomp 829 pnum = rf_allsky.pnum 830 rfac = rf_allsky.rfac 831 832 g2e = array([[-0.054882486, -0.993821033, -0.096476249], [0.494116468, -0.110993846, 0.862281440], [-0.867661702, -0.000346354, 0.497154957]])

833 834 -def findext_drimmel(lng, lat, distance=None, rescaling=True, 835 redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):

836 """ 837 Procedure to retrieve the absorption in V from three-dimensional 838 grids, based on the Galactic dust distribution of Drimmel & Spergel. 839 840 Example usage: 841 842 1. Find the extinction for a star at galactic longitude 10.2 and 843 latitude 9.0. If the distance is not given, we return the 844 complete extinction along the line of sight (i.e. put the star somewhere 845 out of the galaxy). 846 847 >>> lng = 10.2 848 >>> lat = 9.0 849 >>> av = findext_drimmel(lng, lat) 850 >>> print("Av at lng = %.2f, lat = %.2f is %.2f magnitude" %(lng, lat, av)) 851 Av at lng = 10.20, lat = 9.00 is 1.24 magnitude 852 853 2. Find the extinction for a star at galactic lattitude 107.05 and 854 longitude -34.93 and a distance of 144.65 parsecs 855 856 >>> lng = 271.05 857 >>> lat = -4.93 858 >>> dd = 144.65 859 >>> ak = findext_marshall(lng, lat, distance = dd) 860 >>> print("Ak at lng = %.2f, lat = %.2f and distance = %.2f parsecs is %.2f magnitude" %(lng, lat, dd, ak)) 861 Ak at lng = 271.05, lat = -4.93 and distance = 144.65 parsecs is 0.02 magnitude 862 863 864 @param lng: Galactic Longitude (in degrees) 865 @type lng: float 866 @param lat: Galactic Lattitude (in degrees) 867 @type lat: float 868 @param distance: Distance to the source (in parsecs) 869 @type distance: float 870 @param rescaling: Rescaling needed or not? 871 @type rescaling: boolean 872 @return: extinction in V band with/without rescaling 873 @rtype: float 874 """ 875 # Constants 876 deg2rad = pi/180. # convert degrees to rads 877 nsky = 393216 # number of COBE pixels 878 879 # Sun's coordinates (get from dprms) 880 xsun = -8.0 881 zsun = 0.015 882 883 # if distance is not given, make it large and put it in kiloparsec 884 if distance is None: 885 d = 1e10 886 else: 887 d = distance/1e3 888 889 # build skymaps of rescaling parameters for each component 890 dfac = ones(nsky) 891 sfac = ones(nsky) 892 lfac = ones(nsky) 893 indx = where(ncomp == 1) 894 dfac[indx] = rfac[indx] 895 indx = where(ncomp == 2) 896 sfac[indx] = rfac[indx] 897 indx = where(ncomp == 3) 898 lfac[indx] = rfac[indx] 899 900 # define abs 901 num = array(d).size 902 out = zeros(num) 903 avloc = zeros(num) 904 abspir = zeros(num) 905 absdisk = zeros(num) 906 907 l = lng*deg2rad # [radians] 908 b = lat*deg2rad # [radians] 909 910 # Now for UIDL code: 911 # -find the index of the corresponding COBE pixel 912 # dimensions of sixpack = 768 x 512 (= 393216) 913 vectoarr = arange(768*512) 914 vectoarr.shape = (512,768) # necessary because I reduced the arrays to vectors. 915 vectoarr = vectoarr.transpose() 916 res = 9 917 pxindex = _ll2pix(lng, lat, res) 918 xout, yout = _pix2xy(pxindex, res, sixpack=True) 919 tblindex = vectoarr[xout, yout] # calculate the maximum distance in the grid 920 dmax = ones(num)*100. 921 if b != 0.: 922 dmax = .49999/abs(sin(b)) - zsun/sin(b) 923 if cos(l) != 0.: 924 dmax = min([dmax,(14.9999/abs(cos(l)) - xsun/cos(l))]) 925 if sin(l) != 0.: 926 dmax = min([dmax, 14.9999/abs(sin(l))]) 927 928 # replace distance with dmax when greater 929 r = array([d]) 930 indx = where(array([d]) > dmax)[0] 931 if len(indx) > 0: 932 r[indx] = dmax 933 934 # heliocentric cartesian coordinates 935 x = array([r*cos(b)*cos(l)]) 936 y = array([r*cos(b)*sin(l)]) 937 z = array([r*sin(b) + zsun]) 938 939 # for stars in Solar neighborhood 940 i = where(logical_and(abs(x) < 1.,abs(y) < 2.))[0] 941 ni = len(i) 942 j = where(logical_or(abs(x) >= 1., abs(y) >= 2.))[0] 943 nj = len(j) 944 945 if ni > 0: 946 # define the local grid 947 dx = 0.02 948 dy = 0.02 949 dz = 0.02 950 nx = 101 951 ny = 201 952 nz = 51 953 954 # grid indices 955 xi = x[i]/dx + float(nx - 1)/2. 956 yj = y[i]/dy + float(ny - 1)/2. 957 zk = z[i]/dz + float(nz - 1)/2. 958 959 # interpolate 960 avloc[i] = _trint(avori2, xi, yj, zk, missing=0.) 961 962 # for stars in Solar neighborhood 963 k = where(logical_and(abs(x) < 0.75, abs(y) < 0.75))[0] 964 nk = len(k) 965 m = where(logical_or(abs(x) >= 0.75, abs(y) >= 0.75))[0] 966 nm = len(m) 967 968 if nk > 0: 969 970 # define the local grid 971 dx = 0.05 972 dy = 0.05 973 dz = 0.02 974 nx = 31 975 ny = 31 976 nz = 51 977 978 # grid indices 979 xi = x[k]/dx + float(nx - 1)/2. 980 yj = y[k]/dy + float(ny - 1)/2. 981 zk = z[k]/dz + float(nz - 1)/2. 982 983 # trilinear_interpolate 984 absdisk[k] = _trint(avdloc,xi,yj,zk,missing=0.) 985 986 # galacto-centric cartesian 987 x = x + xsun 988 989 #larger orion arm grid 990 if nj > 0: 991 # calculate the allowed maximum distance for larger orion grid 992 dmax = 100. 993 if b != 0.: 994 dmax = .49999/abs(sin(b)) - zsun/sin(b) 995 if cos(l) > 0.: 996 dmax = min([dmax, 2.374999/abs(cos(l))]) 997 if cos(l) < 0.: 998 dmax = min([dmax, 1.374999/abs(cos(l))]) 999 if sin(l) != 0.: 1000 dmax = min([dmax, 3.749999/abs(sin(l))]) 1001 1002 # replace distance with dmax when greater 1003 r1 = array([d]) 1004 indx = where(array([d]) >= dmax)[0] 1005 n = len(indx) 1006 if n > 0: 1007 r1[indx] = dmax 1008 1009 # galactocentric centric cartesian coordinates 1010 x1 = array([r1*cos(b)*cos(l) + xsun]) 1011 y1 = array([r1*cos(b)*sin(l)]) 1012 z1 = array([r1*sin(b) + zsun]) 1013 1014 # define the grid 1015 dx = 0.05 1016 dy = 0.05 1017 dz = 0.02 1018 nx = 76 1019 ny = 151 1020 nz = 51 1021 1022 # grid indices 1023 xi = x1[j]/dx + 2.5*float(nx - 1) 1024 yj = y1[j]/dy + float(ny - 1)/2. 1025 zk = z1[j]/dz + float(nz - 1)/2. 1026 1027 # trilinear_interpolate 1028 avloc[j] = _trint(avori,xi,yj,zk,missing=0.) 1029 1030 # define the grid 1031 dx = 0.2 1032 dy = 0.2 1033 dz = 0.02 1034 nx = 151 1035 ny = 151 1036 nz = 51 1037 1038 # grid indices 1039 xi = x/dx + float(nx - 1)/2. 1040 yj = y/dy + float(ny - 1)/2. 1041 zk = z/dz + float(nz - 1)/2. 1042 1043 # trilinear_interpolate 1044 abspir = _trint(avspir,xi,yj,zk,missing=0.) 1045 if nm > 0: 1046 absdisk[m] = _trint(avdisk,xi[m],yj[m],zk[m],missing=0.) 1047 1048 # apply rescaling factors or not 1049 if rescaling: 1050 out = (dfac[tblindex]*absdisk + sfac[tblindex]*abspir + lfac[tblindex]*avloc).flatten() 1051 else: 1052 out = (absdisk + abspir + avloc).flatten() 1053 1054 #-- Drimmel is standard in Av, but you can change this: 1055 # In case of norm=Av, the conversion factor is just 1. 1056 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'],**kwargs) 1057 1058 return(out/redflux[0])

1059

1060 -def _pix2xy(pixel, resolution, sixpack=0):

1061 """ 1062 _pix2xy creates a raster image (sky cube or face) given a pixelindex and a 1063 resolution The data array can be either a vector or two-dimensional array. 1064 In the latter case, the data for each raster image can be stored in either 1065 the columns or rows. The procedure also returns the x and y raster 1066 coordinates of each pixel. Only right oriented, ecliptic coordinate images 1067 are built. 1068 SMOLDERS SEAL OF APPROVAL 1069 """ 1070 # reform pixel to get rid of all length-1 dimensions 1071 pixel = array(pixel) 1072 resolution = int(resolution) 1073 newshape = [] 1074 for s in pixel.shape: 1075 if s > 1: 1076 newshape.append() 1077 pixel.shape = newshape 1078 pix_sz = (pixel) 1079 1080 if max(pixel) > 6*4**(resolution-1): 1081 raise('Maximum pixel number too large for resolution') 1082 1083 # set up flag values for RASTR 1084 data = [-1] 1085 face = -1 1086 bad_pixval = 0.0 1087 # call rasterization routine 1088 xout, yout = _rastr(pixel,resolution) 1089 return(xout, yout)

1090

1091 -def _rastr(pixel,resolution):

1092 """ 1093 SMOLDERS SEAL OF APPROVAL 1094 """ 1095 pixel = array(pixel) 1096 n = pixel.size 1097 i0 = 3 1098 j0 = 2 1099 offx = [0,0,1,2,2,1] 1100 offy = [1,0,0,0,1,1] 1101 fij = _pix2fij(pixel,resolution) 1102 cube_side = 2**(resolution-1) 1103 lenc = i0*cube_side 1104 if n > 1: 1105 x_out = offx[fij[0,:]] * cube_side + fij[1,:] 1106 x_out = lenc - (x_out+1) 1107 y_out = offy[fij[0,:]] * cube_side + fij[2,:] 1108 else: 1109 x_out = offx[fij[0]] * cube_side + fij[1] 1110 x_out = lenc - (x_out+1) 1111 y_out = offy[fij[0]] * cube_side + fij[2] 1112 return(x_out, y_out)

1113

1114 -def _pix2fij(pixel,resolution):

1115 """ 1116 This function takes an n-element pixel array and generates an n by 3 element 1117 array containing the corresponding face, column, and row number (the latter 1118 two within the face). 1119 SMOLDERS SEAL OF APPROVAL 1120 """ 1121 # get number of pixels 1122 pixel = array(pixel) 1123 n = pixel.size 1124 output = array(zeros((3,n)), int) 1125 res1 = resolution - 1 1126 num_pix_face = 4**res1 1127 # 1128 face = pixel/num_pix_face 1129 fpix = pixel-num_pix_face*face 1130 output[0,:] = face 1131 pow_2 = 2**arange(16) 1132 ii = array(zeros(n), int) 1133 jj = array(zeros(n), int) 1134 # 1135 if n > 1: 1136 for i in arange(n): 1137 for bit in arange(res1): 1138 ii[i] = ii[i] | (pow_2[bit]*(1 & fpix)) 1139 fpix = fpix >> 1 1140 jj[i] = jj[i] | (pow_2[bit]*(1 & fpix)) 1141 fpix = fpix >> 1 1142 else: 1143 for bit in arange(res1): 1144 ii = ii | (pow_2[bit]*(1 & fpix)) 1145 fpix = fpix >> 1 1146 jj = jj | (pow_2[bit]*(1 & fpix)) 1147 fpix = fpix >> 1 1148 output[1,:] = ii 1149 output[2,:] = jj 1150 return output

1151

1152 -def _incube(alpha, beta):

1153 gstar = 1.37484847732 1154 g = -0.13161671474 1155 m = 0.004869491981 1156 w1 = -0.159596235474 1157 c00 = 0.141189631152 1158 c10 = 0.0809701286525 1159 c01 = -0.281528535557 1160 c11 = 0.15384112876 1161 c20 = -0.178251207466 1162 c02 = 0.106959469314 1163 d0 = 0.0759196200467 1164 d1 = -0.0217762490699 1165 r0 = 0.577350269 1166 aa = alpha**2. 1167 bb = beta**2 1168 a4 = aa**2 1169 b4 = bb**2 1170 onmaa = 1.-aa 1171 onmbb = 1.-bb 1172 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)))) 1173 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)))) 1174 return x, y

1175

1176 -def _ll2uv(lng, lat):

1177 """ 1178 Convert longitude, latitude to unit vectors 1179 SMOLDERS SEAL OF APPROVAL 1180 1181 @param lng : galactic longitude 1182 @type lng : float 1183 @param lat : galactic lattitude 1184 @type lat : float 1185 @return : unitvector 1186 @rtype : ndarray 1187 """ 1188 vector = zeros(3) 1189 d2r = pi/180 1190 lng = lng * d2r 1191 lat = lat * d2r 1192 vector[0] = cos(lat) * cos(lng) 1193 vector[1] = cos(lat) * sin(lng) 1194 vector[2] = sin(lat) 1195 return vector

1196

1197 -def _galvec2eclvec(in_uvec):

1198 """ 1199 Unitvector for galactic coordinates to unitvector for ecliptic coordinates 1200 SMOLDERS SEAL OF APPROVAL 1201 """ 1202 out_uvec = dot(in_uvec, g2e) 1203 return out_uvec

1204

1205 -def _axisxy(vector):

1206 """ 1207 converts unitvector into nface number (0-5) and X,Y in range 0-1 1208 SMOLDERS SEAL OF APPROVAL 1209 """ 1210 vector = array(vector) 1211 vs = vector.shape 1212 if len(vs) > 1: 1213 vec0 = array(vector[:,0]) 1214 vec1 = array(vector[:,1]) 1215 vec2 = array(vector[:,2]) 1216 else: 1217 vec0 = array([vector[0]]) 1218 vec1 = array([vector[1]]) 1219 vec2 = array([vector[2]]) 1220 abs_yx = abs(vec1/vec0) 1221 abs_zx = abs(vec2/vec0) 1222 abs_zy = abs(vec2/vec1) 1223 # 1224 nface = zeros(len(vec0)) 1225 for i in arange(len(vec0)): 1226 nface[i] = (0 * (abs_zx[i] >= 1 and abs_zy[i] >= 1 and vec2[i] >= 0) + 1227 5 * (abs_zx[i] >= 1 and abs_zy[i] >= 1 and vec2[i] < 0) + 1228 1 * (abs_zx[i] < 1 and abs_yx[i] < 1 and vec0[i] >= 0) + 1229 3 * (abs_zx[i] < 1 and abs_yx[i] < 1 and vec0[i] < 0) + 1230 2 * (abs_zy[i] < 1 and abs_yx[i] >= 1 and vec1[i] >= 0) + 1231 4 * (abs_zy[i] < 1 and abs_yx[i] >= 1 and vec1[i] < 0)) 1232 # 1233 nface_0 = (nface == 0)*1. 1234 nface_1 = (nface == 1)*1. 1235 nface_2 = (nface == 2)*1. 1236 nface_3 = (nface == 3)*1. 1237 nface_4 = (nface == 4)*1. 1238 nface_5 = (nface == 5)*1. 1239 # 1240 eta = (vec2/vec0)*(nface_1 - nface_3) + (vec2/vec1)*(nface_2 - nface_4) - (vec0/vec2)*(nface_0 + nface_5) 1241 xi = (vec1/vec0)*(nface_1 + nface_3) - (vec0/vec1)*(nface_2 + nface_4) + (vec1/vec2) * (nface_0 - nface_5) 1242 x, y = _incube(xi,eta) 1243 x = (x+1.)/2. 1244 y = (y+1.)/2. 1245 return x,y,array(nface, dtype=int)

1246

1247 -def _fij2pix(fij,res):

1248 """ 1249 This function takes an n by 3 element vector containing the face, column, and 1250 row number (the latter two within the face) of a pixel and converts it into an 1251 n-element pixel array for a given resolution. 1252 """ 1253 #get the number of pixels 1254 fij = array(fij) 1255 n = fij.size/3 1256 # generate output pixel and intermediate pixel arrays 1257 pixel = array(zeros(n), dtype=int) 1258 pixel_1 = array(zeros(n), dtype=int) 1259 # get input column and row numbers 1260 if n > 1: 1261 ff = array(fij[0,:], dtype=int) 1262 ii = array(fij[1,:], dtype=int) 1263 jj = array(fij[2,:], dtype=int) 1264 else: 1265 ff = int(fij[0]) 1266 ii = int(fij[1]) 1267 jj = int(fij[2]) 1268 # calculate the number of pixels in a face 1269 num_pix_face = 4**(res-1) 1270 pow_2 = 2**arange(16) 1271 # if col bit set then set corresponding even bit in pixel_l 1272 # if row bit set then set corresponding odd bit in pixel_l 1273 if n > 1: 1274 for i in arange(n): 1275 for bit in arange(res-1): 1276 pixel_1[i] = pixel_1[i] | ((pow_2[bit] & ii[i]) << bit) 1277 pixel_1[i] = pixel_1[i] | ((pow_2[bit] & jj[i]) << (bit+1)) 1278 else: 1279 for bit in arange(res-1): 1280 pixel_1 = pixel_1 | ((pow_2[bit] & ii) << bit) 1281 pixel_1 = pixel_1 | ((pow_2[bit] & jj) << (bit+1)) 1282 # add face number offset 1283 pixel = ff*num_pix_face + pixel_1 1284 return pixel

1285

1286 -def _uv2pix(vector, resolution):

1287 """ 1288 Routine returns pixel number given unit vector pointing to center of pixel 1289 resolution of the cube. 1290 SMOLDERS SEAL OF APPROVAL 1291 """ 1292 two = 2 1293 vector = array(vector) 1294 n_vec = int(vector.size/3) 1295 res1 = resolution - 1 1296 num_pix_side = int(two**res1) 1297 x, y, face = _axisxy(vector) 1298 ia = array(x*num_pix_side, dtype=int) 1299 ib = array([2.**res1 - 1]*n_vec, dtype=int) 1300 ja = array(y*num_pix_side, dtype=int) 1301 jb = array([2.**res1 - 1]*n_vec, dtype=int) 1302 i = zeros(n_vec) 1303 j = zeros(n_vec) 1304 for k in arange(n_vec): 1305 i[k] = min([ia[k],ib[k]]) 1306 j[k] = min([ja[k],jb[k]]) 1307 pixel = _fij2pix(array([face,i,j]),resolution) 1308 return pixel

1309

1310 -def _ll2pix(lng, lat, res=9):

1311 """ 1312 _ll2pix is a python function to convert galactic coordinates to DIRBE pixels 1313 (coorconv([lng,lat], infmt='L', outfmt='P', inco='G', outco='R9')) 1314 1315 Input 1316 @param lng : galactic longitude 1317 @type lng : float 1318 @param lat : galactic lattitude 1319 @type lat : float 1320 @return: output coordinate array 1321 @rtype: ndarray 1322 1323 (Note: The default coordinate system for pixels is ecliptic.) 1324 """ 1325 # COMMON sky_conv_com, e2g,e2q,g2e,g2q,q2e,q2g,load_flag 1326 uvec = _ll2uv(lng, lat) 1327 uvec = _galvec2eclvec(uvec) 1328 output = _uv2pix(uvec,res) 1329 return output

1330

1331 -def _trint(mm,x,y,z,missing=None):

1332 """ 1333 Pixel-based trilinear interpolation, where mm is a 3d data cube. 1334 """ 1335 xmi = array(floor(x),int) 1336 xd = abs(xmi-x) 1337 xma = array(ceil(x),int) 1338 ymi = array(floor(y), int) 1339 yd = abs(ymi-y) 1340 yma = array(ceil(y),int) 1341 zmi = array(floor(z), int) 1342 zd = abs(zmi-z) 1343 zma = array(ceil(z), int) 1344 i1 = mm[zmi,ymi,xmi]*(1.-zd) + mm[zma,ymi,xmi]*(zd) 1345 i2 = mm[zmi,yma,xmi]*(1.-zd) + mm[zma,yma,xmi]*(zd) 1346 j1 = mm[zmi,ymi,xma]*(1.-zd) + mm[zma,ymi,xma]*(zd) 1347 j2 = mm[zmi,yma,xma]*(1.-zd) + mm[zma,yma,xma]*(zd) 1348 w1 = i1*(1-yd) + i2*yd 1349 w2 = j1*(1-yd) + j2*yd 1350 iv = w1*(1.-xd) + w2*xd 1351 return iv

1352

1353 #} 1354 1355 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1356 #{ Schlegel 3D extinction model presented in Schlegel et al. 1357 # The Astrophysical Journal, 500:525È553, 1998 June 20, 1358 # "MAPS OF DUST INFRARED EMISSION FOR USE IN ESTIMATION 1359 # OF REDDENING AND COSMIC MICROWAVE BACKGROUND RADIATION 1360 # FOREGROUNDS" 1361 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1362 1363 @memoized 1364 -def get_schlegel_data_south():

1365 # Read in the Schlegel data of the southern hemisphere 1366 1367 #dustname = config.get_datafile('schlegel',"SFD_dust_4096_sgp.fits") 1368 #maskname = config.get_datafile('schlegel',"SFD_mask_4096_sgp.fits") 1369 dustname = os.path.join(cc.path.ivsdata,'schlegel',"SFD_dust_4096_sgp.fits") 1370 maskname = os.path.join(cc.path.ivsdata,'schlegel',"SFD_mask_4096_sgp.fits") 1371 data = pf.getdata(dustname) 1372 mask = pf.getdata(maskname) 1373 return data, mask

1374

1375 -def get_schlegel_data_north():

1376 # Read in the Schlegel data of the northern hemisphere 1377 #dustname = config.get_datafile('schlegel',"SFD_dust_4096_ngp.fits") 1378 #maskname = config.get_datafile('schlegel',"SFD_mask_4096_ngp.fits") 1379 dustname = os.path.join(cc.path.ivsdata,'schlegel',"SFD_dust_4096_ngp.fits") 1380 maskname = os.path.join(cc.path.ivsdata,'schlegel',"SFD_mask_4096_ngp.fits") 1381 data = pf.getdata(dustname) 1382 mask = pf.getdata(maskname) 1383 return data, mask

1384

1385 -def _lb2xy_schlegel(ll, bb):

1386 """ 1387 Converts coordinates in lattitude and longitude to coordinates 1388 in x and y pixel coordinates. 1389 1390 The x and y coordinates you find with these formulas are not 1391 the coordinates you can read in ds9, but 1 smaller. Hence 1392 (x+1, y+1) are the coordinates you find in DS9. 1393 1394 Input 1395 @param ll : galactic longitude 1396 @type ll : float 1397 @param bb : galactic lattitude 1398 @type bb : float 1399 @return: output coordinate array 1400 @rtype: ndarray 1401 """ 1402 deg2rad = pi/180. # convert degrees to rads 1403 1404 if bb <= 0 : hs = -1. 1405 elif bb > 0: hs = +1. 1406 1407 yy = 2048 * sqrt(1. - hs * sin(bb*deg2rad)) * cos(ll*deg2rad) + 2047.5 1408 xx = -2048 * hs * sqrt(1 - hs * sin(bb*deg2rad)) * sin(ll*deg2rad) + 2047.5 1409 1410 return xx, yy

1411

1412 -def findext_schlegel(ll, bb, distance = None, redlaw='cardelli1989', Rv=3.1, norm='Av',**kwargs):

1413 """ 1414 Get the "Schlegel" extinction at certain longitude and latitude 1415 1416 This function returns the E(B-V) maps of Schlegel et al. (1998), 1417 depending on wether the distance is given or not, the E(B-V) value 1418 is corrected. If distance is set we use the distance-corrected 1419 values: 1420 1421 E(B-V) = E(B-V)_maps * (1 - exp(-10 * r * sin(|b|))) 1422 where E(B-V) is the value to be used and E(B-V)_maps the value 1423 as found with the Schlegel dust maps 1424 1425 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. 1426 1427 ! WARNING: the schlegel maps are not usefull when |b| < 5 degrees ! 1428 """ 1429 deg2rad = pi/180. # convert degrees to rads 1430 dd = distance 1431 if distance is not None: 1432 dd = distance/1.e3 # convert to kpc 1433 1434 1435 # first get the right pixel coordinates 1436 xx, yy = _lb2xy_schlegel(ll,bb) 1437 1438 # read in the right map: 1439 if bb <= 0: 1440 data, mask = get_schlegel_data_south() 1441 elif bb > 0: 1442 data, mask = get_schlegel_data_north() 1443 1444 if abs(bb) < 10.: 1445 logger.warning("Schlegel is not good for lattitudes > 10 degrees") 1446 1447 # the xy-coordinates are: 1448 xl = floor(xx) 1449 yl = floor(yy) 1450 xh = xl + 1. 1451 yh = yl + 1. 1452 1453 # the weights are just the distances to the points 1454 w1 = (xl-xx)**2 + (yl-yy)**2 1455 w2 = (xl-xx)**2 + (yh-yy)**2 1456 w3 = (xh-xx)**2 + (yl-yy)**2 1457 w4 = (xh-xx)**2 + (yh-yy)**2 1458 1459 # the values of these points are: 1460 v1 = data[xl, yl] 1461 v2 = data[xl, yh] 1462 v3 = data[xh, yl] 1463 v4 = data[xh, yh] 1464 f1 = mask[xl, yl] 1465 f2 = mask[xl, yh] 1466 f3 = mask[xh, yl] 1467 f4 = mask[xh, yh] 1468 1469 ebv = (w1*v1 + w2*v2 + w3*v3 + w4*v4) / (w1 + w2 + w3 + w4) 1470 1471 # Check flags at the right pixels 1472 logger.info("flag of pixel 1 is: %i" %f1) 1473 logger.info("flag of pixel 2 is: %i" %f2) 1474 logger.info("flag of pixel 3 is: %i" %f3) 1475 logger.info("flag of pixel 4 is: %i" %f4) 1476 1477 if dd is not None: 1478 ebv = ebv * (1. - exp(-10. * dd * sin(abs(bb*deg2rad)))) 1479 1480 # if Rv is given, by definition we find Av = Ebv*Rv 1481 av = ebv*Rv 1482 1483 #-- Marshall is standard in Ak, but you can change this: 1484 redwave, redflux = get_law(redlaw,Rv=Rv,norm=norm,photbands=['JOHNSON.V'],**kwargs) 1485 1486 return av/redflux[0]

1487 1488 #} 1489 1490 if __name__ == "__main__": 1491 print findext_marshall(10.2,9.) 1492 print findext_marshall(10.2,9.,norm='Ak') 1493

Source Code for Module ComboCode.cc.ivs.sed.extinctionmodels