Added new black spin model

include/agn_feedback.h

@@ -25,10 +25,13 @@
 #define INCLUDE_AGN_FEEDBACK_H_
 
 #include <memory>
+#include <random>
 #include <utility>
 
 #include "components.h"
 #include "cosmology.h"
+#include "execution.h"
+#include "galaxy.h"
 #include "options.h"
 #include "recycling.h"
 #include "subhalo.h"
@@ -50,22 +53,43 @@ class AGNFeedbackParameters {
 	double tau_fold = 0;
 	double accretion_eff_cooling = 0;
 	double kappa_agn = 0;
-	double nu_smbh = 0;
+	double nu_smbh = 0.1;
 
 	bool qso_feedback = false;
-	bool spin_v07 = false;
-	double kappa_radio = 0;
 	double epsilon_qso = 0;
-	double eta_superedd = 4;
-	double hot_halo_threshold = 0;
+
+	// Parameters relevant to Bravo22 model
+	float kappa_radio = 1;
+	float eta_superedd = 4;
+	float hot_halo_threshold = 1;
+	float alpha_adaf = 0.1;
+	float alpha_td = 0.1;
+	float delta_adaf = 0.2;
+	float mdotcrit_adaf = 0.01;
+	float low_accretion_adaf = 1e-4;
+	float constant_lowlum_adaf = 0;
+	float constant_highlum_adaf = 0;
 
 	enum AGNFeedbackModel {
 		CROTON16 = 0,
 		BOWER06,
-		BRAVO19
+		BRAVO22
+	};
+
+	enum SpinModel{
+		VOLONTERI07 = 0,
+		GRIFFIN20,
+		CONSTANTSPIN
+	};
+
+	enum AccretionDiskModel{
+		WARPEDDISK = 0,
+		SELFGRAVITYDISK
 	};
 
 	AGNFeedbackModel model = CROTON16;
+	SpinModel spin_model = VOLONTERI07;
+	AccretionDiskModel accretion_disk_model;
 };
 
 
@@ -80,26 +104,36 @@ class AGNFeedback {
 
 	void plant_seed_smbh(Subhalo &subhalo);
 	double eddington_luminosity(double mbh);
-	double accretion_rate_hothalo_smbh(double Lcool, double mbh);
+	double accretion_rate_hothalo_smbh(double Lcool, double tacc, BlackHole &smbh);
+	double accretion_rate_hothalo_smbh_limit(double mheatrate, double vvir, BlackHole &smbh);
 	double accretion_rate_ratio(double macc, double mBH);
-	double agn_bolometric_luminosity(double macc, double mBH);
-	double agn_mechanical_luminosity(double macc, double mBH);
-	double smbh_growth_starburst(double mgas, double vvir);
+	double agn_bolometric_luminosity(BlackHole &smbh);
+	double agn_mechanical_luminosity(BlackHole &smbh);
+	double smbh_growth_starburst(double mgas, double vvir, double tacc, BlackHole &smbh);
 	double smbh_accretion_timescale(Galaxy &galaxy, double z);
-	double accretion_rate_hothalo_smbh_limit(double mheatrate, double vvir);
 	double qso_critical_luminosity(double mgas, double m, double r);
 	double salpeter_timescale(double Lbol, double mbh);
 	double qso_outflow_velocity(double Lbol, double mbh, double zgas, double mgas, double mbulge, double rbulge);
-	void qso_outflow_rate(double mgas, double macc, double mBH, double zgas, double vcirc,
+	void qso_outflow_rate(double mgas, BlackHole &smbh, double zgas, double vcirc,
 			double sfr, double mbulge, double rbulge, double &beta_halo, double &beta_ejec);
-
+	void griffin20_spinup_accretion(double delta_mbh, double tau_acc, BlackHole &smbh);
+	void griffin20_spinup_mergers(BlackHole &smbh_primary, const BlackHole &smbh_secondary);
+	void volonteri07_spin(BlackHole &smbh);
+	float efficiency_luminosity_agn(float spin, float r_lso);
 
 	// TODO: move this to private when possible
 	AGNFeedbackParameters parameters;
 
 private:
 	CosmologyPtr cosmology;
 	RecyclingParameters recycle_params;
+	ExecutionParameters exec_params;
+
+	std::uniform_real_distribution<double> distribution;
+
+	double angle_acc_disk(const BlackHole &smbh);
+	double phi_acc_disk(const BlackHole &smbh);
+	double final_spin(const double mbh, const double mfin, const double r_lso);
 
 };
 

include/baryon.h

@@ -160,6 +160,10 @@ class BlackHole : public BaryonBase {
 
 	/** macc_sb: accretion rate onto the black hole during starbursts. */
 	float macc_sb = 0;
+
+	/** spin: black hole spin. */
+	float spin = 0;
+
 };
 
 }  // namespace shark

include/components.h

@@ -29,6 +29,8 @@
 #include <cstdint>
 #include <memory>
 
