fix units in wind mass loading calculation

examples/parameter_example.yml

@@ -541,7 +541,7 @@ SIMBAFeedback:
   wind_scatter:                     0.1
   delay_time:                       0.05
   wind_mass_eta:                    9.
-  wind_mass_spectrum_break:         5.2e9
+  wind_mass_spectrum_break_msun:    5.2e9
   low_mass_power:                   -0.317
   high_mass_power:                  -0.761
   sn_energy_erg:                    1.e51

src/feedback/SIMBA/feedback.h

@@ -65,12 +65,12 @@ inline void compute_kick_speed(struct spart *sp, const struct feedback_props *fe
  */
 inline void compute_mass_loading(struct spart *sp, const struct feedback_props *feedback_props) {
   // ALEXEI: check units of host_galaxy_mass with Romeel: is this a mass or mass/msun?
-  if (sp->mass < feedback_props->simba_mass_spectrum_break) {
+  if (sp->mass < feedback_props->simba_mass_spectrum_break_msun) {
     sp->feedback_data.to_distribute.wind_mass = feedback_props->simba_wind_mass_eta 
-        * sp->feedback_data.host_galaxy_mass * pow(sp->mass/feedback_props->simba_mass_spectrum_break,feedback_props->simba_low_mass_power);
+        * sp->feedback_data.host_galaxy_mass * pow(sp->mass/feedback_props->simba_mass_spectrum_break_msun,feedback_props->simba_low_mass_power);
   } else {
     sp->feedback_data.to_distribute.wind_mass = feedback_props->simba_wind_mass_eta 
-        * sp->feedback_data.host_galaxy_mass * pow(sp->mass/feedback_props->simba_mass_spectrum_break,feedback_props->simba_high_mass_power);
+        * sp->feedback_data.host_galaxy_mass * pow(sp->mass/feedback_props->simba_mass_spectrum_break_msun,feedback_props->simba_high_mass_power);
   }
 };
 
@@ -101,11 +101,12 @@ inline void compute_heating(struct spart *sp, const struct feedback_props *feedb
     u_SN *= 2.61634;
   }
 
-  // ALEXEI: commented out for debugging
-  //if (u_wind > feedback_props->simba_wind_energy_limit) 
-  //  sp->feedback_data.to_distribute.v_kick *= sqrt(feedback_props->simba_wind_energy_limit*u_SN/u_wind);
-  //if (feedback_props->simba_wind_energy_limit < 1.f) 
-  //  u_SN *= feedback_props->simba_wind_energy_limit;
+  if (u_wind > feedback_props->simba_wind_energy_limit) {
+    sp->feedback_data.to_distribute.v_kick *= sqrt(feedback_props->simba_wind_energy_limit*u_SN/u_wind);
+  }
+  if (feedback_props->simba_wind_energy_limit < 1.f) {
+    u_SN *= feedback_props->simba_wind_energy_limit;
+  }
 
   /* Now we can decide if there's any energy left over to distribute */
   sp->feedback_data.to_distribute.u_extra = max(u_SN - u_wind, 0.);

src/feedback/SIMBA/feedback_properties.h

@@ -42,7 +42,7 @@ struct feedback_props {
 
   /* mass loading parameters */
   float simba_wind_mass_eta;
-  float simba_mass_spectrum_break;
+  float simba_mass_spectrum_break_msun;
   float simba_low_mass_power;
   float simba_high_mass_power;
 
@@ -73,6 +73,8 @@ INLINE static void feedback_props_init(struct feedback_props *fp,
                                        const struct cosmology *cosmo) {
   // ALEXEI: double check units, make sure all use one system!!!
 
+  const double msun = 1.989e33 / units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+
   /* Initialize parameters for calculating rotational velocity of galaxy */
   fp->simba_host_galaxy_mass_norm = parser_get_param_float(params, "SIMBAFeedback:galaxy_mass_norm") / units_cgs_conversion_factor(us, UNIT_CONV_MASS); // 102.329 ALEXEI: guide values added in until figured out what are appropriate values.
   fp->simba_v_circ_exp = parser_get_param_float(params, "SIMBAFeedback:v_circ_exp"); // 0.26178;
@@ -88,7 +90,7 @@ INLINE static void feedback_props_init(struct feedback_props *fp,
 
   /* read parameters for mass loading calculation */
   fp->simba_wind_mass_eta = parser_get_param_float(params,"SIMBAFeedback:wind_mass_eta");
-  fp->simba_mass_spectrum_break = parser_get_param_float(params,"SIMBAFeedback:wind_mass_spectrum_break") / units_cgs_conversion_factor(us,UNIT_CONV_MASS);
+  fp->simba_mass_spectrum_break_msun = parser_get_param_float(params,"SIMBAFeedback:wind_mass_spectrum_break_msun") / msun;
   fp->simba_low_mass_power = parser_get_param_float(params,"SIMBAFeedback:low_mass_power");
   fp->simba_high_mass_power = parser_get_param_float(params,"SIMBAFeedback:high_mass_power");