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

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 

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

24 :param cosmo: cosmological model 

25 :param z: redshift of the galaxy 

26 :return: angular size of the galaxy in arcsec 

27 """ 

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

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

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

31 return r50_ang 

32 

33 

34########################################################## 

35# 

36# Sampling galaxy size 

37# 

38########################################################## 

39 

40 

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

42 """ 

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

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

45 distribution with a mean that depends on the absolute magnitude 

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

47 size using the cosmological distance. 

48 

49 :param z: redshift of the galaxy 

50 :param abs_mag: absolute magnitude of the galaxy 

51 :param cosmo: cosmological model 

52 :param par: ucat parameters 

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

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

55 """ 

56 backwards_compatibility(par) 

57 

58 # set defaults 

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

60 par.logr50_sampling_method = "single" 

61 

62 if par.logr50_sampling_method == "single": 

63 r50_phys = sample_r50_phys( 

64 abs_mag_shift=abs_mag - par.logr50_phys_M0, 

65 logr50_phys_std=par.logr50_phys_std, 

66 logr50_phys_mean_slope=par.logr50_phys_mean_slope, 

67 logr50_phys_mean_intcpt=par.logr50_phys_mean_intcpt, 

68 logr50_alpha=par.logr50_alpha, 

69 z=z, 

70 ) 

71 

72 elif par.logr50_sampling_method == "red_blue": 

73 r50_phys = sample_r50_phys( 

74 abs_mag_shift=abs_mag - par.logr50_phys_M0, 

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

76 logr50_phys_mean_slope=getattr( 

77 par, f"logr50_phys_mean_slope_{galaxy_type}" 

78 ), 

79 logr50_phys_mean_intcpt=getattr( 

80 par, f"logr50_phys_mean_intcpt_{galaxy_type}" 

81 ), 

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

83 z=z, 

84 ) 

85 

86 elif par.logr50_sampling_method == "sdss_fit": 

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

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

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

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

91 

92 # Equation (16) 

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

94 

95 if galaxy_type == "red": 

96 # Equation (14) 

97 slope = -0.4 * par.logr50_sdss_fit_a_red 

98 intcp = par.logr50_sdss_fit_b_red 

99 elif galaxy_type == "blue": 

100 alpha = par.logr50_sdss_fit_alpha_blue 

101 beta = par.logr50_sdss_fit_beta_blue 

102 gamma = par.logr50_sdss_fit_gamma_blue 

103 # Equation (15) 

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

105 slope = -0.4 * alpha 

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

107 else: 

108 raise ValueError( 

109 f"unsupported galaxy_type={galaxy_type} for" 

110 " logr50_sampling_method=sdss_fit" 

111 ) 

112 

113 r50_phys = sample_r50_phys( 

114 abs_mag_shift=abs_mag, 

115 logr50_phys_std=std, 

116 logr50_phys_mean_slope=slope, 

117 logr50_phys_mean_intcpt=intcp, 

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

119 z=z, 

120 ) 

121 

122 else: 

123 raise Exception( 

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

125 ) 

126 

127 # convert physical size to observerd angular size 

128 r50_ang = r50_phys_to_ang(r50_phys, cosmo, z) 

129 

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

131 r50_ang_arcsec = r50_ang.astype(par.catalog_precision) 

132 r50_phys = r50_phys.astype(par.catalog_precision) 

133 

134 return r50_ang_pix, r50_ang_arcsec, r50_phys 

135 

136 

137def sample_r50_phys( 

138 abs_mag_shift, 

139 logr50_phys_std, 

140 logr50_phys_mean_slope, 

141 logr50_phys_mean_intcpt, 

142 logr50_alpha=0.0, 

143 z=0.0, 

144): 

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

146 r50_phys += abs_mag_shift * logr50_phys_mean_slope + logr50_phys_mean_intcpt 

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

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

149 

150 return r50_phys 

151 

152 

153def backwards_compatibility(par): 

154 if hasattr(par, "sample_r50_model"): 

155 warn( 

156 DeprecationWarning( 

157 "sample_r50_model is deprecated, " 

158 "define logr50_phys_M0 or logr50_sdss_fit_M0 instead" 

159 ), 

160 stacklevel=1, 

161 ) 

162 if par.sample_r50_model == "base": 

163 par.logr50_phys_M0 = 0.0 

164 elif par.sample_r50_model == "shift20": 

165 par.logr50_phys_M0 = -20.0