diff --git a/colibre/auto_plotter/snia_rates.yml b/colibre/auto_plotter/snia_rates.yml index 9d968293..4d05a617 100644 --- a/colibre/auto_plotter/snia_rates.yml +++ b/colibre/auto_plotter/snia_rates.yml @@ -1,4 +1,4 @@ -stellar_mass_snia_rates_50: +stellar_mass_snia_rates_per_stellar_mass_50: type: "scatter" legend_loc: "upper left" x: @@ -7,10 +7,10 @@ stellar_mass_snia_rates_50: start: 1e6 end: 1e12 y: - quantity: "snia_rates.snia_rates_50_kpc" - units: "1/year" - start: 3e-8 - end: 1e-1 + quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" + units: "1/year/Solar_Mass" + start: 3e-15 + end: 3e-12 mean: plot: true log: true @@ -23,17 +23,19 @@ stellar_mass_snia_rates_50: value: 1e12 units: Solar_Mass metadata: - title: Stellar Mass-SNIa rate relation (50 kpc aperture) + title: Stellar Mass-SNIa rate per stellar mass relation (50 kpc aperture) caption: Uses a 50 kpc 3D aperture. - section: SNIa Rate + section: SNIa Rate (Stellar Mass) observational_data_bracket_width: 10.0 observational_data: - - filename: GalaxyStellarMassSNIaRate/Kistler2014.hdf5 - - filename: GalaxyStellarMassSNIaRate/Graur2015.hdf5 - - filename: GalaxyStellarMassSNIaRate/Graur2017.hdf5 - - filename: GalaxyStellarMassSNIaRate/Wiseman2021.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Kistler2014.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2015.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2017.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Brown2019_vol_limited.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Brown2019_full.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Wiseman2021.hdf5 -stellar_mass_snia_rates_active_only_50: +stellar_mass_snia_rates_per_stellar_mass_active_only_50: type: "scatter" comment: "Active only" comment_loc: "lower left" @@ -45,10 +47,10 @@ stellar_mass_snia_rates_active_only_50: start: 1e6 end: 1e12 y: - quantity: "snia_rates.snia_rates_50_kpc" - units: "1/year" - start: 3e-8 - end: 1e-1 + quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" + units: "1/year/Solar_Mass" + start: 3e-15 + end: 3e-12 mean: plot: true log: true @@ -61,15 +63,15 @@ stellar_mass_snia_rates_active_only_50: value: 1e12 units: Solar_Mass metadata: - title: Stellar Mass-SNIa rate relation (active only, 50 kpc aperture) - caption: Uses a 50 kpc 3D aperture and only active galaxies. - section: SNIa Rate + title: Stellar Mass-SNIa rate per stellar mass relation (active only, 50 kpc aperture) + caption: Uses a 50 kpc 3D aperture. + section: SNIa Rate (Stellar Mass) observational_data: - - filename: GalaxyStellarMassSNIaRate/Smith2012_active.hdf5 - - filename: GalaxyStellarMassSNIaRate/Graur2015_active.hdf5 - - filename: GalaxyStellarMassSNIaRate/Graur2017_active.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Smith2012_active.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2015_active.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2017_active.hdf5 -stellar_mass_snia_rates_passive_only_50: +stellar_mass_snia_rates_per_stellar_mass_passive_only_50: type: "scatter" comment: "Passive only" comment_loc: "lower left" @@ -81,10 +83,10 @@ stellar_mass_snia_rates_passive_only_50: start: 1e6 end: 1e12 y: - quantity: "snia_rates.snia_rates_50_kpc" - units: "1/year" - start: 3e-8 - end: 1e-1 + quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" + units: "1/year/Solar_Mass" + start: 3e-15 + end: 3e-12 mean: plot: true log: true @@ -97,15 +99,15 @@ stellar_mass_snia_rates_passive_only_50: value: 1e12 units: Solar_Mass metadata: - title: Stellar Mass-SNIa rate relation (passive only, 50 kpc aperture) - caption: Uses a 50 kpc 3D aperture and only passive galaxies. - section: SNIa Rate + title: Stellar Mass-SNIa rate per stellar mass relation (passive only, 50 kpc aperture) + caption: Uses a 50 kpc 3D aperture. + section: SNIa Rate (Stellar Mass) observational_data: - - filename: GalaxyStellarMassSNIaRate/Smith2012_passive.hdf5 - - filename: GalaxyStellarMassSNIaRate/Graur2015_passive.hdf5 - - filename: GalaxyStellarMassSNIaRate/Graur2017_passive.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Smith2012_passive.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2015_passive.hdf5 + - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2017_passive.hdf5 -stellar_mass_snia_rates_per_stellar_mass_50: +stellar_mass_snia_rates_50: type: "scatter" legend_loc: "upper left" x: @@ -114,10 +116,10 @@ stellar_mass_snia_rates_per_stellar_mass_50: start: 1e6 end: 1e12 y: - quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" - units: "1/year/Solar_Mass" - start: 3e-15 - end: 3e-12 + quantity: "snia_rates.snia_rates_50_kpc" + units: "1/year" + start: 3e-8 + end: 1e-1 mean: plot: true log: true @@ -130,17 +132,17 @@ stellar_mass_snia_rates_per_stellar_mass_50: value: 1e12 units: Solar_Mass metadata: - title: Stellar Mass-SNIa rate per stellar mass relation (50 kpc aperture) + title: Stellar Mass-SNIa rate relation (50 kpc aperture) caption: Uses a 50 kpc 3D aperture. - section: SNIa Rate + section: SNIa Rate (Stellar Mass) observational_data_bracket_width: 10.0 observational_data: - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Kistler2014.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2015.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2017.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Wiseman2021.hdf5 + - filename: GalaxyStellarMassSNIaRate/Kistler2014.hdf5 + - filename: GalaxyStellarMassSNIaRate/Graur2015.hdf5 + - filename: GalaxyStellarMassSNIaRate/Graur2017.hdf5 + - filename: GalaxyStellarMassSNIaRate/Wiseman2021.hdf5 -stellar_mass_snia_rates_per_stellar_mass_active_only_50: +stellar_mass_snia_rates_active_only_50: type: "scatter" comment: "Active only" comment_loc: "lower left" @@ -152,10 +154,10 @@ stellar_mass_snia_rates_per_stellar_mass_active_only_50: start: 1e6 end: 1e12 y: - quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" - units: "1/year/Solar_Mass" - start: 3e-15 - end: 3e-12 + quantity: "snia_rates.snia_rates_50_kpc" + units: "1/year" + start: 3e-8 + end: 1e-1 mean: plot: true log: true @@ -168,15 +170,15 @@ stellar_mass_snia_rates_per_stellar_mass_active_only_50: value: 1e12 units: Solar_Mass metadata: - title: Stellar Mass-SNIa rate per stellar mass relation (active only, 50 kpc aperture) - caption: Uses a 50 kpc 3D aperture. - section: SNIa Rate + title: Stellar Mass-SNIa rate relation (active only, 50 kpc aperture) + caption: Uses a 50 kpc 3D aperture and only active galaxies. + section: SNIa Rate (Stellar Mass) observational_data: - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Smith2012_active.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2015_active.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2017_active.hdf5 + - filename: GalaxyStellarMassSNIaRate/Smith2012_active.hdf5 + - filename: GalaxyStellarMassSNIaRate/Graur2015_active.hdf5 + - filename: GalaxyStellarMassSNIaRate/Graur2017_active.hdf5 -stellar_mass_snia_rates_per_stellar_mass_passive_only_50: +stellar_mass_snia_rates_passive_only_50: type: "scatter" comment: "Passive only" comment_loc: "lower left" @@ -188,10 +190,10 @@ stellar_mass_snia_rates_per_stellar_mass_passive_only_50: start: 1e6 end: 1e12 y: - quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" - units: "1/year/Solar_Mass" - start: 3e-15 - end: 3e-12 + quantity: "snia_rates.snia_rates_50_kpc" + units: "1/year" + start: 3e-8 + end: 1e-1 mean: plot: true log: true @@ -204,13 +206,13 @@ stellar_mass_snia_rates_per_stellar_mass_passive_only_50: value: 1e12 units: Solar_Mass metadata: - title: Stellar Mass-SNIa rate per stellar mass relation (passive only, 50 kpc aperture) - caption: Uses a 50 kpc 3D aperture. - section: SNIa Rate + title: Stellar Mass-SNIa rate relation (passive only, 50 kpc aperture) + caption: Uses a 50 kpc 3D aperture and only passive galaxies. + section: SNIa Rate (Stellar Mass) observational_data: - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Smith2012_passive.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2015_passive.hdf5 - - filename: GalaxyStellarMassSNIaRatePerStellarMass/Graur2017_passive.hdf5 + - filename: GalaxyStellarMassSNIaRate/Smith2012_passive.hdf5 + - filename: GalaxyStellarMassSNIaRate/Graur2015_passive.hdf5 + - filename: GalaxyStellarMassSNIaRate/Graur2017_passive.hdf5 gas_metallicity_snia_rates_per_stellar_mass_active_only_50: comment: "$M_\\star > 10^{10}$ $M_\\odot$ and active only" @@ -243,7 +245,7 @@ gas_metallicity_snia_rates_per_stellar_mass_active_only_50: metadata: title: Gas metallicity-SNIa rate per stellar mass relation (active only, 50 kpc aperture) caption: Uses a 50 kpc 3D aperture, active galaxies only with a stellar mass above 1e10 Msun. - section: SNIa Rate + section: SNIa Rate (Gas Metallicity) observational_data: - filename: GalaxyGasMetallicitySNIaRatePerStellarMass/Graur2017.hdf5 @@ -278,11 +280,11 @@ gas_metallicity_snia_rates_per_stellar_mass_active_only_50_Mstar5e10: metadata: title: Gas metallicity-SNIa rate per stellar mass relation (active only, 50 kpc aperture) caption: Uses a 50 kpc 3D aperture, active galaxies only with a stellar mass above 5e10 Msun. - section: SNIa Rate + section: SNIa Rate (Gas Metallicity) observational_data: - filename: GalaxyGasMetallicitySNIaRatePerStellarMass/Graur2017.hdf5 -star_formation_rates_snia_rates_per_stellar_mass_50: +star_formation_rates_snia_rates_50: comment: "$M_\\star > 10^{10}$ $M_\\odot$" comment_loc: "lower left" type: "scatter" @@ -294,10 +296,10 @@ star_formation_rates_snia_rates_per_stellar_mass_50: start: 1e-3 end: 1e2 y: - quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" - units: "1/year/Solar_Mass" - start: 1e-14 - end: 3e-12 + quantity: "snia_rates.snia_rates_50_kpc" + units: "1/year" + start: 1e-6 + end: 1e-1 mean: plot: true log: true @@ -310,19 +312,19 @@ star_formation_rates_snia_rates_per_stellar_mass_50: value: 1e2 units: "Solar_Mass/year" metadata: - title: SFR-SNIa rate per stellar mass relation (50 kpc aperture) + title: SFR-SNIa rate relation (50 kpc aperture) caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above 1e10 Msun. - section: SNIa Rate + section: SNIa Rate (SFR and sSFR) observational_data: - - filename: StarFormationRateSNIaRatePerStellarMass/Graur2015.hdf5 - - filename: StarFormationRateSNIaRatePerStellarMass/Graur2017.hdf5 + - filename: StarFormationRateSNIaRate/Sullivan2006.hdf5 + - filename: StarFormationRateSNIaRate/Smith2012.hdf5 -star_formation_rates_snia_rates_per_stellar_mass_50_Mstar5e10: - comment: "$M_\\star > 5 \\times 10^{10}$ $M_\\odot$" +star_formation_rates_snia_rates_per_stellar_mass_50: + comment: "$M_\\star > 10^{10}$ $M_\\odot$" comment_loc: "lower left" type: "scatter" legend_loc: "upper left" - selection_mask: "derived_quantities.stellar_mass_is_bigger_than_5e10_msun_50_kpc" + selection_mask: "derived_quantities.stellar_mass_is_bigger_than_1e10_msun_50_kpc" x: quantity: "apertures.sfr_gas_50_kpc" units: "Solar_Mass/year" @@ -346,15 +348,56 @@ star_formation_rates_snia_rates_per_stellar_mass_50_Mstar5e10: units: "Solar_Mass/year" metadata: title: SFR-SNIa rate per stellar mass relation (50 kpc aperture) - caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above 5e10 Msun. - section: SNIa Rate + caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above $10^{10} \mathrm{M_\odot}$. + section: SNIa Rate (SFR and sSFR) observational_data: - filename: StarFormationRateSNIaRatePerStellarMass/Graur2015.hdf5 - filename: StarFormationRateSNIaRatePerStellarMass/Graur2017.hdf5 -star_formation_rates_snia_rates_50: +specific_star_formation_rates_snia_rates_per_stellar_mass_50: + comment: "$M_\\star > 10^{10}$ $M_\\odot$" + comment_loc: "lower left" type: "scatter" legend_loc: "upper left" + selection_mask: "derived_quantities.stellar_mass_is_bigger_than_1e10_msun_50_kpc" + x: + quantity: "derived_quantities.specific_sfr_gas_50_kpc" + units: "1/gigayear" + start: 1e-3 + end: 1e1 + y: + quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" + units: "1/year/Solar_Mass" + start: 1e-14 + end: 3e-12 + mean: + plot: true + log: true + adaptive: true + number_of_bins: 30 + start: + value: 1e-4 + units: "1/gigayear" + end: + value: 1e1 + units: "1/gigayear" + metadata: + title: sSFR-SNIa rate per stellar mass relation (50 kpc aperture) + caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above $10^{10} \mathrm{M_\odot}$. + section: SNIa Rate (SFR and sSFR) + observational_data: + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Mannucci2005.hdf5 + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Sullivan2006.hdf5 + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Smith2012.hdf5 + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Graur2015.hdf5 + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Graur2017.hdf5 + +star_formation_rates_snia_rates_50_Mstar5e10: + comment: "$M_\\star > 5 \\times 10^{10}$ $M_\\odot$" + comment_loc: "lower left" + type: "scatter" + legend_loc: "upper left" + selection_mask: "derived_quantities.stellar_mass_is_bigger_than_5e10_msun_50_kpc" x: quantity: "apertures.sfr_gas_50_kpc" units: "Solar_Mass/year" @@ -378,22 +421,23 @@ star_formation_rates_snia_rates_50: units: "Solar_Mass/year" metadata: title: SFR-SNIa rate relation (50 kpc aperture) - caption: Uses a 50 kpc 3D aperture. - section: SNIa Rate + caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above 1e10 Msun. + section: SNIa Rate (SFR and sSFR) observational_data: + - filename: StarFormationRateSNIaRate/Sullivan2006.hdf5 - filename: StarFormationRateSNIaRate/Smith2012.hdf5 -specific_star_formation_rates_snia_rates_per_stellar_mass_50: - comment: "$M_\\star > 10^{10}$ $M_\\odot$" +star_formation_rates_snia_rates_per_stellar_mass_50_Mstar5e10: + comment: "$M_\\star > 5 \\times 10^{10}$ $M_\\odot$" comment_loc: "lower left" type: "scatter" legend_loc: "upper left" - selection_mask: "derived_quantities.stellar_mass_is_bigger_than_1e10_msun_50_kpc" + selection_mask: "derived_quantities.stellar_mass_is_bigger_than_5e10_msun_50_kpc" x: - quantity: "derived_quantities.specific_sfr_gas_50_kpc" - units: "1/gigayear" + quantity: "apertures.sfr_gas_50_kpc" + units: "Solar_Mass/year" start: 1e-3 - end: 1e1 + end: 1e2 y: quantity: "derived_quantities.snia_rate_per_stellar_mass_50_kpc" units: "1/year/Solar_Mass" @@ -406,17 +450,17 @@ specific_star_formation_rates_snia_rates_per_stellar_mass_50: number_of_bins: 30 start: value: 1e-3 - units: "1/gigayear" + units: "Solar_Mass/year" end: - value: 1e1 - units: "1/gigayear" + value: 1e2 + units: "Solar_Mass/year" metadata: - title: sSFR-SNIa rate per stellar mass relation (50 kpc aperture) - caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above 1e10 Msun. - section: SNIa Rate + title: SFR-SNIa rate per stellar mass relation (50 kpc aperture) + caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above $5 \times 10^{10} \mathrm{M_\odot}$. + section: SNIa Rate (SFR and sSFR) observational_data: - - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Graur2015.hdf5 - - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Graur2017.hdf5 + - filename: StarFormationRateSNIaRatePerStellarMass/Graur2015.hdf5 + - filename: StarFormationRateSNIaRatePerStellarMass/Graur2017.hdf5 specific_star_formation_rates_snia_rates_per_stellar_mass_50_Mstar5e10: comment: "$M_\\star > 5 \\times 10^{10}$ $M_\\odot$" @@ -440,7 +484,7 @@ specific_star_formation_rates_snia_rates_per_stellar_mass_50_Mstar5e10: adaptive: true number_of_bins: 30 start: - value: 1e-3 + value: 1e-4 units: "1/gigayear" end: value: 1e1 @@ -448,8 +492,11 @@ specific_star_formation_rates_snia_rates_per_stellar_mass_50_Mstar5e10: metadata: title: sSFR-SNIa rate per stellar mass relation (50 kpc aperture) caption: Uses a 50 kpc 3D aperture, only select galaxies with a stellar mass above 5e10 Msun. - section: SNIa Rate + section: SNIa Rate (SFR and sSFR) observational_data: + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Mannucci2005.hdf5 + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Sullivan2006.hdf5 + - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Smith2012.hdf5 - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Graur2015.hdf5 - filename: SpecificStarFormationRateSNIaRatePerStellarMass/Graur2017.hdf5 @@ -480,7 +527,7 @@ stellar_mass_snia_rates_30: metadata: title: Stellar Mass-SNIa rate relation (50 kpc aperture) caption: Uses a 30 kpc 3D aperture. - section: SNIa Rate + section: SNIa Rate (Stellar Mass) show_on_webpage: false observational_data: - filename: GalaxyStellarMassSNIaRate/Wiseman2021.hdf5 @@ -512,7 +559,7 @@ stellar_mass_snia_rates_per_stellar_mass_30: metadata: title: Stellar Mass-SNIa rate per stellar mass relation (30 kpc aperture) caption: Uses a 30 kpc 3D aperture. - section: SNIa Rate + section: SNIa Rate (Stellar Mass) show_on_webpage: false observational_data: - filename: GalaxyStellarMassSNIaRatePerStellarMass/Wiseman2021.hdf5 @@ -546,7 +593,7 @@ stellar_mass_snia_rates_100: metadata: title: Stellar Mass-SNIa rate relation (100 kpc aperture) caption: Uses a 30 kpc 3D aperture. - section: SNIa Rate + section: SNIa Rate (Stellar Mass) show_on_webpage: false observational_data: - filename: GalaxyStellarMassSNIaRate/Wiseman2021.hdf5 @@ -578,7 +625,7 @@ stellar_mass_snia_rates_per_stellar_mass_100: metadata: title: Stellar Mass-SNIa rate per stellar mass relation (100 kpc aperture) caption: Uses a 30 kpc 3D aperture. - section: SNIa Rate + section: SNIa Rate (Stellar Mass) show_on_webpage: false observational_data: - filename: GalaxyStellarMassSNIaRatePerStellarMass/Wiseman2021.hdf5 diff --git a/colibre/config.yml b/colibre/config.yml index d11c30e6..958ea354 100644 --- a/colibre/config.yml +++ b/colibre/config.yml @@ -20,6 +20,66 @@ description_template: description.html # Location and description of additional figures scripts: + - filename: scripts/cosmic_log_SNIa_rate_CSFH_based.py + title: 'Cosmic SNIa rate' + caption: SNIa rate history calculated from the SFR.txt file produced by SWIFT. + section: SN Rate (cosmic - CSFH based) + output_file: log_SNIa_rate_history_based_on_CSFH.png + - filename: scripts/cosmic_SNIa_rate_CSFH_based.py + title: 'Cosmic SNIa rate' + caption: SNIa rate history calculated from the SFR.txt file produced by SWIFT. + section: SN Rate (cosmic - CSFH based) + output_file: SNIa_rate_history_based_on_CSFH.png + - filename: scripts/cosmic_log_CC_SN_rate_CSFH_based.py + title: 'Cosmic CC SN rate' + caption: CC SN rate history calculated from the SFR.txt file produced by SWIFT. + section: SN Rate (cosmic - CSFH based) + output_file: log_CC_SN_rate_history_based_on_CSFH.png + - filename: scripts/cosmic_CC_SN_rate_CSFH_based.py + title: 'Cosmic CC SN rate' + caption: CC SN rate history calculated from the SFR.txt file produced by SWIFT. + section: SN Rate (cosmic - CSFH based) + output_file: CC_SN_rate_history_based_on_CSFH.png + - filename: scripts/cosmic_SNIa_over_CC_SN_rate_CSFH_based.py + title: 'Cosmic SNIa rate / CC SN rate' + caption: SNIa rate over CC SN rate history calculated from the SFR.txt file produced by SWIFT. + section: SN Rate (cosmic - CSFH based) + output_file: SNIa_over_CC_rate_history_based_on_CSFH.png + - filename: scripts/cosmic_SNIa_rate.py + title: 'Cosmic SNIa rate' + caption: SNIa rate history plotted directly from the SNIa.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: SNIa_rate_history.png + - filename: scripts/cosmic_log_SNIa_rate.py + title: 'Cosmic SNIa rate' + caption: SNIa rate history plotted directly from the SNIa.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: log_SNIa_rate_history.png + - filename: scripts/cosmic_CC_SN_rate.py + title: 'Cosmic CC SN rate' + caption: CC SN rate history plotted directly from the SNII.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: CC_SN_rate_history.png + - filename: scripts/cosmic_log_CC_SN_rate.py + title: 'Cosmic CC SN rate' + caption: CC SN rate history plotted directly from the SNII.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: log_CC_SN_rate_history.png + - filename: scripts/cosmic_log_NSM_rate.py + title: 'Cosmic NS-NS merger rate' + caption: NS-NS merger rate or short GRB rate plotted directly from the r_processes.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: log_NSM_rate_history.png + - filename: scripts/cosmic_log_CEJSN_rate.py + title: 'Cosmic common-envelop jets SN rate' + caption: Common-envelop (CE) jets SN rate plotted directly from the r_processes.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: log_CEJSN_rate_history.png + - filename: scripts/cosmic_log_collapsar_rate.py + title: 'Cosmic collapsar rate' + caption: Collapsar rate plotted directly from the r_processes.txt file produced by SWIFT. + section: SN Rate (cosmic - stochastic) + output_file: log_collapsar_rate_history.png - filename: scripts/density_temperature.py caption: Density-temperature diagram. If present, dashed line represents the entropy floor (equation of state). output_file: density_temperature.png @@ -503,8 +563,3 @@ scripts: section: Column Densities additional_arguments: parallel: True - - filename: scripts/cosmic_SNIa_rate.py - title: 'Cosmic SNIa rate' - caption: SNIa rate history plotted directly from the SNIa.txt file produced by SWIFT. - section: SNIa Rate - output_file: SNIa_rate_history.png diff --git a/colibre/description.html b/colibre/description.html index eb1ce987..c7124f03 100644 --- a/colibre/description.html +++ b/colibre/description.html @@ -379,6 +379,99 @@

