Coverage for src/galsbi/ucat/galaxy_population_models/galaxy_size.py: 94%

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

2 

3""" 

4Created 2021 

5author: Tomasz Kacprzak 

6""" 

7 

8from warnings import warn 

9 

10import numpy as np 

11 

12 

13def apply_pycosmo_distfun(dist_fun, z): 

14 z_unique, inv_ind = np.unique(z, return_inverse=True) 

15 dist = dist_fun(a=1 / (1 + z_unique))[inv_ind] # Unit: Mpc 

16 return dist 

17 

18 

19def r50_phys_to_ang(r50_phys, cosmo, z): 

20 """ 

21 Convert physical size to angular size 

22 :param r50_phys: physical size of the galaxy in kpc 

23 :param cosmo: cosmological model 

24 :param z: redshift of the galaxy 

25 :return: angular size of the galaxy in arcsec 

26 """ 

27 d_a = apply_pycosmo_distfun(cosmo.background.dist_ang_a, z) # unit: Mpc 

28 r50_ang = r50_phys / (1e3 * d_a) # unit: rad 

29 r50_ang = np.rad2deg(r50_ang) * 60**2 # unit: arcsec 

30 return r50_ang 

31 

32 

33########################################################## 

34# 

35# Sampling galaxy size 

36# 

37########################################################## 

38 

39 

40def sample_r50_for_galaxy_type(z, abs_mag, cosmo, par, galaxy_type): 

41 """ 

42 Sample the physical size of a galaxy given its absolute magnitude 

43 and redshift. The physical size is sampled from a log-normal 

44 distribution with a mean that depends on the absolute magnitude 

45 and redshift. The physical size is then converted to an angular 

46 size using the cosmological distance. 

47 :param z: redshift of the galaxy 

48 :param abs_mag: absolute magnitude of the galaxy 

49 :param cosmo: cosmological model 

50 :param par: ucat parameters 

51 :param galaxy_type: type of the galaxy (e.g. "red", "blue") 

52 :return: angular size of the galaxy in pixels, in arcsec and physical size in kpc 

53 """ 

54 backwards_compatibility(par) 

55 

56 # set defaults 

57 if not hasattr(par, "logr50_sampling_method"): 

58 par.logr50_sampling_method = "single" 

59 

60 if par.logr50_sampling_method == "single": 

61 r50_phys = sample_r50_phys( 

62 abs_mag_shift=abs_mag - par.logr50_phys_M0, 

63 logr50_phys_std=par.logr50_phys_std, 

64 logr50_phys_mean_slope=par.logr50_phys_mean_slope, 

65 logr50_phys_mean_intcpt=par.logr50_phys_mean_intcpt, 

66 logr50_alpha=par.logr50_alpha, 

67 z=z, 

68 ) 

69 

70 elif par.logr50_sampling_method == "red_blue": 

71 r50_phys = sample_r50_phys( 

72 abs_mag_shift=abs_mag - par.logr50_phys_M0, 

73 logr50_phys_std=getattr(par, f"logr50_phys_std_{galaxy_type}"), 

74 logr50_phys_mean_slope=getattr( 

75 par, f"logr50_phys_mean_slope_{galaxy_type}" 

76 ), 

77 logr50_phys_mean_intcpt=getattr( 

78 par, f"logr50_phys_mean_intcpt_{galaxy_type}" 

79 ), 

80 logr50_alpha=getattr(par, f"logr50_alpha_{galaxy_type}"), 

81 z=z, 

82 ) 

83 

84 elif par.logr50_sampling_method == "sdss_fit": 

85 # this sampling is based on https://arxiv.org/abs/astro-ph/0301527 

86 sigma1 = getattr(par, f"logr50_sdss_fit_sigma1_{galaxy_type}") 

87 sigma2 = getattr(par, f"logr50_sdss_fit_sigma2_{galaxy_type}") 

88 M0 = getattr(par, f"logr50_sdss_fit_M0_{galaxy_type}") 

89 

90 # Equation (16) 

91 std = sigma2 + ((sigma1 - sigma2) / (1 + 10 ** (-0.8 * (abs_mag - M0)))) 

92 

93 if galaxy_type == "red": 

94 # Equation (14) 

95 slope = -0.4 * par.logr50_sdss_fit_a_red 

96 intcp = par.logr50_sdss_fit_b_red 

97 elif galaxy_type == "blue": 

98 alpha = par.logr50_sdss_fit_alpha_blue 

99 beta = par.logr50_sdss_fit_beta_blue 

100 gamma = par.logr50_sdss_fit_gamma_blue 

101 # Equation (15) 

102 # (intcp is not a classic intercept, but we can reuse sample_r50_phys) 

103 slope = -0.4 * alpha 

104 intcp = (beta - alpha) * np.log10(1 + 10 ** (-0.4 * (abs_mag - M0))) + gamma 

105 else: 

106 raise ValueError( 

107 f"unsupported galaxy_type={galaxy_type} for" 

108 " logr50_sampling_method=sdss_fit" 

109 ) 

110 

111 r50_phys = sample_r50_phys( 

112 abs_mag_shift=abs_mag, 

113 logr50_phys_std=std, 

114 logr50_phys_mean_slope=slope, 

115 logr50_phys_mean_intcpt=intcp, 

116 logr50_alpha=getattr(par, f"logr50_alpha_{galaxy_type}"), 

117 z=z, 

118 ) 

119 

120 else: 

121 raise Exception( 

122 f"unsupported logr50_sampling_method={par.logr50_sampling_method}" 

123 ) 

124 

125 # convert physical size to observerd angular size 

126 r50_ang = r50_phys_to_ang(r50_phys, cosmo, z) 

127 

128 r50_ang_pix = r50_ang.astype(par.catalog_precision) / par.pixscale 

129 r50_ang_arcsec = r50_ang.astype(par.catalog_precision) 

130 r50_phys = r50_phys.astype(par.catalog_precision) 

131 

132 return r50_ang_pix, r50_ang_arcsec, r50_phys 

133 

134 

135def sample_r50_phys( 

136 abs_mag_shift, 

137 logr50_phys_std, 

138 logr50_phys_mean_slope, 

139 logr50_phys_mean_intcpt, 

140 logr50_alpha=0.0, 

141 z=0.0, 

142): 

143 r50_phys = np.random.normal(loc=0, scale=logr50_phys_std, size=len(abs_mag_shift)) 

144 r50_phys += abs_mag_shift * logr50_phys_mean_slope + logr50_phys_mean_intcpt 

145 r50_phys = np.exp(r50_phys) # unit: kpc 

146 r50_phys *= (1 + z) ** logr50_alpha 

147 

148 return r50_phys 

149 

150 

151def backwards_compatibility(par): 

152 if hasattr(par, "sample_r50_model"): 

153 warn( 

154 DeprecationWarning( 

155 "sample_r50_model is deprecated, " 

156 "define logr50_phys_M0 or logr50_sdss_fit_M0 instead" 

157 ), 

158 stacklevel=1, 

159 ) 

160 if par.sample_r50_model == "base": 

161 par.logr50_phys_M0 = 0.0 

162 elif par.sample_r50_model == "shift20": 

163 par.logr50_phys_M0 = -20.0