Coverage for src / galsbi / ucat_sps / galaxy_population_models / galaxy_gas_ionization.py: 100%

1# Copyright (C) 2025 LMU Munich 

2# Author: Luca Tortorelli 

3# created: Apr 2025 

4 

5import numpy as np 

6from scipy.stats import truncnorm 

7 

8from galsbi.ucat_sps.config.common_sps import SOLAR_GAS_METALLICITY, SOLAR_LOG_OH_PLUS12 

9 

10 

11def sample_gas_ionization_Kashino19( 

12 Zgas_final, star_formation_rates, log10_stellar_masses, par, seed=None 

13): 

14 """ 

15 This function computes the gas ionization parameter using the relation in 

16 equation 12 of Kashino+19. We keep the electron density fixed to the value with 

17 which the CLOUDY templates of FSPS have been generated, i.e. ne~100 cm^-3. We sample 

18 logU from a Normal distribution with mean given by the Kashino+19 relation and 

19 scatter of 0.1 dex, as highlighted in Kashino+19. The Normal distribution is 

20 truncated within the range [-4, -1]. 

21 

22 :param Zgas_final: (array_like[n_gal,]) final gas-phase metallicity in absolute 

23 units. 

24 :param star_formation_rates: (array_like[n_gal,]) SFR averaged over the last 100 Myr 

25 of the galaxy SFH. Units are Msun/yr. 

26 :param log10_stellar_masses: (array_like[n_gal,]) log10 of the stellar masses 

27 sampled from the stellar mass function. 

28 :param par: (obj) par objects containing the Ucat parameters. 

29 :param seed: (int) random seed for the truncnorm sampling. 

30 :return logU: (array_like[n_gal,]) gas ionization parameters. 

31 """ 

32 if seed is None: 

33 seed = 42 

34 rng = np.random.default_rng(seed=seed) 

35 

36 log_sSFR = np.log10(star_formation_rates / (10**log10_stellar_masses)) 

37 logOH = np.log10(Zgas_final / SOLAR_GAS_METALLICITY) + SOLAR_LOG_OH_PLUS12 - 12 

38 log_ne = 2 # 100 cm^-3 

39 logU_mean = ( 

40 par.gas_ionization_upsilon0 

41 + par.gas_ionization_upsilon1 * (logOH + 4) 

42 + par.gas_ionization_upsilon2 * log_ne 

43 + par.gas_ionization_upsilon3 * (log_sSFR + 9) 

44 ) 

45 logUmin = (par.logUmin - logU_mean) / par.gas_ionization_scatter 

46 logUmax = (par.logUmax - logU_mean) / par.gas_ionization_scatter 

47 logU = truncnorm.rvs( 

48 logUmin, 

49 logUmax, 

50 loc=logU_mean, 

51 scale=par.gas_ionization_scatter, 

52 size=len(Zgas_final), 

53 random_state=rng, 

54 ) 

55 

56 return logU