Chemistry

+

SNIa DTD

+{% if data.metadata.parameters.get("SNIaDTD:normalization_timescale_Gyr", False) %} + + + + + + {% if data.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", False) %} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + {% else %} + + + + + + + + + + + + + + + + + + + {% if data.metadata.parameters.get("SNIaDTD:power_law_slope", False) %} + + {% else %} + + {% endif %} + + {% endif %} +
ParameterValue
ModelGaussian
SNIa efficiency per Msun (Gaussian part){{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_efficiency_gauss_p_Msun") }}
SNIa efficiency per Msun (constant part){{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_efficiency_const_p_Msun") }}
SNIa delay time{{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_delay_time_Gyr") }} Gyr
Normalisation time scale (constant part){{ data.metadata.parameters | get_if_present_float("SNIaDTD:normalization_timescale_Gyr") }} Gyr
Characteristic time of the Guassian{{ data.metadata.parameters | get_if_present_float("SNIaDTD:characteristic_time_Gyr") }} Gyr
Standard deviation time of the Guassian{{ data.metadata.parameters | get_if_present_float("SNIaDTD:STD_characteristic_time_Gyr") }} Gyr
ModelPower law
SNIa efficiency per Msun{{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_efficiency_p_Msun") }}
SNIa delay time{{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_delay_time_Gyr") }} Gyr
Normalisation time scale{{ data.metadata.parameters | get_if_present_float("SNIaDTD:normalization_timescale_Gyr") }} Gyr
Power law slope{{ data.metadata.parameters | get_if_present_float("SNIaDTD:power_law_slope") }}1.0
+ +{% else %} +{% if data.metadata.parameters.get("SNIaDTD:SNIa_timescale_Gyr", False) %} + + + + + + + + + + + + + + + + + + + + + +
ParameterValue
Modelexponential
SNIa efficiency per Msun{{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_efficiency_p_Msun") }}
SNIa delay time{{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_delay_time_Gyr") }} Gyr
Exponential delay time{{ data.metadata.parameters | get_if_present_float("SNIaDTD:SNIa_timescale_Gyr") }} Gyr
+{% else %} +

DTD model is not found

+{% endif %} +{% endif %} +

AGN feedback

{% if data.metadata.parameters.get("COLIBREAGN:AGN_delta_T_K", False) %} diff --git a/colibre/scripts/cosmic_CC_SN_rate.py b/colibre/scripts/cosmic_CC_SN_rate.py new file mode 100644 index 00000000..047878f4 --- /dev/null +++ b/colibre/scripts/cosmic_CC_SN_rate.py @@ -0,0 +1,162 @@ +""" +Plots the cosmic CCSN rate history. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +from swiftpipeline.argumentparser import ScriptArgumentParser + +arguments = ScriptArgumentParser( + description="Creates a CC SN rate history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +CCSN_filenames = [f"{directory}/SNII.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.semilogx() + +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +CC_SN_rate_output_units = 1.0 / (unyt.yr * unyt.Mpc ** 3) + +for idx, (snapshot_filename, CCSN_filename, name) in enumerate( + zip(snapshot_filenames, CCSN_filenames, names) +): + data = np.loadtxt( + CCSN_filename, + usecols=(4, 6, 11), + dtype=[("a", np.float32), ("z", np.float32), ("CC SN rate", np.float32)], + ) + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + CC_SN_rate_units = 1.0 / (units.time * units.length ** 3) + + # a, Redshift, SFR + scale_factor = data["a"] + CC_SN_rate = (data["CC SN rate"] * CC_SN_rate_units).to(CC_SN_rate_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factor, CC_SN_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicCCSNRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + ) + ) + observation_labels.append(f"{obs_data.citation}") + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc="center right", fontsize="xx-small" +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(0.0, 26.0) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel( + f"CC SN rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]" +) +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/CC_SN_rate_history.png") + diff --git a/colibre/scripts/cosmic_CC_SN_rate_CSFH_based.py b/colibre/scripts/cosmic_CC_SN_rate_CSFH_based.py new file mode 100644 index 00000000..6199251c --- /dev/null +++ b/colibre/scripts/cosmic_CC_SN_rate_CSFH_based.py @@ -0,0 +1,191 @@ +""" +Plots the star formation history. Modified version of the script in the +github.com/swiftsim/swiftsimio-examples repository. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from load_sfh_data import read_obs_data + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +sfr_output_units = unyt.msun / (unyt.year * unyt.Mpc ** 3) +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +SNIa_output_units = 1. / (unyt.year * unyt.Mpc ** 3) + +from swiftpipeline.argumentparser import ScriptArgumentParser + +def CC_SN_DTD(t, t_min, t_max): + mask = (t > t_min) & (t < t_max) + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1./(t_max - t_min) * (np.ones(len(t[mask]))) + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +arguments = ScriptArgumentParser( + description="Creates a CC SN history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +sfr_filenames = [f"{directory}/SFR.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.set_xscale("log") + +for idx, (snapshot_filename, sfr_filename, name) in enumerate( + zip(snapshot_filenames, sfr_filenames, names) +): + data = np.genfromtxt(sfr_filename).T + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + boxsize = snapshot.metadata.boxsize + box_volume = boxsize[0] * boxsize[1] * boxsize[2] + + sfr_units = snapshot.gas.star_formation_rates.units + + # a, Redshift, SFR + scale_factor = data[2] + redshift = data[3] + star_formation_rate = (data[7] * sfr_units / box_volume).to(sfr_output_units) + times = data[1] * units.time + dM = data[4] * units.mass + + # Calculate the times we plot + times_desired = np.linspace(0.05,13.8,500) * unyt.Gyr + scale_factors_use = np.interp(times_desired, times, scale_factor) + + SNIa_rate = np.zeros(len(times_desired)) + + for i in range(1,len(times_desired)): + time_consider = times[times < times_desired[i]] + time_since_formation = times_desired[i] - time_consider + dM_consider = dM[times < times_desired[i]] + SNIa_rate_individual = 1.180e-2 * CC_SN_DTD(time_since_formation, 3*unyt.Myr, 0.04*unyt.Gyr) * dM_consider / unyt.Msun + + SNIa_rate_sum = np.sum(SNIa_rate_individual) + SNIa_rate[i] = SNIa_rate_sum + + SNIa_rate = (SNIa_rate / box_volume / unyt.Gyr).to(SNIa_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factors_use, SNIa_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicCCSNRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + capsize=1.0, + ) + ) + observation_labels.append(f"{obs_data.citation}") + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc=4, fontsize=4, ncol=2 +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(0, 26.0) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel(f"CC SN rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]") +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/CC_SN_rate_history_based_on_CSFH.png") diff --git a/colibre/scripts/cosmic_SNIa_over_CC_SN_rate_CSFH_based.py b/colibre/scripts/cosmic_SNIa_over_CC_SN_rate_CSFH_based.py new file mode 100644 index 00000000..ca15dce7 --- /dev/null +++ b/colibre/scripts/cosmic_SNIa_over_CC_SN_rate_CSFH_based.py @@ -0,0 +1,254 @@ +""" +Plots the cosmic SNIa rate over the core-collapse SNe rate based on the star formation history. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from load_sfh_data import read_obs_data + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +sfr_output_units = unyt.msun / (unyt.year * unyt.Mpc ** 3) +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +SNIa_output_units = 1. / (unyt.year * unyt.Mpc ** 3) + +from swiftpipeline.argumentparser import ScriptArgumentParser + +def power_law_beta_one_DTD(t, t_delay, tH): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1./(np.log(tH) - np.log(t_delay)) / t[mask] + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def power_law_DTD(t, t_delay, tH, beta): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = (1.-beta)/(tH**(1-beta) - (t_delay)**(1-beta)) / t[mask]**beta + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + + +def exponential_law_DTD(t, t_decay, t_delay): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = np.exp(-(t[mask]-t_delay)/t_decay)/t_decay + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def Gaussian_law_DTD(t, t_delay, tmean, tsigma): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1/np.sqrt(2*tsigma**2) * np.exp(-0.5*(t[mask]-tmean)**2 / (2*tsigma)**2) + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def CC_SN_DTD(t, t_min, t_max): + mask = (t > t_min) & (t < t_max) + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1./(t_max - t_min) * (np.ones(len(t[mask]))) + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +arguments = ScriptArgumentParser( + description="Creates a star formation history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +sfr_filenames = [f"{directory}/SFR.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.set_xscale("log") + + +for idx, (snapshot_filename, sfr_filename, name) in enumerate( + zip(snapshot_filenames, sfr_filenames, names) +): + data = np.genfromtxt(sfr_filename).T + + snapshot = load(snapshot_filename) + + # find out which DTD model we are currently using in the code + normalization_timescale = snapshot.metadata.parameters.get("SNIaDTD:normalization_timescale_Gyr", None) + if normalization_timescale is None: + have_normalization_timescale = False + else: + normalization_timescale = float(normalization_timescale) * unyt.Gyr + have_normalization_timescale = True + + exponential_decay = snapshot.metadata.parameters.get("SNIaDTD:SNIa_timescale_Gyr", None) + if exponential_decay is None: + have_exponential_decay = False + else: + have_exponential_decay = True + exponential_decay = float(exponential_decay) * unyt.Gyr + + Gaussian_SNIa_efficiency = snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", None) + if Gaussian_SNIa_efficiency is None: + have_Gaussian = False + else: + have_Gaussian = True + Gaussian_SNIa_efficiency = float(Gaussian_SNIa_efficiency) / unyt.Msun + + power_law_slope = snapshot.metadata.parameters.get("SNIaDTD:power_law_slope", None) + if power_law_slope is None: + have_slope = False + else: + have_slope = True + power_law_slope = float(power_law_slope) + + if have_Gaussian: + used_DTD = "Gaussian" + elif have_slope and have_normalization_timescale: + used_DTD = "power-law" + elif have_normalization_timescale: + used_DTD = "power-law-beta-one" + else: + used_DTD = "exponential" + + delay_time = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_delay_time_Gyr", False) ) * unyt.Gyr + + if used_DTD in ["power-law", "exponential", "power-law-beta-one"]: + SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_p_Msun", False) ) / unyt.Msun + elif used_DTD == "Gaussian": + Gaussian_SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", False)) / unyt.Msun + Gaussian_const_SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_const_p_Msun", False)) / unyt.Msun + Gaussian_characteristic_time = float(snapshot.metadata.parameters.get("SNIaDTD:characteristic_time_Gyr", False)) * unyt.Gyr + Gaussian_std_time = float(snapshot.metadata.parameters.get("SNIaDTD:STD_characteristic_time_Gyr", False)) * unyt.Gyr + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + boxsize = snapshot.metadata.boxsize + box_volume = boxsize[0] * boxsize[1] * boxsize[2] + + sfr_units = snapshot.gas.star_formation_rates.units + + # a, Redshift, SFR + scale_factor = data[2] + redshift = data[3] + star_formation_rate = (data[7] * sfr_units / box_volume).to(sfr_output_units) + times = data[1] * units.time + dM = data[4] * units.mass + + # Calculate the times we plot + times_desired = np.linspace(0.05,13.8,500) * unyt.Gyr + scale_factors_use = np.interp(times_desired, times, scale_factor) + + SNIa_rate = np.zeros(len(times_desired)) + CC_SN_rate = np.zeros(len(times_desired)) + + for i in range(1,len(times_desired)): + time_consider = times[times < times_desired[i]] + time_since_formation = times_desired[i] - time_consider + dM_consider = dM[times < times_desired[i]] + if used_DTD == "power-law": + SNIa_rate_individual = SNIa_efficiency * power_law_DTD(time_since_formation, delay_time, normalization_timescale,power_law_slope) * dM_consider + elif used_DTD == "power-law-beta-one": + SNIa_rate_individual = SNIa_efficiency * power_law_beta_one_DTD(time_since_formation, delay_time, normalization_timescale) * dM_consider + elif used_DTD == "exponential": + SNIa_rate_individual = SNIa_efficiency * exponential_law_DTD(time_since_formation, exponential_decay, delay_time) * dM_consider + elif used_DTD == "Gaussian": + SNIa_rate_individual = Gaussian_SNIa_efficiency * Gaussian_law_DTD(time_since_formation, delay_time, Gaussian_characteristic_time, Gaussian_std_time) * dM_consider + + SNIa_rate_sum = np.sum(SNIa_rate_individual) + SNIa_rate[i] = SNIa_rate_sum + + CC_SN_rate_individual = 1.180e-2 * CC_SN_DTD(time_since_formation, 3*unyt.Myr, 0.04*unyt.Gyr) * dM_consider / unyt.Msun + + CC_SN_rate_sum = np.sum(CC_SN_rate_individual) + CC_SN_rate[i] = CC_SN_rate_sum + + + SNIa_rate = (SNIa_rate / box_volume / unyt.Gyr).to(SNIa_output_units) + CC_SN_rate = (CC_SN_rate / box_volume / unyt.Gyr).to(SNIa_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factors_use, SNIa_rate.value / CC_SN_rate.value, zorder=10000)[0] + ) + simulation_labels.append(name) + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(0, 0.6) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel("SNIa rate / CC SN rate") +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/SNIa_over_CC_rate_history_based_on_CSFH.png") diff --git a/colibre/scripts/cosmic_SNIa_rate.py b/colibre/scripts/cosmic_SNIa_rate.py index b440b3d1..0f5c77fa 100644 --- a/colibre/scripts/cosmic_SNIa_rate.py +++ b/colibre/scripts/cosmic_SNIa_rate.py @@ -93,6 +93,7 @@ markeredgecolor="none", markersize=2, zorder=-10, + capsize=1.0, ) ) observation_labels.append(f"{obs_data.citation}") @@ -117,7 +118,7 @@ ax.set_xticklabels(redshift_labels) observation_legend = ax.legend( - observation_lines, observation_labels, markerfirst=True, loc="center right" + observation_lines, observation_labels, markerfirst=True, loc="lower right", fontsize="xx-small", ncol=2 ) simulation_legend = ax.legend( diff --git a/colibre/scripts/cosmic_SNIa_rate_CSFH_based.py b/colibre/scripts/cosmic_SNIa_rate_CSFH_based.py new file mode 100644 index 00000000..84570cf3 --- /dev/null +++ b/colibre/scripts/cosmic_SNIa_rate_CSFH_based.py @@ -0,0 +1,273 @@ +""" +Plots the cosmic SNIa rate based on the star formation history. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from load_sfh_data import read_obs_data + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +sfr_output_units = unyt.msun / (unyt.year * unyt.Mpc ** 3) +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +SNIa_output_units = 1. / (unyt.year * unyt.Mpc ** 3) + +from swiftpipeline.argumentparser import ScriptArgumentParser + +def power_law_beta_one_DTD(t, t_delay, tH): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1./(np.log(tH) - np.log(t_delay)) / t[mask] + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def power_law_DTD(t, t_delay, tH, beta): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = (1.-beta)/(tH**(1-beta) - (t_delay)**(1-beta)) / t[mask]**beta + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + + +def exponential_law_DTD(t, t_decay, t_delay): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = np.exp(-(t[mask]-t_delay)/t_decay)/t_decay + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def Gaussian_law_DTD(t, t_delay, tmean, tsigma): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + norm = 1./(tsigma* np.sqrt(2*np.pi)) + delta_t = t[mask] - tmean + value_new = norm * np.exp(-0.5 * delta_t**2 / tsigma**2) + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +arguments = ScriptArgumentParser( + description="Creates a cosmic SNIa rate history plot based on the star formation history, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +sfr_filenames = [f"{directory}/SFR.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.set_xscale("log") + + +for idx, (snapshot_filename, sfr_filename, name) in enumerate( + zip(snapshot_filenames, sfr_filenames, names) +): + data = np.genfromtxt(sfr_filename).T + + snapshot = load(snapshot_filename) + + # find out which DTD model we are currently using in the code + normalization_timescale = snapshot.metadata.parameters.get("SNIaDTD:normalization_timescale_Gyr", None) + if normalization_timescale is None: + have_normalization_timescale = False + else: + normalization_timescale = float(normalization_timescale) * unyt.Gyr + have_normalization_timescale = True + + exponential_decay = snapshot.metadata.parameters.get("SNIaDTD:SNIa_timescale_Gyr", None) + if exponential_decay is None: + have_exponential_decay = False + else: + have_exponential_decay = True + exponential_decay = float(exponential_decay) * unyt.Gyr + + Gaussian_SNIa_efficiency = snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", None) + if Gaussian_SNIa_efficiency is None: + have_Gaussian = False + else: + have_Gaussian = True + Gaussian_SNIa_efficiency = float(Gaussian_SNIa_efficiency) / unyt.Msun + + power_law_slope = snapshot.metadata.parameters.get("SNIaDTD:power_law_slope", None) + if power_law_slope is None: + have_slope = False + else: + have_slope = True + power_law_slope = float(power_law_slope) + + if have_Gaussian: + used_DTD = "Gaussian" + elif have_slope and have_normalization_timescale: + used_DTD = "power-law" + elif have_normalization_timescale: + used_DTD = "power-law-beta-one" + else: + used_DTD = "exponential" + + delay_time = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_delay_time_Gyr", False) ) * unyt.Gyr + + if used_DTD in ["power-law", "exponential", "power-law-beta-one"]: + SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_p_Msun", False) ) / unyt.Msun + elif used_DTD == "Gaussian": + Gaussian_SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", False)) / unyt.Msun + Gaussian_const_SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_const_p_Msun", False)) / unyt.Msun + Gaussian_characteristic_time = float(snapshot.metadata.parameters.get("SNIaDTD:characteristic_time_Gyr", False)) * unyt.Gyr + Gaussian_std_time = float(snapshot.metadata.parameters.get("SNIaDTD:STD_characteristic_time_Gyr", False)) * unyt.Gyr + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + boxsize = snapshot.metadata.boxsize + box_volume = boxsize[0] * boxsize[1] * boxsize[2] + + sfr_units = snapshot.gas.star_formation_rates.units + + # a, Redshift, SFR + scale_factor = data[2] + redshift = data[3] + star_formation_rate = (data[7] * sfr_units / box_volume).to(sfr_output_units) + times = data[1] * units.time + dM = data[4] * units.mass + + # Calculate the times we plot + times_desired = np.linspace(0.05,13.8,500) * unyt.Gyr + scale_factors_use = np.interp(times_desired, times, scale_factor) + + SNIa_rate = np.zeros(len(times_desired)) + + for i in range(1,len(times_desired)): + time_consider = times[times < times_desired[i]] + time_since_formation = times_desired[i] - time_consider + dM_consider = dM[times < times_desired[i]] + if used_DTD == "power-law": + SNIa_rate_individual = SNIa_efficiency * power_law_DTD(time_since_formation, delay_time, normalization_timescale,power_law_slope) * dM_consider + elif used_DTD == "power-law-beta-one": + SNIa_rate_individual = SNIa_efficiency * power_law_beta_one_DTD(time_since_formation, delay_time, normalization_timescale) * dM_consider + elif used_DTD == "exponential": + SNIa_rate_individual = SNIa_efficiency * exponential_law_DTD(time_since_formation, exponential_decay, delay_time) * dM_consider + elif used_DTD == "Gaussian": + SNIa_rate_individual = Gaussian_SNIa_efficiency * Gaussian_law_DTD(time_since_formation, delay_time, Gaussian_characteristic_time, Gaussian_std_time) * dM_consider + + SNIa_rate_sum = np.sum(SNIa_rate_individual) + SNIa_rate[i] = SNIa_rate_sum + + SNIa_rate = (SNIa_rate / box_volume / unyt.Gyr).to(SNIa_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factors_use, SNIa_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicSNIaRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + capsize=1.0, + ) + ) + observation_labels.append(f"{obs_data.citation}") + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc=4, fontsize=4, ncol=2 +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(0, 1.6) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel(f"SNIa rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]") +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/SNIa_rate_history_based_on_CSFH.png") diff --git a/colibre/scripts/cosmic_log_CC_SN_rate.py b/colibre/scripts/cosmic_log_CC_SN_rate.py new file mode 100644 index 00000000..cd232a58 --- /dev/null +++ b/colibre/scripts/cosmic_log_CC_SN_rate.py @@ -0,0 +1,163 @@ +""" +Plots the cosmic log CCSN rate history. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +from swiftpipeline.argumentparser import ScriptArgumentParser + +arguments = ScriptArgumentParser( + description="Creates a log CC SN rate history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +CCSN_filenames = [f"{directory}/SNII.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.semilogx() + +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +CC_SN_rate_output_units = 1.0 / (unyt.yr * unyt.Mpc ** 3) + +for idx, (snapshot_filename, CCSN_filename, name) in enumerate( + zip(snapshot_filenames, CCSN_filenames, names) +): + data = np.loadtxt( + CCSN_filename, + usecols=(4, 6, 11), + dtype=[("a", np.float32), ("z", np.float32), ("CC SN rate", np.float32)], + ) + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + CC_SN_rate_units = 1.0 / (units.time * units.length ** 3) + + # a, Redshift, SFR + scale_factor = data["a"] + CC_SN_rate = (data["CC SN rate"] * CC_SN_rate_units).to(CC_SN_rate_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factor, CC_SN_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicCCSNRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + ) + ) + observation_labels.append(f"{obs_data.citation}") + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc="center right", fontsize="xx-small" +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") +ax.set_yscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(1e-1, 3e1) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel( + f"CC SN rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]" +) +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_CC_SN_rate_history.png") + diff --git a/colibre/scripts/cosmic_log_CC_SN_rate_CSFH_based.py b/colibre/scripts/cosmic_log_CC_SN_rate_CSFH_based.py new file mode 100644 index 00000000..e8b38fc9 --- /dev/null +++ b/colibre/scripts/cosmic_log_CC_SN_rate_CSFH_based.py @@ -0,0 +1,192 @@ +""" +Plots the cosmic core-collapse SNe rate based on the star formation history and use a logged y-axis. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from load_sfh_data import read_obs_data + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +sfr_output_units = unyt.msun / (unyt.year * unyt.Mpc ** 3) +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +SNIa_output_units = 1. / (unyt.year * unyt.Mpc ** 3) + +from swiftpipeline.argumentparser import ScriptArgumentParser + +def CC_SN_DTD(t, t_min, t_max): + mask = (t > t_min) & (t < t_max) + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1./(t_max - t_min) * (np.ones(len(t[mask]))) + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +arguments = ScriptArgumentParser( + description="Creates a CC SN history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +sfr_filenames = [f"{directory}/SFR.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.set_xscale("log") +ax.set_yscale("log") + + +for idx, (snapshot_filename, sfr_filename, name) in enumerate( + zip(snapshot_filenames, sfr_filenames, names) +): + data = np.genfromtxt(sfr_filename).T + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + boxsize = snapshot.metadata.boxsize + box_volume = boxsize[0] * boxsize[1] * boxsize[2] + + sfr_units = snapshot.gas.star_formation_rates.units + + # a, Redshift, SFR + scale_factor = data[2] + redshift = data[3] + star_formation_rate = (data[7] * sfr_units / box_volume).to(sfr_output_units) + times = data[1] * units.time + dM = data[4] * units.mass + + # Calculate the times we plot + times_desired = np.linspace(0.05,13.8,500) * unyt.Gyr + scale_factors_use = np.interp(times_desired, times, scale_factor) + + SNIa_rate = np.zeros(len(times_desired)) + + for i in range(1,len(times_desired)): + time_consider = times[times < times_desired[i]] + time_since_formation = times_desired[i] - time_consider + dM_consider = dM[times < times_desired[i]] + SNIa_rate_individual = 1.180e-2 * CC_SN_DTD(time_since_formation, 3*unyt.Myr, 0.04*unyt.Gyr) * dM_consider / unyt.Msun + + SNIa_rate_sum = np.sum(SNIa_rate_individual) + SNIa_rate[i] = SNIa_rate_sum + + SNIa_rate = (SNIa_rate / box_volume / unyt.Gyr).to(SNIa_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factors_use, SNIa_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicCCSNRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + capsize=1.0, + ) + ) + observation_labels.append(f"{obs_data.citation}") + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc=4, fontsize=4, ncol=2 +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(1e-1, 3e1) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel(f"CC SN rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]") +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_CC_SN_rate_history_based_on_CSFH.png") diff --git a/colibre/scripts/cosmic_log_CEJSN_rate.py b/colibre/scripts/cosmic_log_CEJSN_rate.py new file mode 100644 index 00000000..4263dbcf --- /dev/null +++ b/colibre/scripts/cosmic_log_CEJSN_rate.py @@ -0,0 +1,138 @@ +""" +Plots the cosmic CEJSN rate history based on the chemical enrichment logger. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +from swiftpipeline.argumentparser import ScriptArgumentParser + +arguments = ScriptArgumentParser( + description="Creates a CEJSN rate history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +r_process_filenames = [f"{directory}/r_processes.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.semilogx() + +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +r_process_rate_output_units = 1.0 / (unyt.yr * unyt.Mpc ** 3) + +for idx, (snapshot_filename, r_process_filename, name) in enumerate( + zip(snapshot_filenames, r_process_filenames, names) +): + data = np.loadtxt( + r_process_filename, + usecols=(2,3, 4, 6, 14, 15, 16), + dtype=[("t2", np.float32), ("t1", np.float32), ("a", np.float32), ("z", np.float32), ("NSM", np.float32), ("CEJSN", np.float32), ("collapsar", np.float32)], + ) + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + r_process_rate_units = 1.0 / (units.time * units.length ** 3) + + volume = snapshot.metadata.boxsize[0] * snapshot.metadata.boxsize[1]* snapshot.metadata.boxsize[2] + volume.convert_to_units("Mpc**3") + + dt = (data["t2"] - data["t1"]) + + # a, Redshift, SFR + scale_factor = data["a"] + NSM_rate = (data["NSM"]/volume.value /dt * r_process_rate_units).to(r_process_rate_output_units) + CEJSN_rate = (data["CEJSN"]/volume.value/dt * r_process_rate_units).to(r_process_rate_output_units) + collapsar_rate = (data["collapsar"]/volume.value/dt * r_process_rate_units).to(r_process_rate_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factor, CEJSN_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") +ax.set_yscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(1e-5, 1e-1) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel( + f"Common-envelop jets SN rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]" +) +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_CEJSN_rate_history.png") diff --git a/colibre/scripts/cosmic_log_NSM_rate.py b/colibre/scripts/cosmic_log_NSM_rate.py new file mode 100644 index 00000000..b33bcb3f --- /dev/null +++ b/colibre/scripts/cosmic_log_NSM_rate.py @@ -0,0 +1,138 @@ +""" +Plots the cosmic neutron star merger rate history based on the chemical enrichment logger. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +from swiftpipeline.argumentparser import ScriptArgumentParser + +arguments = ScriptArgumentParser( + description="Creates a neutral star merger rate history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +r_process_filenames = [f"{directory}/r_processes.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.semilogx() + +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +r_process_rate_output_units = 1.0 / (unyt.yr * unyt.Mpc ** 3) + +for idx, (snapshot_filename, r_process_filename, name) in enumerate( + zip(snapshot_filenames, r_process_filenames, names) +): + data = np.loadtxt( + r_process_filename, + usecols=(2,3, 4, 6, 14, 15, 16), + dtype=[("t2", np.float32), ("t1", np.float32), ("a", np.float32), ("z", np.float32), ("NSM", np.float32), ("CEJSN", np.float32), ("collapsar", np.float32)], + ) + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + r_process_rate_units = 1.0 / (units.time * units.length ** 3) + + volume = snapshot.metadata.boxsize[0] * snapshot.metadata.boxsize[1]* snapshot.metadata.boxsize[2] + volume.convert_to_units("Mpc**3") + + dt = (data["t2"] - data["t1"]) + + # a, Redshift, SFR + scale_factor = data["a"] + NSM_rate = (data["NSM"]/volume.value /dt * r_process_rate_units).to(r_process_rate_output_units) + CEJSN_rate = (data["CEJSN"]/volume.value/dt * r_process_rate_units).to(r_process_rate_output_units) + collapsar_rate = (data["collapsar"]/volume.value/dt * r_process_rate_units).to(r_process_rate_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factor, NSM_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") +ax.set_yscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(1e-5, 1e-1) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel( + f"NS-NS merger rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]" +) +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_NSM_rate_history.png") diff --git a/colibre/scripts/cosmic_log_SNIa_rate.py b/colibre/scripts/cosmic_log_SNIa_rate.py new file mode 100644 index 00000000..05317b48 --- /dev/null +++ b/colibre/scripts/cosmic_log_SNIa_rate.py @@ -0,0 +1,163 @@ +""" +Plots the cosmic log SNIa rate history. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +from swiftpipeline.argumentparser import ScriptArgumentParser + +arguments = ScriptArgumentParser( + description="Creates a log SNIa rate history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +SNIa_filenames = [f"{directory}/SNIa.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.semilogx() + +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +SNIa_rate_output_units = 1.0 / (unyt.yr * unyt.Mpc ** 3) + +for idx, (snapshot_filename, SNIa_filename, name) in enumerate( + zip(snapshot_filenames, SNIa_filenames, names) +): + data = np.loadtxt( + SNIa_filename, + usecols=(4, 6, 11), + dtype=[("a", np.float32), ("z", np.float32), ("SNIa rate", np.float32)], + ) + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + SNIa_rate_units = 1.0 / (units.time * units.length ** 3) + + # a, Redshift, SFR + scale_factor = data["a"] + SNIa_rate = (data["SNIa rate"] * SNIa_rate_units).to(SNIa_rate_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factor, SNIa_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicSNIaRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + capsize=1.0, + ) + ) + observation_labels.append(f"{obs_data.citation}") + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc="lower right", fontsize="xx-small", ncol=2 +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") +ax.set_yscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(3e-2, 2.0) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel( + f"SNIa rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]" +) +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_SNIa_rate_history.png") diff --git a/colibre/scripts/cosmic_log_SNIa_rate_CSFH_based.py b/colibre/scripts/cosmic_log_SNIa_rate_CSFH_based.py new file mode 100644 index 00000000..737e4f11 --- /dev/null +++ b/colibre/scripts/cosmic_log_SNIa_rate_CSFH_based.py @@ -0,0 +1,304 @@ +""" +Plots the cosmic SNIa rate history based on the star formation history and use a logged y-axis. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from load_sfh_data import read_obs_data + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +sfr_output_units = unyt.msun / (unyt.year * unyt.Mpc ** 3) +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +SNIa_output_units = 1. / (unyt.year * unyt.Mpc ** 3) + +from swiftpipeline.argumentparser import ScriptArgumentParser + +def power_law_beta_one_DTD(t, t_delay, tH): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = 1./(np.log(tH) - np.log(t_delay)) / t[mask] + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def power_law_DTD(t, t_delay, tH, beta): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = (1.-beta)/(tH**(1-beta) - (t_delay)**(1-beta)) / t[mask]**beta + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + + +def exponential_law_DTD(t, t_decay, t_delay): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + value_new = np.exp(-(t[mask]-t_delay)/t_decay)/t_decay + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def Gaussian_law_DTD(t, t_delay, tmean, tsigma): + mask = t > 0.04 + value = np.zeros(len(t)) / unyt.Gyr + norm = 1./(tsigma* np.sqrt(2*np.pi)) + delta_t = t[mask] - tmean + value_new = norm * np.exp(-0.5 * delta_t**2 / tsigma**2) + value_new.convert_to_units("Gyr**-1") + value[mask] = value_new + return value + +def broken_power_law_DTD(t, slope_short_time, slope_long_time, break_time_Gyr, delay_time, normalization_timescale): + mask1 = (t > 0.04) & (t < break_time_Gyr) + mask2 = (t>= break_time_Gyr) + value = np.zeros(len(t)) / unyt.Gyr + norm1 = (1- slope_short_time) * (1 - slope_long_time) + norm2a = (slope_short_time - slope_long_time) * break_time_Gyr + norm2b = (1 - slope_short_time) * normalization_timescale**(1-slope_long_time) * break_time_Gyr**slope_long_time + norm2c = (1 - slope_long_time) * normalization_timescale**(1-slope_short_time) * break_time_Gyr**slope_short_time + norm2 = norm2a + norm2b + norm2c + norm = norm1 / norm2 + value[mask1] = norm * (t[mask1]/break_time_Gyr)**(-slope_short_time) + value[mask2] = norm * (t[mask2]/break_time_Gyr)**(-slope_long_time) + return value + + +arguments = ScriptArgumentParser( + description="Creates a star formation history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +sfr_filenames = [f"{directory}/SFR.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.set_xscale("log") +ax.set_yscale("log") + + +for idx, (snapshot_filename, sfr_filename, name) in enumerate( + zip(snapshot_filenames, sfr_filenames, names) +): + data = np.genfromtxt(sfr_filename).T + + snapshot = load(snapshot_filename) + + # find out which DTD model we are currently using in the code + normalization_timescale = snapshot.metadata.parameters.get("SNIaDTD:normalization_timescale_Gyr", None) + if normalization_timescale is None: + have_normalization_timescale = False + else: + normalization_timescale = float(normalization_timescale) * unyt.Gyr + have_normalization_timescale = True + + exponential_decay = snapshot.metadata.parameters.get("SNIaDTD:SNIa_timescale_Gyr", None) + if exponential_decay is None: + have_exponential_decay = False + else: + have_exponential_decay = True + exponential_decay = float(exponential_decay) * unyt.Gyr + + Gaussian_SNIa_efficiency = snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", None) + if Gaussian_SNIa_efficiency is None: + have_Gaussian = False + else: + have_Gaussian = True + Gaussian_SNIa_efficiency = float(Gaussian_SNIa_efficiency) / unyt.Msun + + power_law_slope = snapshot.metadata.parameters.get("SNIaDTD:power_law_slope", None) + if power_law_slope is None: + have_slope = False + else: + have_slope = True + power_law_slope = float(power_law_slope) + + power_law_slope_short_time = snapshot.metadata.parameters.get("SNIaDTD:power_law_slope_short_time", None) + if power_law_slope_short_time is None: + have_broken_slope = False + else: + have_broken_slope = True + power_law_slope_short_time = float(power_law_slope_short_time) + + if have_Gaussian: + used_DTD = "Gaussian" + elif have_broken_slope: + used_DTD = "broken-power-law" + elif have_slope and have_normalization_timescale: + used_DTD = "power-law" + elif have_normalization_timescale: + used_DTD = "power-law-beta-one" + else: + used_DTD = "exponential" + + delay_time = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_delay_time_Gyr", False) ) * unyt.Gyr + + if used_DTD in ["power-law", "exponential", "power-law-beta-one"]: + SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_p_Msun", False) ) / unyt.Msun + elif used_DTD == "Gaussian": + Gaussian_SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_gauss_p_Msun", False)) / unyt.Msun + Gaussian_const_SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_const_p_Msun", False)) / unyt.Msun + Gaussian_characteristic_time = float(snapshot.metadata.parameters.get("SNIaDTD:characteristic_time_Gyr", False)) * unyt.Gyr + Gaussian_std_time = float(snapshot.metadata.parameters.get("SNIaDTD:STD_characteristic_time_Gyr", False)) * unyt.Gyr + elif used_DTD == "broken-power-law": + SNIa_efficiency = float(snapshot.metadata.parameters.get("SNIaDTD:SNIa_efficiency_p_Msun", False) ) / unyt.Msun + power_law_slope_long_time = float(snapshot.metadata.parameters.get("SNIaDTD:power_law_slope_long_time", False)) + break_time_Gyr = float(snapshot.metadata.parameters.get("SNIaDTD:break_time_Gyr", False) ) * unyt.Gyr + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + boxsize = snapshot.metadata.boxsize + box_volume = boxsize[0] * boxsize[1] * boxsize[2] + + sfr_units = snapshot.gas.star_formation_rates.units + + # a, Redshift, SFR + scale_factor = data[2] + redshift = data[3] + star_formation_rate = (data[7] * sfr_units / box_volume).to(sfr_output_units) + times = data[1] * units.time + dM = data[4] * units.mass + + # Calculate the times we plot + times_desired = np.linspace(0.05,13.8,500) * unyt.Gyr + scale_factors_use = np.interp(times_desired, times, scale_factor) + + SNIa_rate = np.zeros(len(times_desired)) + + for i in range(1,len(times_desired)): + time_consider = times[times < times_desired[i]] + time_since_formation = times_desired[i] - time_consider + dM_consider = dM[times < times_desired[i]] + if used_DTD == "power-law": + SNIa_rate_individual = SNIa_efficiency * power_law_DTD(time_since_formation, delay_time, normalization_timescale,power_law_slope) * dM_consider + elif used_DTD == "power-law-beta-one": + SNIa_rate_individual = SNIa_efficiency * power_law_beta_one_DTD(time_since_formation, delay_time, normalization_timescale) * dM_consider + elif used_DTD == "exponential": + SNIa_rate_individual = SNIa_efficiency * exponential_law_DTD(time_since_formation, exponential_decay, delay_time) * dM_consider + elif used_DTD == "Gaussian": + SNIa_rate_individual = Gaussian_SNIa_efficiency * Gaussian_law_DTD(time_since_formation, delay_time, Gaussian_characteristic_time, Gaussian_std_time) * dM_consider + elif used_DTD == "broken-power-law": + SNIa_rate_individual = SNIa_efficiency * broken_power_law_DTD(time_since_formation, power_law_slope_short_time, power_law_slope_long_time, break_time_Gyr, delay_time, normalization_timescale) * dM_consider + + SNIa_rate_sum = np.sum(SNIa_rate_individual) + SNIa_rate[i] = SNIa_rate_sum + + SNIa_rate = (SNIa_rate / box_volume / unyt.Gyr).to(SNIa_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factors_use, SNIa_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + +observation_lines = [] +observation_labels = [] + +path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}" +observational_data = load_observations( + sorted(glob.glob(f"{path_to_obs_data}/data/CosmicSNIaRate/*.hdf5")) +) + +for obs_data in observational_data: + observation_lines.append( + ax.errorbar( + obs_data.x.value, + obs_data.y.value * multiplicative_factor, + yerr=obs_data.y_scatter.value * multiplicative_factor, + label=obs_data.citation, + linestyle="none", + marker="o", + elinewidth=0.5, + markeredgecolor="none", + markersize=2, + zorder=-10, + capsize=1.0, + ) + ) + observation_labels.append(f"{obs_data.citation}") + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +observation_legend = ax.legend( + observation_lines, observation_labels, markerfirst=True, loc=4, fontsize=4, ncol=2 +) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +ax.add_artist(observation_legend) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(3e-2, 2.0) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel(f"SNIa rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]") +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_SNIa_rate_history_based_on_CSFH.png") diff --git a/colibre/scripts/cosmic_log_collapsar_rate.py b/colibre/scripts/cosmic_log_collapsar_rate.py new file mode 100644 index 00000000..7ffdfc2b --- /dev/null +++ b/colibre/scripts/cosmic_log_collapsar_rate.py @@ -0,0 +1,138 @@ +""" +Plots the cosmic collapser rate history based on the chemical enrichment logger. +""" +import unyt + +import matplotlib.pyplot as plt +import numpy as np +import glob + +from swiftsimio import load + +from velociraptor.observations import load_observations + +from astropy.cosmology import z_at_value +from astropy.units import Gyr + +from swiftpipeline.argumentparser import ScriptArgumentParser + +arguments = ScriptArgumentParser( + description="Creates a collapser rate history plot, with added observational data." +) + +snapshot_filenames = [ + f"{directory}/{snapshot}" + for directory, snapshot in zip(arguments.directory_list, arguments.snapshot_list) +] + +r_process_filenames = [f"{directory}/r_processes.txt" for directory in arguments.directory_list] + +names = arguments.name_list +output_path = arguments.output_directory + +plt.style.use(arguments.stylesheet_location) + +simulation_lines = [] +simulation_labels = [] + +fig, ax = plt.subplots() + +ax.semilogx() + +log_multiplicative_factor = 4 +multiplicative_factor = 10 ** log_multiplicative_factor +r_process_rate_output_units = 1.0 / (unyt.yr * unyt.Mpc ** 3) + +for idx, (snapshot_filename, r_process_filename, name) in enumerate( + zip(snapshot_filenames, r_process_filenames, names) +): + data = np.loadtxt( + r_process_filename, + usecols=(2,3, 4, 6, 14, 15, 16), + dtype=[("t2", np.float32), ("t1", np.float32), ("a", np.float32), ("z", np.float32), ("NSM", np.float32), ("CEJSN", np.float32), ("collapsar", np.float32)], + ) + + snapshot = load(snapshot_filename) + + # Read cosmology from the first run in the list + if idx == 0: + cosmology = snapshot.metadata.cosmology + + units = snapshot.units + r_process_rate_units = 1.0 / (units.time * units.length ** 3) + + volume = snapshot.metadata.boxsize[0] * snapshot.metadata.boxsize[1]* snapshot.metadata.boxsize[2] + volume.convert_to_units("Mpc**3") + + dt = (data["t2"] - data["t1"]) + + # a, Redshift, SFR + scale_factor = data["a"] + NSM_rate = (data["NSM"]/volume.value /dt * r_process_rate_units).to(r_process_rate_output_units) + CEJSN_rate = (data["CEJSN"]/volume.value/dt * r_process_rate_units).to(r_process_rate_output_units) + collapsar_rate = (data["collapsar"]/volume.value/dt * r_process_rate_units).to(r_process_rate_output_units) + + # High z-order as we always want these to be on top of the observations + simulation_lines.append( + ax.plot(scale_factor, collapsar_rate.value * multiplicative_factor, zorder=10000)[0] + ) + simulation_labels.append(name) + + +redshift_ticks = np.array([0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, 50.0, 100.0]) +redshift_labels = [ + "$0$", + "$0.2$", + "$0.5$", + "$1$", + "$2$", + "$3$", + "$5$", + "$10$", + "$20$", + "$50$", + "$100$", +] +a_ticks = 1.0 / (redshift_ticks + 1.0) + +ax.set_xticks(a_ticks) +ax.set_xticklabels(redshift_labels) + +simulation_legend = ax.legend( + simulation_lines, simulation_labels, markerfirst=False, loc="upper right" +) + +# Create second X-axis (to plot cosmic time alongside redshift) +ax2 = ax.twiny() +ax2.set_xscale("log") +ax.set_yscale("log") + +# Cosmic-time ticks (in Gyr) along the second X-axis +t_ticks = np.array([0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, cosmology.age(1.0e-5).value]) + +# To place the new ticks onto the X-axis we need to know the corresponding scale factors +a_ticks_2axis = [ + 1.0 / (1.0 + z_at_value(cosmology.age, t_tick * Gyr)) for t_tick in t_ticks +] + +# Attach the ticks to the second X-axis +ax2.set_xticks(a_ticks_2axis) + +# Format the ticks' labels +ax2.set_xticklabels(["$%2.1f$" % t_tick for t_tick in t_ticks]) + +# Final adjustments +ax.tick_params(axis="x", which="minor", bottom=False) +ax2.tick_params(axis="x", which="minor", top=False) + +ax.set_ylim(1e-5, 2.0) +ax.set_xlim(1.02, 0.07) +ax2.set_xlim(1.02, 0.07) + +ax.set_xlabel("Redshift $z$") +ax.set_ylabel( + f"Collapsar rate [$10^{{-{log_multiplicative_factor}}}$ yr$^{{-1}}$ cMpc$^{{-3}}$]" +) +ax2.set_xlabel("Cosmic time [Gyr]") + +fig.savefig(f"{output_path}/log_collapsar_rate_history.png") diff --git a/observational_data b/observational_data index a3fdc9fe..aec5c536 160000 --- a/observational_data +++ b/observational_data @@ -1 +1 @@ -Subproject commit a3fdc9feadf6f7ad9c024a4e51da197afde21230 +Subproject commit aec5c536f1ead5cb13257610246e5ce8d449ea4d