+#include "galaxy.h"
+
 namespace shark {
 
 using galaxy_id_t = int;

include/galaxy_mergers.h

@@ -88,6 +88,7 @@ class GalaxyMergers {
 			CosmologyPtr cosmology,
 			CosmologicalParameters cosmo_params,
 			ExecutionParameters execparams,
+			AGNFeedbackParameters agn_params,
 			SimulationParameters simparams,
 			DarkMatterHalosPtr darkmatterhalo,
 			std::shared_ptr<BasicPhysicalModel> physicalmodel,
@@ -155,6 +156,7 @@ class GalaxyMergers {
 	std::shared_ptr<Cosmology> cosmology;
 	CosmologicalParameters cosmo_params;
 	ExecutionParameters exec_params;
+	AGNFeedbackParameters agn_params;
 	SimulationParameters simparams;
 	DarkMatterHalosPtr darkmatterhalo;
 	std::shared_ptr<BasicPhysicalModel> physicalmodel;

include/numerical_constants.h

@@ -52,6 +52,9 @@ namespace constants {
 	/** Square root of 2 */
 	constexpr double SQRT2 = 1.4142135623730951;
 
+	/** Square root of 3 */
+	constexpr double SQRT3 = 1.7320508075688772;
+
 	/** pi */
 	constexpr double PI = 3.14159265358979323846;
 

include/physical_model.h

@@ -63,8 +63,7 @@ class PhysicalModel {
 	 * redshift: current redshift.
 	 * vsubh: subhalo virial velocity.
 	 * vgal: galaxy velocity at the half-mass radius of disk or bulge.
-	 * mBHacc: BH accretion rate in the case of starbursts.
-	 * mBH: supermassive black hole mass.
+	 * smbh: supermassive black hole mass.
 	 * burst: whether this is a starburst or not.
 	 */
 	struct solver_params {
@@ -79,16 +78,15 @@ class PhysicalModel {
 		double redshift;
 		double vsubh;
 		double vgal;
-		double mBHacc;
-		double mBH;
+		BlackHole &smbh;
 	};
 
 	PhysicalModel(
 			double ode_solver_precision,
 			ODESolver::ode_evaluator evaluator,
 			GasCooling gas_cooling) :
-		params {*this, false, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
-		starburst_params {*this, true, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
+		params {*this, false, 0., 0., 0., 0., 0., 0., 0., 0., 0., nullptr},
+		starburst_params {*this, true, 0., 0., 0., 0., 0., 0., 0., 0., 0., nullptr},
 		ode_solver(evaluator, NC, ode_solver_precision, &params),
 		starburst_ode_solver(evaluator, NC, ode_solver_precision, &starburst_params),
 		ode_values(NC), starburst_ode_values(NC),
@@ -112,8 +110,7 @@ class PhysicalModel {
 		 * rstar: half-stellar mass radius of the disk [Mpc/h]
 		 * vsubh: virial velocity of the host subhalo [km/s]
 		 * jcold_halo: specific angular momentum of the cooling gas [Msun/h Mpc/h km/s]
-		 * mBHacc: BH accretion rate due to starbursts; \equiv 0 in the case of star formation in disks.
-		 * mBH: supermassive black hole mass.
+		 * smbh: supermassive black hole class.
 		 * burst: boolean parameter indicating if this is a starburst or not.
 		 */
 
@@ -140,7 +137,7 @@ class PhysicalModel {
 		params.vsubh      = subhalo.Vvir;
 		params.jcold_halo = subhalo.cold_halo_gas.sAM;
 		params.delta_t = delta_t;
-		params.mBH = galaxy.smbh.mass;
+		params.smbh = galaxy.smbh;
 		params.redshift = z;
 
 		from_galaxy(ode_values, subhalo, galaxy);
@@ -161,19 +158,17 @@ class PhysicalModel {
 		 * rstar: half-stellar mass radius of the bulge [Mpc/h]
 		 * vsubh: virial velocity of the host subhalo [km/s]
 		 * jcold_halo: specific angular momentum of the cooling gas [Msun/h Mpc/h km/s]
-		 * mBHacc: BH accretion rate in [Msun/Gyr/h]
-		 * mBH: supermassive black hole mass [Msun/h]
+		 * smbh: supermassive black hole class.
 		 * burst: boolean parameter indicating if this is a starburst or not.
 		 */
 
 		starburst_params.rgas = galaxy.bulge_gas.rscale; //gas scale radius.
 		starburst_params.rstar = galaxy.bulge_stars.rscale; //stellar scale radius.
 		starburst_params.vsubh = subhalo.Vvir;
 		starburst_params.vgal = galaxy.bulge_gas.sAM / galaxy.bulge_gas.rscale;
-		starburst_params.mBHacc = galaxy.smbh.macc_sb;
-		starburst_params.mBH = galaxy.smbh.mass;
 		starburst_params.delta_t = delta_t;
 		starburst_params.redshift = z;
+		starburst_params.smbh = galaxy.smbh;
 
 		from_galaxy_starburst(starburst_ode_values, subhalo, galaxy);
 		starburst_ode_solver.evolve(starburst_ode_values, delta_t);