Update the time step calculating

src/runner.c

@@ -475,6 +475,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
   const struct hydro_props *restrict hydro_props = e->hydro_properties;
   const struct unit_system *restrict us = e->internal_units;
   struct cooling_function_data *restrict cooling = e->cooling_func;
+  const double time_base = e->time_base;
+  const integertime_t ti_current = e->ti_current;
 
   TIMER_TIC;
 
@@ -496,12 +498,27 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
 
       if (part_is_active(p, e)) {
 
+        /* Calculate the time step of the current particle*/
+        double dt_star;
+        if (with_cosmology) {
+          const integertime_t ti_step = get_integer_timestep(p->time_bin);
+          const integertime_t ti_begin =
+              get_integer_time_begin(ti_current - 1, p->time_bin);
+
+          dt_star =
+              cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
+
+        } else {
+          dt_star = get_timestep(p->time_bin, time_base);
+        }
+
         // const float rho = hydro_get_physical_density(p, cosmo);
         if (star_formation_convert_to_star(e, starform, p, xp, constants, cosmo,
-                                           hydro_props, us, cooling)) {
+                                           hydro_props, us, cooling, dt_star)) {
           /* Convert your particle to a star */
           struct spart* sp = cell_convert_part_to_spart(e, c, p, xp);
 
+          /* Copy the properties of the gas particle to the star particle */
           star_formation_copy_properties(e, p, xp, sp, starform, constants,
                                          cosmo, with_cosmology);
           // struct spart *sp =        cell_conert_part_to_spart(c, p, ...);

src/starformation/schaye08/starformation.h

@@ -156,7 +156,6 @@ INLINE static int star_formation_potential_to_become_star(
   if (particle_density < rho_crit_times_min_over_den)
     return 0;
 
-
   /* In this case there are actually multiple possibilities
    * because we also need to check if the physical density exceeded
    * the appropriate limit */
@@ -165,14 +164,14 @@ INLINE static int star_formation_potential_to_become_star(
   double density_threshold_metal_dep =
           starform->density_threshold * pow(Z * starform->Z0_inv, starform->n_Z0);
 
-  double density_threshold_current = min(density_threshold_metal_dep, starform->density_threshold_max);
+  /* Calculate the maximum between both and convert to mass density instead of number density*/
+  double density_threshold_current = min(density_threshold_metal_dep, starform->density_threshold_max) * phys_const->const_proton_mass;
   
   
   /* Check if it exceeded the maximum density */
-  if (particle_density < density_threshold_current)
+  if (particle_density*p->chemistry_data.smoothed_metal_mass_fraction[0] < density_threshold_current)
     return 0;
     
-      
   /* Calculate the temperature */
   const double temperature = cooling_get_temperature(phys_const, hydro_props, us, cosmo,
 						     cooling, p, xp);
@@ -200,18 +199,19 @@ INLINE static int star_formation_convert_to_star(
     const struct phys_const* const phys_const, const struct cosmology* cosmo,
     const struct hydro_props* restrict hydro_props,
     const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling) {
+    const struct cooling_function_data* restrict cooling,
+    const double dt_star) {
 
   if (star_formation_potential_to_become_star(
           starform, p, xp, phys_const, cosmo, hydro_props, us, cooling)) {
     /* Get the pressure */
     const double pressure =
         starform->EOS_pressure_norm *
-        pow(p->rho / starform->EOS_density_norm, starform->polytropic_index);
+        pow(hydro_get_physical_density(p,cosmo) / starform->EOS_density_norm, starform->polytropic_index);
 
     /* Calculate the propability of forming a star */
     const double prop = starform->SF_normalization *
-                        pow(pressure, starform->SF_power_law) * p->time_bin;
+                        pow(pressure, starform->SF_power_law) * dt_star;
 
     /* Calculate the seed */
     unsigned int seed = (p->id + e->ti_current) % 8191;
@@ -219,6 +219,8 @@ INLINE static int star_formation_convert_to_star(
     /* Generate a random number between 0 and 1. */
     const double randomnumber = rand_r(&seed) * starform->inv_RAND_MAX;
 
+    message("Passed whole boundary thing! random number = %e, prop = %e time_bin %d", randomnumber, prop, p->time_bin);
+
     /* Calculate if we form a star */
     return (prop > randomnumber);
 
@@ -304,7 +306,7 @@ INLINE static void starformation_init_backend(
 
   /* Conversion of number density from cgs */
   static const float dimension_numb_den[5] = {0, -3, 0, 0, 0};
-  const double conversion_numb_density =
+  const double conversion_numb_density = 1/
        units_general_cgs_conversion_factor(us, dimension_numb_den);
 
   /* Quantities that have to do with the Normal Kennicutt-