Adding ram pressure stripping of ISM

Included the option of stripping the ISM of galaxies due to ram
pressure. This can be computed with or without the halo stripping on,
but stripping as an option needs to be switched on.

There remains one issue to consider. Currently type2 satellites don't
experience any ram pressure stripping of the ISM after becoming type 2
(only type 1 experience ram pressure stripping). This is because it is
not clear how to compute a consistent satellite galaxy position and
velocity within the halo. That's to be refined in the future.

include/dark_matter_halos.h

@@ -93,7 +93,7 @@ class DarkMatterHalos {
 
 	double subhalo_dynamical_time (Subhalo &subhalo, double z);
 
-	double halo_virial_radius(HaloPtr &halo, double z);
+	double halo_virial_radius(const HaloPtr &halo, double z);
 
 	double halo_virial_velocity (double mvir, double redshift);
 
@@ -109,7 +109,7 @@ class DarkMatterHalos {
 
 	void cooling_gas_sAM(Subhalo &subhalo, double z);
 
-	float enclosed_total_mass(Subhalo &subhalo, double z, float r);
+	float enclosed_total_mass(const Subhalo &subhalo, double z, float r);
 
 	void disk_sAM(Subhalo &subhalo, Galaxy &galaxy, double z);
 

include/environment.h

@@ -49,6 +49,7 @@ class EnvironmentParameters{
 	bool stripping = true;
 	bool tidal_stripping = false;
 	float minimum_halo_mass_fraction = 0.01;
+	float alpha_rps_halo = 1;
 	float Accuracy_RPS = 0.05;
 
 };
@@ -64,12 +65,27 @@ class Environment{
 
 	void process_satellite_subhalo_environment (Subhalo &satellite_subhalo, SubhaloPtr &central_subhalo, double z);
 	BaryonBase remove_tidal_stripped_stars(SubhaloPtr &subhalo, Galaxy &galaxy, BaryonBase lost_stellar);
-	double process_ram_pressure_stripping(const SubhaloPtr &primary,
+
+	double process_ram_pressure_stripping_gas(const SubhaloPtr &primary,
 			Subhalo &secondary,
-			double z);
+			double z,
+			double ram_press,
+			bool halo_strip,
+			bool ism_strip);
+
 	double ram_pressure_stripping_hot_gas(const SubhaloPtr &primary,
-			Subhalo &secondary,
+			const Subhalo &secondary,
 			double r,
+			double z,
+			double ram_press);
+
+	double ram_pressure_stripping_galaxy_gas(const GalaxyPtr &galaxy,
+			double r,
+			double z,
+			double ram_press);
+
+	double ram_pressure(const SubhaloPtr &primary,
+			const Subhalo &secondary,
 			double z);
 
 	using func_x = double (*)(double x, void *);
@@ -83,6 +99,8 @@ class Environment{
 	SimulationParameters simparams;
 	Root_Solver root_solver;
 
+	float remove_gas(BaryonBase &component, double m_removed, float f_gas);
+
 };
 
 using EnvironmentPtr = std::shared_ptr<Environment>;

include/galaxy.h

@@ -32,6 +32,7 @@
 #include "baryon.h"
 #include "components.h"
 #include "mixins.h"
+#include "numerical_constants.h"
 
 namespace shark {
 
@@ -129,6 +130,10 @@ class Galaxy : public Identifiable<galaxy_id_t> {
 	/// maximum circular velocity.
 	float vmax = 0;
 
+	// ram-pressure stripping radius
+	float r_rps = 0;
+	BaryonBase ram_pressure_stripped_gas;
+
 	/// star formation and gas history of this galaxy across snapshots
 	std::vector<HistoryItem> history;
 
@@ -295,15 +300,21 @@ class Galaxy : public Identifiable<galaxy_id_t> {
 
 	}
 
+	double enclosed_mass_gas(float r) const
+	{
+
+		return enclosed_mass_exponential(r, disk_gas.mass, disk_gas.rscale) +
+				enclosed_mass_exponential(r, bulge_gas.mass, bulge_gas.rscale);
+	}
+
 	double enclosed_mass_exponential(float r, float m, float r50) const
 	{
 		if(m == 0){
 			return 0;
 		}
 		else{
 			auto rnorm = r/(r50/1.67);
-			auto mass = m * (1 - (1 + rnorm) * std::exp(-rnorm));
-			return mass;
+			return m * (1 - (1 + rnorm) * std::exp(-rnorm));
 		}
 	}
 
@@ -318,6 +329,49 @@ class Galaxy : public Identifiable<galaxy_id_t> {
 		}
 	}
 
+	double surface_density_disk(float r){
+		return surface_density_exponential(r, disk_gas.mass, disk_gas.rscale) +
+				surface_density_exponential(r, disk_stars.mass, disk_stars.rscale);
+	}
+
+	double surface_density_bulge(float r){
+		return surface_density_exponential(r, bulge_gas.mass, bulge_gas.rscale) +
+				surface_density_plummer(r, bulge_stars.mass, bulge_stars.rscale);
+	}
+
+	double surface_density_gas(float r){
+		return surface_density_exponential(r, disk_gas.mass, disk_gas.rscale) +
+				surface_density_exponential(r, bulge_gas.mass, bulge_gas.rscale);
+	}
+
+	double surface_density_stars(float r){
+		return surface_density_exponential(r, disk_stars.mass, disk_stars.rscale) +
+				surface_density_plummer(r, bulge_stars.mass, bulge_stars.rscale);
+	}
+
+	double surface_density_exponential(float r, float m, float r50) const
+	{
+		if(m == 0){
+			return 0;
+		}
+		else{
+			auto re = r50/1.67;
+			auto rnorm = r/re;
+			return m / (constants::PI2 * std::pow(re,2)) * std::exp(-rnorm);
+		}
+	}
+
+	double surface_density_plummer(float r, float m, float r50) const
+	{
+		if(m ==0){
+			return 0;
+		}
+		else{
+			auto re = r50/1.3;
+			return m * std::pow(re, 2) / (constants::PI * std::pow( std::pow(r,2) + std::pow(re,2), 2));
+		}
+	}
+
 };
 
 // Support for less-based comparison of Galaxy objects. This allows them to be

include/mixins.h

@@ -100,6 +100,38 @@ class xyz {
 		return xyz(sum);
 	}
 
+	xyz &operator-=(T scalar)
+	{
+		x -= scalar;
+		y -= scalar;
+		z -= scalar;
+		return *this;
+	}
+
+	xyz operator-(T scalar) const
+	{
+		xyz sum(*this);
+		sum -= scalar;
+		return sum;
+	}
+
+	template <typename U>
+	xyz &operator-=(const xyz<U> &lhs)
+	{
+		x -= lhs.x;
+		y -= lhs.y;
+		z -= lhs.z;
+		return *this;
+	}
+
+	template <typename U>
+	xyz operator-(const xyz<U> &lhs) const
+	{
+		xyz<U> sum(*this);
+		sum -= lhs;
+		return xyz(sum);
+	}
+
 	xyz &operator*=(T scalar)
 	{
 		x *= scalar;

src/dark_matter_halos.cpp

@@ -155,7 +155,7 @@ double DarkMatterHalos::subhalo_dynamical_time (Subhalo &subhalo, double z){
 	return constants::MPCKM2GYR * cosmology->comoving_to_physical_size(r, z) / v;
 }
 
-double DarkMatterHalos::halo_virial_radius(HaloPtr &halo, double z){
+double DarkMatterHalos::halo_virial_radius(const HaloPtr &halo, double z){
 
 	/**
 	 * Function to calculate the halo virial radius. Returns virial radius in physical Mpc/h.
@@ -293,9 +293,9 @@ void DarkMatterHalos::cooling_gas_sAM(Subhalo &subhalo, double z){
 
 }
 
-float DarkMatterHalos::enclosed_total_mass(Subhalo &subhalo, double z, float r){
+float DarkMatterHalos::enclosed_total_mass(const Subhalo &subhalo, double z, float r){
 
-	GalaxyPtr galaxy;
+	ConstGalaxyPtr galaxy;
 	if(subhalo.subhalo_type == Subhalo::SATELLITE){
 		galaxy = subhalo.type1_galaxy();
 	}