diff --git a/.github/workflows/build-docs.yml b/.github/workflows/build-docs.yml index c9e56dc9..56a099e0 100644 --- a/.github/workflows/build-docs.yml +++ b/.github/workflows/build-docs.yml @@ -9,7 +9,10 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: ammaraskar/sphinx-action@master - with: - docs-folder: "doc" - pre-build-command: 'pip install -U sphinx_rtd_theme docutils' + - uses: actions/setup-python@v5 + - name: Install dependencies + run: pip install -r doc/requirements.txt + - name: Build documentation + env: + SPHINXOPTS: -W --keep-going + run: make html -C doc diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 71532212..6f79128e 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -10,3 +10,7 @@ submodules: sphinx: configuration: doc/conf.py fail_on_warning: true + +python: + install: + - requirements: doc/requirements.txt diff --git a/doc/hdf5_properties/galaxies.rst b/doc/hdf5_properties/galaxies.rst index f520a5a7..c1385cdb 100644 --- a/doc/hdf5_properties/galaxies.rst +++ b/doc/hdf5_properties/galaxies.rst @@ -180,4 +180,3 @@ * ``l_y``: total angular momentum component y of galaxy [Msun pMpc km/s]. From VELOCIraptor. * ``l_z``: total angular momentum component z of galaxy [Msun pMpc km/s]. From VELOCIraptor. * ``main_progenitor``: =1 if subhalo is the main progenitor' =0 otherwise. - diff --git a/doc/requirements.txt b/doc/requirements.txt new file mode 100644 index 00000000..4e09d9c3 --- /dev/null +++ b/doc/requirements.txt @@ -0,0 +1,2 @@ +Sphinx==7.3.7 +sphinx-rtd-theme==2.0.0 diff --git a/src/dark_matter_halos.cpp b/src/dark_matter_halos.cpp index 776b24d5..15119626 100644 --- a/src/dark_matter_halos.cpp +++ b/src/dark_matter_halos.cpp @@ -363,7 +363,7 @@ void DarkMatterHalos::cooling_gas_sAM(Subhalo &subhalo, double z){ float DarkMatterHalos::enclosed_total_mass(const Subhalo &subhalo, double z, float r){ - ConstGalaxyPtr galaxy; + ConstGalaxyPtr galaxy = nullptr; double mvir = 0; double rvir = 0; diff --git a/standard_plots/calculate_co_emission.py b/standard_plots/calculate_co_emission.py index 43dbc3db..136f90ea 100644 --- a/standard_plots/calculate_co_emission.py +++ b/standard_plots/calculate_co_emission.py @@ -40,7 +40,18 @@ sigma_gas = 20.0 #km/s for CO thresh_thin_disk = 0.01 -thresh_super_edd = 1.0 +thresh_super_edd = 4.0 + + +delta_adaf = 0.2 +mdotcrit_adaf = 0.01 + +alpha_adaf = 0.1 +beta = 1 - alpha_adaf / 0.55 +low_accretion_adaf = 0.001 * (delta_adaf / 0.0005) * (1 - beta) / beta * np.power(alpha_adaf, 2.0) +constant_lowlum_adaf = (delta_adaf / 0.0005) * (1 - beta) / 0.5 * 6 +constant_highlum_adaf = beta / 0.5 / np.power(alpha_adaf, 2) * 6 + gammaSFR = 1.0 alphaFx = 1.0 @@ -48,7 +59,70 @@ zsun = 0.0189 -def prepare_data(hdf5_data, index, model_dir, snapshot, subvol, obsdir): +def plot_co_sled(plt, LCO, snap, outdir): + + xj = np.array([1,2,3,4,5,6,7,8,9,10]) + colors = ['Orange', 'SeaGreen', 'Purple'] + #plot evolution of SMGs in the UVJ plane + xtit="$\\rm J_{\\rm upper}$" + ytit="$\\rm S_{\\rm CO}/S_{\\rm CO(1-0)}$" + + xmin, xmax, ymin, ymax = 0.5, 10, 0, 20 + xleg = xmax - 0.18 * (xmax-xmin) + yleg = ymin + 0.1 * (ymax-ymin) + + fig = plt.figure(figsize=(5,4)) + ax = fig.add_subplot(111) + common.prepare_ax(ax, xmin, xmax, ymin, ymax, xtit, ytit, locators=(1, 1, 1, 1)) + + med_co = np.zeros(shape = (3,10)) + med_co[0,0] = 1 + med_co[1,0] = 1 + med_co[2,0] = 1 + + for i in range(1,10): + ind=np.where(LCO[:,0] > 0) + med_co[0,i] = np.median(LCO[ind,i]/LCO[ind,0]) + per = np.percentile(LCO[ind,i]/LCO[ind,0], [16, 50, 84]) + med_co[0,i] = per[1] + med_co[1,i] = per[0] + med_co[2,i] = per[2] + + ax.fill_between(xj, med_co[1,:], med_co[2,:], color='orange') + ax.plot(xj, med_co[0,:], linestyle='solid', color='darkred') + + namefig = "example_co_sleds_massive_sf_" + snap + ".pdf" + common.savefig(outdir, fig, namefig) + + +def radiation_efficiency(spin, efficiency): + + eff = efficiency + a = abs(spin); + a2 = a**2 + + z1 = 1 + np.power(1 - a2, 0.333) * (np.power(1 + a, 0.333) + np.power(1 - a, 0.333)) + z2 = np.sqrt(3 * a2 + np.power(z1, 2)) + + r_lso = np.zeros(shape = len(spin)) + + ind = np.where(spin >= 0) + r_lso[ind] = 3 + z2[ind] - np.sqrt((3 - z1[ind] ) * (3 + z1[ind] + 2 * z2[ind] )) + ind = np.where(spin < 0) + r_lso[ind] = 3 + z2[ind] + np.sqrt((3 - z1[ind] ) * (3 + z1[ind] + 2 * z2[ind] )) + + eff[0,:] = 1 - np.sqrt(1 - 2 / (3 * r_lso)) + + ind = np.where(eff[0,:] < 0) + eff[0,ind] = 0.07 + ind = np.where(eff[0,:] > 0.5) + eff[0,ind] = 0.5 + eff[1,:] = r_lso + + efficiency = eff + + +def prepare_data(hdf5_data, index, model_dir, snapshot, subvol, obsdir, read_spin, test_co_sleds): #read ascii table from Bayet et al. (2011) datanu = np.genfromtxt(os.path.join(obsdir, 'Models', 'CO','emlines_carbonmonoxide_restnu.data'), delimiter=' ', skip_header=11) @@ -83,8 +157,19 @@ def prepare_data(hdf5_data, index, model_dir, snapshot, subvol, obsdir): MaxCRs = np.max(CRs) #define maximum CRs probed by the models. # read galaxy information in lightcone - (h0, _, idgal, mmol_b, mmol_d, rd, rb, mzd, mzb, sfr_d, sfr_b, mgd, mgb, - mbh, mbh_acc_hh, mbh_acc_sb, mdisk, mbulge) = hdf5_data + if(read_spin): + (h0, _, idgal, mmol_b, mmol_d, rd, rb, mzd, mzb, sfr_d, sfr_b, mgd, mgb, + mbh, mbh_acc_hh, mbh_acc_sb, mdisk, mbulge, spin) = hdf5_data + else: + (h0, _, idgal, mmol_b, mmol_d, rd, rb, mzd, mzb, sfr_d, sfr_b, mgd, mgb, + mbh, mbh_acc_hh, mbh_acc_sb, mdisk, mbulge) = hdf5_data + + efficiency = np.zeros(shape=(2,len(mbh))) + if(read_spin): + radiation_efficiency(spin, efficiency) + else: + efficiency[0,:] = 0.1 + efficiency[1,:] = 6 #define metallicities for disk and bulge zd = np.zeros(shape = len(mzd)) @@ -109,7 +194,7 @@ def prepare_data(hdf5_data, index, model_dir, snapshot, subvol, obsdir): # define bolometric luminosities using Eqs 1 in Amarantidis et al. (2019) ind = np.where(mbh_acc_hh+mbh_acc_sb > 0) - Lthin_disk[ind] = 0.1*pow(c_light*1e2, 2.0) *(mbh_acc_hh[ind]+mbh_acc_sb[ind])/h0 * 6.329113924050633e-24 #in 1e40 ergs/s + Lthin_disk[ind] = efficiency[0,ind] * pow(c_light*1e2, 2.0) *(mbh_acc_hh[ind]+mbh_acc_sb[ind])/h0 * 6.329113924050633e-24 #in 1e40 ergs/s # thin disks ind = np.where((mnorm > thresh_thin_disk) & (mnorm < thresh_super_edd)) @@ -117,10 +202,12 @@ def prepare_data(hdf5_data, index, model_dir, snapshot, subvol, obsdir): # super-eddington ind = np.where(mnorm > thresh_super_edd) Lbol[ind] = thresh_super_edd * (1.0 + np.log(mnorm[ind]/thresh_super_edd)) * Ledd[ind] #in 1e40 ergs/s - # ADAfs adopting alpa_ADAF=0.1 and beta = 0.81 and r_lso = 6 (spin =1) - ind = np.where((mnorm < thresh_thin_disk) & (mnorm > 0)) - Lbol[ind] = 0.2 * Lthin_disk[ind] * (mnorm[ind]/1e-2) * 1.63 - + # ADAfs + ind = np.where((mnorm < thresh_thin_disk) & (mnorm > 0) & (mnorm > low_accretion_adaf)) + Lbol[ind] = 0.2 * efficiency[0,ind] * pow(c_light*1e2, 2.0) *(mbh_acc_hh[ind]+mbh_acc_sb[ind])/h0 * 6.329113924050633e-24 * mnorm[ind] * constant_highlum_adaf / efficiency[1,ind] + ind = np.where((mnorm < thresh_thin_disk) & (mnorm > 0) & (mnorm < low_accretion_adaf)) + Lbol[ind] = 0.0002 * efficiency[0,ind] * pow(c_light*1e2, 2.0) *(mbh_acc_hh[ind]+mbh_acc_sb[ind])/h0 * 6.329113924050633e-24 * constant_lowlum_adaf / efficiency[1,ind] + # L-hard-xrays Eq 34 from Griffin et al. (2018) Lrat = np.log10(Lbol/Lsun) Lx = pow(10.0, -1.54 - 0.24*Lrat - 0.012 * pow(Lrat, 2.0) + 0.0015 * pow(Lrat, 3.0)) * Lbol #in 1e40 erg/s @@ -212,21 +299,47 @@ def get_co_emissions(mmol, zcoldg, guv, fx): hf.create_dataset('frequency_co_rest', data=nuco) hf.close() + if(test_co_sleds): + ssfr_cut = 10**(-1 + 0.5 * redshift) + ind = np.where( ((mdisk + mbulge)/h0 > 1e10) & ((sfr_d + sfr_b)/(mdisk + mbulge) > ssfr_cut)) + return LCOd[ind] + LCOb[ind] + def main(model_dir, output_dir, redshift_table, subvols, obs_dir): - snapshots = range(30,200) + plt = common.load_matplotlib() + test_co_sleds = False + read_spin = True + + if(test_co_sleds): + zlist = [0, 0.5, 1, 1.5, 2, 2.5, 3.0] + snapshots = redshift_table[zlist] + else: + snapshots = range(61,200) + # Loop over redshift and subvolumes plt = common.load_matplotlib() - fields = {'galaxies': ('id_galaxy', 'mmol_bulge', 'mmol_disk', 'rgas_disk', - 'rgas_bulge', 'mgas_metals_disk', 'mgas_metals_bulge', 'sfr_disk', 'sfr_burst', - 'mgas_disk', 'mgas_bulge','m_bh','bh_accretion_rate_hh','bh_accretion_rate_sb', - 'mstars_disk', 'mstars_bulge')} + + + if(read_spin): + fields = {'galaxies': ('id_galaxy', 'mmol_bulge', 'mmol_disk', 'rgas_disk', + 'rgas_bulge', 'mgas_metals_disk', 'mgas_metals_bulge', 'sfr_disk', 'sfr_burst', + 'mgas_disk', 'mgas_bulge','m_bh','bh_accretion_rate_hh','bh_accretion_rate_sb', + 'mstars_disk', 'mstars_bulge', 'bh_spin')} + else: + fields = {'galaxies': ('id_galaxy', 'mmol_bulge', 'mmol_disk', 'rgas_disk', + 'rgas_bulge', 'mgas_metals_disk', 'mgas_metals_bulge', 'sfr_disk', 'sfr_burst', + 'mgas_disk', 'mgas_bulge','m_bh','bh_accretion_rate_hh','bh_accretion_rate_sb', + 'mstars_disk', 'mstars_bulge')} for index,snapshot in enumerate(snapshots): for subv in subvols: hdf5_data = common.read_data(model_dir, snapshot, fields, [subv]) - prepare_data(hdf5_data, index, model_dir, snapshot, subv, obs_dir) + if(test_co_sleds): + LCO = prepare_data(hdf5_data, index, model_dir, snapshot, subv, obs_dir, read_spin, test_co_sleds) + plot_co_sled(plt, LCO, str(zlist[index]), output_dir) + else: + prepare_data(hdf5_data, index, model_dir, snapshot, subv, obs_dir, read_spin, test_co_sleds) if __name__ == '__main__': main(*common.parse_args())