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 

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

19 """ 

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

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

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

23 return spline 

24 

25 

26def extinction_coefficient(lam, excess_b_v, spline): 

27 """ 

28 Calculate the extinction coefficient as described in 

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

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

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

32 

33 :param lam: Wavelength in micrometre 

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

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

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

37 (1999, DOI: 10.1086/316293) 

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

39 """ 

40 

41 def uv_extinction(x): 

42 """ 

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

44 Fitzpatrick (1999, DOI: 10.1086/316293) 

45 

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

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

48 """ 

49 

50 def uv_curvature(x): 

51 f = np.zeros_like(x) 

52 mask = x >= 5.9 

53 y = x[mask] - 5.9 

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

55 return f 

56 

57 r_v = 3.1 

58 c_1 = 4.50777 - 14.184 / r_v 

59 c_2 = -0.824 + 4.717 / r_v 

60 c_3 = 3.23 

61 c_4 = 0.41 

62 gam = 0.99 

63 x_0 = 4.596 

64 

65 x_sq = x**2 

66 

67 k = ( 

68 c_1 

69 + c_2 * x 

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

71 + c_4 * uv_curvature(x) 

72 ) 

73 

74 return k + r_v 

75 

76 lam_inv = 1 / lam 

77 uv_mask = ( 

78 lam_inv >= 1 / 0.27 

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

80 

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

82 np.newaxis, ... 

83 ] 

84 ext_coeff = ( 

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

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

87 

88 return ext_coeff 

89 

90 

91def construct_intrinsic_spectrum(coeff, templates_amp): 

92 """ 

93 :param templates_amp: (n_templates, n_lambda) 

94 :param coeff: (n_gal, n_templates) 

95 :return: (n_gal, n_lambda) 

96 """ 

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

98 return specrum 

99 

100 

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

102 """ 

103 :param spectrum: (n_gal, n_lambda) 

104 :param lam_obs: (n_lambda,) 

105 :param excess_b_v: (n_gal,) 

106 :param extinction_spline: 

107 """ 

108 spectrum *= 10 ** ( 

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

110 ) 

111 

112 

113def construct_reddened_spectrum( 

114 lam_obs, templates_amp, coeff, excess_b_v, extinction_spline 

115): 

116 """ 

117 :param lam_obs: (n_lambda, ) 

118 :param templates_amp: (n_templates, n_lambda) 

119 :param coeff: (n_gal, n_templates) 

120 :param excess_b_v: (n_gal,) 

121 :param extinction_spline: 

122 :return: (n_gal, n_lambda) 

123 """ 

124 spectrum = construct_intrinsic_spectrum(coeff, templates_amp) 

125 apply_extinction(spectrum, lam_obs, excess_b_v, extinction_spline) 

126 return spectrum