Coverage for src/galsbi/ucat/spectrum_util.py: 100%

1# Copyright (C) 2018 ETH Zurich, Institute for Particle Physics and Astrophysics 

2 

3""" 

4Created on Mar 6, 2018 

5author: Joerg Herbel 

6""" 

7 

8import numpy as np 

9import scipy.interpolate 

10 

11 

12def spline_ext_coeff(): 

13 """ 

14 Perform the cubic spline interpolation described by 

15 Fitzpatrick (1999, DOI: 10.1086/316293) to obtain the extinction law in the optical 

16 and IR region (for wavelengths > 2700 Angstrom) 

17 :return: spline object than can be evaluated at arbitrary inverse wavelengths 

18 """ 

19 spline_x = (0.0, 0.377, 0.820, 1.667, 1.828, 2.141, 2.433, 3.704, 3.846) 

20 spline_y = (0.0, 0.265, 0.829, 2.688, 3.055, 3.806, 4.315, 6.265, 6.591) 

21 spline = scipy.interpolate.InterpolatedUnivariateSpline(spline_x, spline_y, k=3) 

22 return spline 

23 

24 

25def extinction_coefficient(lam, excess_b_v, spline): 

26 """ 

27 Calculate the extinction coefficient as described in 

28 Schlafly et. al. (2011, DOI: 10.1088/0004-637X/737/2/103). The extinction law is the 

29 one given by Fitzpatrick (1999, DOI: 10.1086/316293), the extinction map is the one 

30 presented in Schlegel et. al. (1998, DOI: 10.1086/305772). 

31 :param lam: Wavelength in micrometre 

32 :param excess_b_v: Excess E(B-V) in B-V color from the map provided by Schlegel et. 

33 al. (1998, DOI: 10.1086/305772) 

34 :param spline: Cubic spline for the optical and IR-region according to Fitzpatrick 

35 (1999, DOI: 10.1086/316293) 

36 :return: Extinction coefficient A_lambda evaluated at the input wavelengths. 

37 """ 

38 

39 def uv_extinction(x): 

40 """ 

41 Extinction law for the UV-region (wavelengths < 2700 Angstrom) according to 

42 Fitzpatrick (1999, DOI: 10.1086/316293) 

43 :param x: Inverse wavelength in micrometre^-1 

44 :return: A_lambda/E(B-V) evaluated at input inverse wavelengthss 

45 """ 

46 

47 def uv_curvature(x): 

48 f = np.zeros_like(x) 

49 mask = x >= 5.9 

50 y = x[mask] - 5.9 

51 f[mask] = 0.5392 * y**2 + 0.05644 * y**3 

52 return f 

53 

54 r_v = 3.1 

55 c_1 = 4.50777 - 14.184 / r_v 

56 c_2 = -0.824 + 4.717 / r_v 

57 c_3 = 3.23 

58 c_4 = 0.41 

59 gam = 0.99 

60 x_0 = 4.596 

61 

62 x_sq = x**2 

63 

64 k = ( 

65 c_1 

66 + c_2 * x 

67 + c_3 * x_sq / ((x_sq - x_0**2) ** 2 + x_sq * gam**2) 

68 + c_4 * uv_curvature(x) 

69 ) 

70 

71 return k + r_v 

72 

73 lam_inv = 1 / lam 

74 uv_mask = ( 

75 lam_inv >= 1 / 0.27 

76 ) # UV region is defined as all wavelengths <= 2700 Angstrom 

77 

78 ext_coeff = np.append(uv_extinction(lam_inv[uv_mask]), spline(lam_inv[~uv_mask]))[ 

79 np.newaxis, ... 

80 ] 

81 ext_coeff = ( 

82 ext_coeff * 0.78 * 1.32 * excess_b_v[..., np.newaxis] / spline(1) 

83 ) # See Schlafly et. al. (2011, DOI: 10.1088/0004-637X/737/2/103), eq. (A1) 

84 

85 return ext_coeff 

86 

87 

88def construct_intrinsic_spectrum(coeff, templates_amp): 

89 """ 

90 :param templates_amp: (n_templates, n_lambda) 

91 :param coeff: (n_gal, n_templates) 

92 :return: (n_gal, n_lambda) 

93 """ 

94 specrum = np.sum(coeff[..., np.newaxis] * templates_amp[np.newaxis, ...], axis=1) 

95 return specrum 

96 

97 

98def apply_extinction(spectrum, lam_obs, excess_b_v, extinction_spline): 

99 """ 

100 :param spectrum: (n_gal, n_lambda) 

101 :param lam_obs: (n_lambda,) 

102 :param excess_b_v: (n_gal,) 

103 :param extinction_spline: 

104 """ 

105 spectrum *= 10 ** ( 

106 -2 / 5 * extinction_coefficient(lam_obs, excess_b_v, extinction_spline) 

107 ) 

108 

109 

110def construct_reddened_spectrum( 

111 lam_obs, templates_amp, coeff, excess_b_v, extinction_spline 

112): 

113 """ 

114 :param lam_obs: (n_lambda, ) 

115 :param templates_amp: (n_templates, n_lambda) 

116 :param coeff: (n_gal, n_templates) 

117 :param excess_b_v: (n_gal,) 

118 :param extinction_spline: 

119 :return: (n_gal, n_lambda) 

120 """ 

121 spectrum = construct_intrinsic_spectrum(coeff, templates_amp) 

122 apply_extinction(spectrum, lam_obs, excess_b_v, extinction_spline) 

123 return spectrum