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, pixscale): 

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

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

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

23 r50_ang /= pixscale # unit: pixel 

24 return r50_ang 

25 

26 

27########################################################## 

28# 

29# Sampling galaxy size 

30# 

31########################################################## 

32 

33 

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

35 backwards_compatibility(par) 

36 

37 # set defaults 

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

39 par.logr50_sampling_method = "single" 

40 

41 if par.logr50_sampling_method == "single": 

42 r50_phys = sample_r50_phys( 

43 abs_mag_shift=abs_mag - par.logr50_phys_M0, 

44 logr50_phys_std=par.logr50_phys_std, 

45 logr50_phys_mean_slope=par.logr50_phys_mean_slope, 

46 logr50_phys_mean_intcpt=par.logr50_phys_mean_intcpt, 

47 logr50_alpha=par.logr50_alpha, 

48 z=z, 

49 ) 

50 

51 elif par.logr50_sampling_method == "red_blue": 

52 r50_phys = sample_r50_phys( 

53 abs_mag_shift=abs_mag - par.logr50_phys_M0, 

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

55 logr50_phys_mean_slope=getattr( 

56 par, f"logr50_phys_mean_slope_{galaxy_type}" 

57 ), 

58 logr50_phys_mean_intcpt=getattr( 

59 par, f"logr50_phys_mean_intcpt_{galaxy_type}" 

60 ), 

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

62 z=z, 

63 ) 

64 

65 elif par.logr50_sampling_method == "sdss_fit": 

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

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

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

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

70 

71 # Equation (16) 

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

73 

74 if galaxy_type == "red": 

75 # Equation (14) 

76 slope = -0.4 * par.logr50_sdss_fit_a_red 

77 intcp = par.logr50_sdss_fit_b_red 

78 elif galaxy_type == "blue": 

79 alpha = par.logr50_sdss_fit_alpha_blue 

80 beta = par.logr50_sdss_fit_beta_blue 

81 gamma = par.logr50_sdss_fit_gamma_blue 

82 # Equation (15) 

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

84 slope = -0.4 * alpha 

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

86 else: 

87 raise ValueError( 

88 f"unsupported galaxy_type={galaxy_type} for" 

89 " logr50_sampling_method=sdss_fit" 

90 ) 

91 

92 r50_phys = sample_r50_phys( 

93 abs_mag_shift=abs_mag, 

94 logr50_phys_std=std, 

95 logr50_phys_mean_slope=slope, 

96 logr50_phys_mean_intcpt=intcp, 

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

98 z=z, 

99 ) 

100 

101 else: 

102 raise Exception( 

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

104 ) 

105 

106 # convert physical size to observerd angular size 

107 r50_ang = r50_phys_to_ang(r50_phys, cosmo, z, par.pixscale) 

108 

109 return r50_ang.astype(par.catalog_precision) 

110 

111 

112def sample_r50_phys( 

113 abs_mag_shift, 

114 logr50_phys_std, 

115 logr50_phys_mean_slope, 

116 logr50_phys_mean_intcpt, 

117 logr50_alpha=0.0, 

118 z=0.0, 

119): 

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

121 r50_phys += abs_mag_shift * logr50_phys_mean_slope + logr50_phys_mean_intcpt 

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

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

124 

125 return r50_phys 

126 

127 

128def backwards_compatibility(par): 

129 if hasattr(par, "sample_r50_model"): 

130 warn( 

131 DeprecationWarning( 

132 "sample_r50_model is deprecated, " 

133 "define logr50_phys_M0 or logr50_sdss_fit_M0 instead" 

134 ), 

135 stacklevel=1, 

136 ) 

137 if par.sample_r50_model == "base": 

138 par.logr50_phys_M0 = 0.0 

139 elif par.sample_r50_model == "shift20": 

140 par.logr50_phys_M0 = -20.0