diff --git a/standard_plots/angularmomentum.py b/standard_plots/angularmomentum.py index d726432..af976e7 100644 --- a/standard_plots/angularmomentum.py +++ b/standard_plots/angularmomentum.py @@ -34,6 +34,14 @@ MpcToKpc = 1e3 G = 4.299e-9 #Gravity constant in units of (km/s)^2 * Mpc/Msun XH = 1.33 +MSOLAR = 1.9891e30 +MSOLAR_g = MSOLAR * 1e3 +MPC2M = 3.0856775807e22 +MPC2CM = MPC2M * 1e2 +k_Boltzmann = 1.3806503e-23 +k_Boltzmann_erg = k_Boltzmann * 1.0e7 +Pressure_Conv = np.pi/ 2 * G * MSOLAR_g / (MPC2CM)**3 * 1e10 / k_Boltzmann_erg + #stellar mass bins mlow = 6.35 @@ -66,6 +74,164 @@ xlf = lbins + dl/2.0 +def v2halo (x, m, c, r): + + cpx = 1 + c * x + cx = c * x + cp = 1 + c + + v = G * m / r * (np.log(cpx) - cx/cpx) / (x * (np.log(cp) - c / cp)) + + return v + + +def v2disk (x, m, c, r): + + cx = c * x + + nom = c + 4.8 * c * np.exp(-0.35 * cx - 3.5 /cx) + denom = cx + cx**(-2) + 2.0 * cx**(-0.5) + + v = G * m / r * nom/denom + + return v + +def v2bulge (x, m, c, r): + + if(m <= 0): + return 0 + else: + + cx = c * x; + + nom = cx**(2.0) * c; + denom = (1 + cx**(2.0))**(1.5); + + v = G * m / r * nom/denom; + + return v + + +def compute_v_r(r, msdisk, rsdisk, mgdisk, rgdisk, mbulge, rbulge, mdm, rvir, c): + + #disk properties. Use composite size of disk. + cd_s = 0 + if(rsdisk > 0): + cd_s = rvir / (rsdisk / 1.67) + cd_g = 0 + if(rgdisk > 0): + cd_g = rvir / (rgdisk / 1.67) + + cb = 0 + if(rbulge > 0): + cb = rvir/(rbulge) + + #halo properties. + ch = c + mvir = mdm + + x = r / rvir; + + #Rotational velocity at the half-mass radius of the disk. + v2tot_d = v2halo(x, mvir, ch, rvir) + v2disk(x, msdisk, cd_s, rvir) + v2disk(x, mgdisk, cd_g, rvir) + v2bulge(x, mbulge, cb, rvir) + return np.sqrt(v2tot_d) + + +def midplane_pressure(Sigma_gas, Sigma_stars, r): + + #/** + # * This function calculate the midplane pressure of the disk, and returns it in units of K/cm^-3. + # */ + + hstar = 0.14 * r #scaleheight of the stellar disk; from Kregel et al. (2002). + veldisp_star = np.sqrt(np.pi * G * hstar * Sigma_stars) #stellar velocity dispersion in km/s. + + star_comp = 0 + + if ((Sigma_stars > 0) & (veldisp_star > 0)): + star_comp = (10.0 / veldisp_star) * Sigma_stars + + pressure = Pressure_Conv * Sigma_gas * (Sigma_gas + star_comp) #in units of K/cm^3. + + return pressure + +def fmol(Sigma_gas, Sigma_stars, r): + + Po = 34673.0 + beta_press = 0.92 + rmol = (midplane_pressure(Sigma_gas,Sigma_stars,r)/Po)**beta_press; + fmol = rmol/(1+rmol); + + if(fmol > 1): + return 1 + elif((fmol > 0) & (fmol < 1)): + return fmol + else: + return 0 + + + +def molecular_surface_density(r, mgas, rgas, mstar, rstar): + + re = rgas / 1.67 + sigma_gas0 = mgas / (2 * np.pi * re**2) + Sigma_gas = sigma_gas0 * np.exp(-r / re) + + if(Sigma_gas <= 0): + Sigma_gas = 0 + + sigma_HI_crit = 0.0 + # Check for low surface densities. + sigma_HI_crit = 0.5 * 1e12 #in Msun/Mpc^2 + if(Sigma_gas <= sigma_HI_crit): + return 0,0 + + Sigma_stars = 0 + # Define Sigma_stars only if stellar mass and radius are positive. + if((rstar > 0) & (mstar > 0)): + rse = rstar / 1.67 + sigma_star0 = mstar / (2 * np.pi * rse**2) + Sigma_stars = sigma_star0 * np.exp(-r / rse) + + fmol_in = fmol(Sigma_gas, Sigma_stars, r) + sigmas_mol = 2 * np.pi * fmol_in * Sigma_gas * r #Add the 2PI*r to Sigma_gas to make integration. + sigmas_atom = 2 * np.pi * (1 - fmol_in) * Sigma_gas * r #Add the 2PI*r to Sigma_gas to make integration. + + return sigmas_mol, sigmas_atom + +def compute_numerical_jmol(msdisk, rsdisk, mgdisk, rgdisk, mbulge, rbulge, mdm, rvir, c, mmol): + + #integrate from r=0 to 10*r50_disk + ngals = len(msdisk) + nrbins = 100 + jmol_comp = np.zeros(shape = ngals) + jatom_comp = np.zeros(shape = ngals) + + mmol_comp = np.zeros(shape = ngals) + + #Msun/pc^2 Msun/(1e6 pc)^2, + for i in range(0,ngals): + if((mgdisk[i] > 0) & (rgdisk[i] > 0) & (mmol[i] > 0)): + sum_int = 0 + sum_int_hi = 0 + sum_int_mol = 0 + sum_int_atomass = 0 + dr = 5.0 * rgdisk[i] / nrbins + for j in range(1,nrbins): + r = dr * j + (sigma_mol_bin, sigma_atom_bin) = molecular_surface_density(r, mgdisk[i], rgdisk[i], msdisk[i], rsdisk[i]) + v_at_r = compute_v_r(r, msdisk[i], rsdisk[i], mgdisk[i], rgdisk[i], mbulge[i], rbulge[i], mdm[i], rvir[i], c[i]) + sum_int += sigma_mol_bin * v_at_r * r * dr + sum_int_hi += sigma_atom_bin * v_at_r * r * dr + sum_int_mol += sigma_mol_bin * dr + sum_int_atomass += sigma_atom_bin * dr + mmol_comp[i] = sum_int_mol + jmol_comp[i] = sum_int / mmol_comp[i] + jatom_comp[i] = sum_int_hi / sum_int_atomass + + return jmol_comp, jatom_comp + + def prepare_data(hdf5_data, index, sam_stars_disk, sam_gas_disk_atom, sam_gas_disk_atom2, sam_gas_disk_mol, sam_gas_disk_atom_ms, sam_halo, sam_ratio_halo_disk, sam_ratio_halo_gal, sam_ratio_halo_disk_gas, disk_size_sat, disk_size_cen, bulge_size, sam_vs_sam_halo_disk, sam_vs_sam_halo_gal, sam_vs_sam_halo_disk_gas, sam_bar, sam_stars, sam_stars2, vmax_halo_gal, sam_barv2, disk_size, rcomb): #, phot_data, phot_data_nod, nbands): @@ -80,8 +246,24 @@ def prepare_data(hdf5_data, index, sam_stars_disk, sam_gas_disk_atom, sam_gas_di mBH, rdisk, rbulge, rg_disk, rg_bulge, typeg, specific_angular_momentum_disk_star, specific_angular_momentum_bulge_star, specific_angular_momentum_disk_gas, specific_angular_momentum_bulge_gas, specific_angular_momentum_disk_gas_atom, specific_angular_momentum_disk_gas_mol, lambda_sub, mvir_s, mvir, matom_disk, mmol_disk, mgas_disk, - matom_bulge, mmol_bulge, mgas_bulge, sfr_disk, sfr_bulge, vmax, lx, ly, lz) = hdf5_data + matom_bulge, mmol_bulge, mgas_bulge, sfr_disk, sfr_bulge, vmax, lx, ly, lz, cnfw, vvir_sub) = hdf5_data + + rvir_sb = G * mvir_s / vvir_sub**2 + jmol = (specific_angular_momentum_disk_gas_mol * mmol_disk + mmol_bulge * specific_angular_momentum_bulge_gas) / (mmol_disk + mmol_bulge) + jatom = (specific_angular_momentum_disk_gas_atom * matom_disk + matom_bulge * specific_angular_momentum_bulge_gas) / (matom_disk + matom_bulge) + + massive_gals = np.where(((mdisk + mbulge)/h0 >= 1e8) & (mmol_disk > 0)) + jmol_approx = (specific_angular_momentum_disk_gas_mol[massive_gals] * mmol_disk[massive_gals] + mmol_bulge[massive_gals] * specific_angular_momentum_bulge_gas[massive_gals]) / (mmol_disk[massive_gals] + mmol_bulge[massive_gals]) + jmol_comp_detail, jatom_comp_detail = compute_numerical_jmol(mdisk[massive_gals]/h0, rdisk[massive_gals]/h0, mgas_disk[massive_gals], rg_disk[massive_gals]/h0, mbulge[massive_gals]/h0 + mgas_bulge[massive_gals]/h0, rbulge[massive_gals]/h0, mvir_s[massive_gals]/h0, rvir_sb[massive_gals]/h0, cnfw[massive_gals], mmol_disk[massive_gals]/h0) + jmol_detailed = (jmol_comp_detail * mmol_disk[massive_gals] + mmol_bulge[massive_gals] * specific_angular_momentum_bulge_gas[massive_gals]) / (mmol_disk[massive_gals] + mmol_bulge[massive_gals]) + jatom_detailed = (jatom_comp_detail * matom_disk[massive_gals] + matom_bulge[massive_gals] * specific_angular_momentum_bulge_gas[massive_gals]) / (matom_disk[massive_gals] + matom_bulge[massive_gals]) + jmol[massive_gals] = jmol_detailed + jatom[massive_gals] = jatom_detailed + #for a,b,c,d,e,f in zip(jmol_approx, jmol_comp_detail, (mmol_disk[massive_gals] + mmol_bulge[massive_gals])/h0, (mdisk[massive_gals] + mbulge[massive_gals])/h0, mbulge[massive_gals] / (mdisk[massive_gals] + mbulge[massive_gals]), (sfr_disk[massive_gals] + sfr_bulge[massive_gals])/h0/1e9): + # print(a*1e3,b*1e3,c,d,e,f) + #print("detailed jmol calculation", jmol_comp_detail) + ind = np.where(mmol_disk > 0) print("sAM mol when mgas_mol>0", specific_angular_momentum_disk_gas_mol[ind]) print("number of galaxies with mmol_disk > 0", len(mmol_disk[ind])) @@ -131,11 +313,11 @@ def prepare_data(hdf5_data, index, sam_stars_disk, sam_gas_disk_atom, sam_gas_di jbar = (specific_angular_momentum_disk_star * mdisk + specific_angular_momentum_disk_gas * mgas_disk ) / (mdisk + mbulge + mgas_disk + mgas_bulge) #specific_angular_momentum_bulge_star * mbulge + specific_angular_momentum_bulge_gas * mgas_bulge) / (mdisk + mbulge + mgas_disk + mgas_bulge) mbar = (mdisk + mbulge + mgas_disk + mgas_bulge) - jbarv2 = (specific_angular_momentum_disk_star * mdisk + specific_angular_momentum_disk_gas * mgas_disk + specific_angular_momentum_bulge_star * mbulge + specific_angular_momentum_bulge_gas * mgas_bulge) / (mdisk + mbulge + mgas_disk + mgas_bulge) + jbarv2 = (specific_angular_momentum_disk_star * mdisk + jmol * (mmol_disk + mmol_bulge) + jatom * (matom_disk + matom_bulge)) / (mdisk + mbulge + mmol_disk + mmol_bulge + matom_disk + matom_bulge) #(specific_angular_momentum_disk_gas_atom * matom_disk + specific_angular_momentum_disk_star * mdisk + specific_angular_momentum_bulge_star * mbulge + matom_bulge * specific_angular_momentum_bulge_gas) / (mdisk + mbulge + matom_disk + matom_bulge) - mbarv2 = (mdisk + mbulge + mgas_disk + mgas_bulge) - jatom = (specific_angular_momentum_disk_gas_atom * matom_disk + matom_bulge * specific_angular_momentum_bulge_gas) / (matom_disk + matom_bulge) - jmol = (specific_angular_momentum_disk_gas_mol * mmol_disk + mmol_bulge * specific_angular_momentum_bulge_gas) / (mmol_disk + mmol_bulge) + mbarv2 = (mdisk + mbulge + mmol_disk + mmol_bulge + matom_disk + matom_bulge) + + jgas = (specific_angular_momentum_disk_gas * mgas_disk + specific_angular_momentum_bulge_gas * mgas_bulge) / (mgas_bulge + mgas_disk) mgas = (mgas_bulge + mgas_disk) @@ -1452,7 +1634,8 @@ def main(modeldir, outdir, redshift_table, subvols, obsdir): 'specific_angular_momentum_disk_gas', 'specific_angular_momentum_bulge_gas', 'specific_angular_momentum_disk_gas_atom', 'specific_angular_momentum_disk_gas_mol', 'lambda_subhalo', 'mvir_subhalo', 'mvir_hosthalo', 'matom_disk', 'mmol_disk', 'mgas_disk', - 'matom_bulge', 'mmol_bulge', 'mgas_bulge','sfr_disk', 'sfr_burst','vmax_subhalo', 'l_x', 'l_y', 'l_z')} + 'matom_bulge', 'mmol_bulge', 'mgas_bulge','sfr_disk', 'sfr_burst','vmax_subhalo', 'l_x', 'l_y', 'l_z', + 'cnfw_subhalo', 'vvir_subhalo')} fields_sed = {'SED/ab_dust': ('bulge_d','bulge_m','bulge_t','disk','total'),} fields_sed_nod = {'SED/ab_nodust': ('bulge_d','bulge_m','bulge_t','disk','total')}