Cleanup and added a new condition for MFOFvir

Clean up from all the useless comments.

Also added in read_halos-velociraptor a new condition to read the FOFMvir instead of M200crit when the last one is not properly defined

src/core/BrightnessTemperature.c

@@ -16,10 +16,13 @@
  * a prescription for line-of-sight redshift-space distortions, taken from
  * 21CMMC (Greig & Mesinger 2018).
  * Written by Bradley Greig.
+ *
+ * Added the computation of 21cm Tb without Pop. III (useful to better
+ * highlight the contribution of Pop III stars). By Manu Ventura.
  */
 
-void ComputeBrightnessTemperatureBox(int snapshot) // Added the computation of 21cm Tb without Pop. III (useful to
-                                                   // better highlight the contribution of Pop III stars)
+void ComputeBrightnessTemperatureBox(int snapshot) // 
+                                                   // 
 {
 
   int ii, jj, kk, i_real, i_padded, iii;

src/core/ComputePowerSpectrum.c

@@ -41,7 +41,7 @@ void Initialise_PowerSpectrum()
   mlog("Initialise_PowerSpectrum set PS_Length to %d.", MLOG_MESG, run_globals.params.PS_Length);
 }
 
-void Compute_PS(int snapshot) // Adding the 21cm PS if only Pop II are present! Still not 100% sure. Still not
+void Compute_PS(int snapshot) // Adding the 21cm PS if only Pop II are present! Still to be tested. Still not
                               // disentangling reionization.
 {
 

src/core/ComputeTs.c

@@ -20,7 +20,7 @@
  */
 
 /*
- * The minihalo feature was written by Emanuele M. Ventura, which includes an
+ * The minihalo feature was written by Manu Ventura, which includes an
  * amalgamation of the requisite functions for LW background using formulation of Qin2020a
  *
  * The output with II are computed accounting only for Pop II galaxies. In particular we compute Lyman-alpha, Xray, LW
@@ -552,8 +552,8 @@ void _ComputeTs(int snapshot)
 
     Luminosity_converstion_factor_GAL *= (SEC_PER_YEAR) / (PLANCK);
 
-    // Do the same for Pop III (Once again, there are a few things you are not 100% sure)
-
+    // Do the same for Pop III. 
+    
 #if USE_MINI_HALOS
     if (fabs(run_globals.params.physics.SpecIndexXrayIII - 1.0) < 0.000001) {
       Luminosity_converstion_factor_III =
@@ -907,4 +907,4 @@ void ComputeTs(int snapshot, timer_info* timer_total)
        snapshot,
        timer_gpu,
        timer_delta(*timer_total));
-}
+}

src/core/PopIII.c

@@ -80,8 +80,6 @@ void initialize_PopIII() // Initialize PopIII quantities that are easily compute
       break;
   }
 
-  // double MminIMF = run_globals.params.physics.MminIMF;
-  // double MmaxIMF = run_globals.params.physics.MmaxIMF;
   assert(MminIMF < MmaxIMF);
 
   run_globals.params.physics.MminIMF = MminIMF;
@@ -576,4 +574,4 @@ double NLBias(double Dist_Radius, double Halo_Mass, double redshift)
 
   return (Psi_Norm * pow(Dist_Radius / 0.01, Alpha_ind) * pow(Halo_Mass / 1e6, Beta_ind) *
           pow(redshift / 20.0, Gamma_ind));
-}
+}

src/core/dracarys.c

@@ -15,7 +15,7 @@
 #include "metal_evo.h"
 #include "read_grids.h"
 #include "stellar_feedback.h"
-#include "virial_properties.h" //For the test
+#include "virial_properties.h" 
 #endif
 #include "save.h"
 #include "tree_flags.h"

src/core/find_HII_bubbles.c

@@ -19,6 +19,8 @@
  *
  * Inclusion of electron fraction (X-ray heating) and inhomogeneous recombinations
  * by Bradley Greig. Relevant functions taken from public version of 21cmFAST.
+ *
+ * Inclusion od Pop. III galaxies when accounting for MINIHALOS by Manu Ventura.
  */
 
 double RtoM(double R)
@@ -529,4 +531,4 @@ void find_HII_bubbles(int snapshot, timer_info* timer_total)
        snapshot,
        timer_gpu,
        timer_delta(*timer_total));
-}
+}

src/core/init.c

@@ -283,7 +283,7 @@ void init_meraxes()
   run_globals.RequestedMassRatioModifier = 1;
   run_globals.RequestedBaryonFracModifier = 1;
 
-  // read in the mean Mvir_crit table (if needed 1 for Reio 2 for LW)
+  // read in the mean Mvir_crit table (if needed, 1 for Reio 2 for LW)
   read_Mcrit_table(1);
 #ifdef USE_MINI_HALOS
   read_Mcrit_table(2);
@@ -316,4 +316,4 @@ void init_meraxes()
 
   // This will be set by Mhysa
   run_globals.mhysa_self = NULL;
-}
+}

src/core/metal_evo.c

@@ -141,7 +141,7 @@ void construct_metal_grids(int snapshot, int local_ngals)
             case prop_prob:
 
               if (gal->RmetalBubble >= 3 * gal->Rvir) { // A bubble can actually pollute the IGM only if it's bigger
-                                                        // than its virial radius (Try with 3 atm)
+                                                        // than its virial radius.
                 buffer_metals[ind] +=
                   (4.0 / 3.0 * M_PI *
                    pow((gal->RmetalBubble) * (1 + redshift), 3.0)); // cMpc/h (same units of cell volume)
@@ -550,8 +550,7 @@ int map_galaxies_to_slabs_metals(int ngals)
 
 void assign_probability_to_galaxies(int ngals_in_metal_slabs,
                                     int snapshot,
-                                    int flag_property) // Right now you are assigning a probability to all galaxies!
-                                                       // Actually you need only to the newly formed
+                                    int flag_property) 
 {
   // Same way in which we assing Mcrit due to Reio and LW feedback in reionization.c
 
@@ -844,10 +843,8 @@ void assign_probability_to_galaxies(int ngals_in_metal_slabs,
           gal->Metals_IGM = (double)buffer_metals[grid_index(ix, iy, iz, MetalGridDim, INDEX_REAL)];
 
         if (flag_property == 2) {
-          // gal->Gas_IGM = (double)buffer_metals[grid_index(ix, iy, iz, MetalGridDim, INDEX_REAL)] + cell_gas;
           gal->Gas_IGM = (double)buffer_metals[grid_index(ix, iy, iz, MetalGridDim, INDEX_REAL)];
-          gal->Metallicity_IGM = calc_metallicity(gal->Gas_IGM, gal->Metals_IGM); // Maybe this can be put somewhere
-                                                                                  // else
+          gal->Metallicity_IGM = calc_metallicity(gal->Gas_IGM, gal->Metals_IGM); 
         }
 
         if (flag_property == 3)
@@ -869,4 +866,4 @@ void assign_probability_to_galaxies(int ngals_in_metal_slabs,
   mlog("...done.", MLOG_CLOSE);
 }
 
-#endif
+#endif

src/core/read_grids.c

@@ -56,7 +56,7 @@ double calc_resample_factor(int n_cell[3])
 
 #if USE_MINI_HALOS
 void smooth_Densitygrid_real(
-  int snapshot) // Need this to put the overdensity in the metal grid (added by Manu, we are working in real space)
+  int snapshot) // Need this to put the overdensity in the metal grid.
 {
   mlog("Smoothing the overdensity of the reionization grid into the metal grid...", MLOG_MESG);
   reion_grids_t* reiogrids = &(run_globals.reion_grids);
@@ -354,4 +354,4 @@ void free_grids_cache()
     free(snapshot_vel);
     free(snapshot_deltax);
   }
-}
+}

src/core/read_halos-velociraptor.c

@@ -319,8 +319,16 @@ void read_trees__velociraptor(int snapshot,
             // fof_group->Mvir = tree_entry.Mass_FOF;
             // fof_group->Rvir = -1;
           } else {
-            fof_group->Mvir = tree_entry.Mass_200crit;
-            fof_group->Rvir = tree_entry.R_200crit;
+              if ((tree_entry.Mass_200crit >= tree_entry.Mass_FOF) || (tree_entry.Mass_200crit < tree_entry.Mass_tot)) {
+                // Adding this since we found an issue in the N-body
+                fof_group->Mvir = tree_entry.Mass_FOF;
+                fof_group->Rvir = -1;
+              }
+
+              else {
+                fof_group->Mvir = tree_entry.Mass_200crit;
+                fof_group->Rvir = tree_entry.R_200crit;
+              }
           }
           fof_group->Vvir = -1;
           fof_group->FOFMvirModifier = 1.0;

src/core/read_params.c

@@ -802,15 +802,6 @@ void read_parameter_file(char* fname, int mode)
       required_tag[n_param] = 1;
       params_type[n_param++] = PARAM_TYPE_DOUBLE;
 
-      /*strncpy(params_tag[n_param], "ReionNionPhotPerBary_III", tag_length);
-      params_addr[n_param] = &(run_params->physics).ReionNionPhotPerBaryIII;
-#if USE_MINI_HALOS
-      required_tag[n_param] = 1;
-#else
-      required_tag[n_param] = 0;
-#endif
-      params_type[n_param++] = PARAM_TYPE_DOUBLE;*/
-
       strncpy(params_tag[n_param], "BlackHoleMassLimitReion", tag_length);
       params_addr[n_param] = &(run_params->physics).BlackHoleMassLimitReion;
       required_tag[n_param] = 1;
@@ -1154,33 +1145,6 @@ void read_parameter_file(char* fname, int mode)
 #endif
       params_type[n_param++] = PARAM_TYPE_DOUBLE;
 
-      /*    strcpy(params_tag[n_param], "AlphaIMF");
-            params_addr[n_param] = &(run_params->physics).AlphaIMF;
-      #if USE_MINI_HALOS
-            required_tag[n_param] = 1;
-      #else
-            required_tag[n_param] = 0;
-      #endif
-            params_type[n_param++] = PARAM_TYPE_DOUBLE;
-
-            strcpy(params_tag[n_param], "MminIMF");
-            params_addr[n_param] = &(run_params->physics).MminIMF;
-      #if USE_MINI_HALOS
-            required_tag[n_param] = 1;
-      #else
-            required_tag[n_param] = 0;
-      #endif
-            params_type[n_param++] = PARAM_TYPE_DOUBLE;
-
-            strcpy(params_tag[n_param], "MmaxIMF");
-            params_addr[n_param] = &(run_params->physics).MmaxIMF;
-      #if USE_MINI_HALOS
-            required_tag[n_param] = 1;
-      #else
-            required_tag[n_param] = 0;
-      #endif
-            params_type[n_param++] = PARAM_TYPE_DOUBLE;*/
-
       strncpy(params_tag[n_param], "PopIII_IMF", tag_length);
       params_addr[n_param] = &(run_params->physics).PopIII_IMF;
 #if USE_MINI_HALOS
@@ -1270,4 +1234,4 @@ void read_parameter_file(char* fname, int mode)
   // If running mpi then broadcast the run parameters to all cores
   MPI_Bcast(run_params, sizeof(run_params_t), MPI_BYTE, 0, run_globals.mpi_comm);
   MPI_Bcast(&(run_globals.units), sizeof(run_units_t), MPI_BYTE, 0, run_globals.mpi_comm);
-}
+}

src/core/reionization.c

@@ -1150,7 +1150,8 @@ int map_galaxies_to_slabs(int ngals)
   return gal_counter;
 }
 
-void assign_Mvir_crit_to_galaxies(int ngals_in_slabs, int flag_feed) // flag = 1 Reio feedback, flag = 2 LW feedback
+void assign_Mvir_crit_to_galaxies(int ngals_in_slabs, int flag_feed) 
+// flag = 1 Reio feedback, flag = 2 LW feedback
 {
   // N.B. We are assuming here that the galaxy_to_slab mapping has been sorted
   // by slab index...
@@ -2133,4 +2134,4 @@ void velocity_gradient(fftwf_complex* box, int local_ix_start, int slab_nx, int
       }
     }
   } // End looping through k box
-}
+}

src/core/stellar_feedback.c

@@ -183,7 +183,8 @@ double get_total_SN_energy(void)
 double get_SN_energy_PopIII(
   int i_burst,
   int snapshot,
-  int SN_type) // SN_type = 0 -> CC, 1 -> PISN (Pop. III have higher masses so we need to account also for PISN!)
+  int SN_type) // SN_type = 0 -> CC, 1 -> PISN 
+               // Pop. III have higher masses so we need to account also for PISN!
 {
   double NumberPISN = run_globals.NumberPISN;
   double NumberSNII = run_globals.NumberSNII;
@@ -194,7 +195,7 @@ double get_SN_energy_PopIII(
     double CC_Fraction = CCSN_PopIII_Fraction(i_burst, snapshot, 0);
     return Enova * CC_Fraction * NumberSNII * 1e10 /
            run_globals.params
-             .Hubble_h; // result in erg * (1e10 Msol / h). You will need to multiply this per mass in internal units.
+             .Hubble_h; // result in erg * (1e10 Msol / h)
   }
   // PISN (feedback here is contemporaneous).
   else if (SN_type == 1) {
@@ -210,4 +211,4 @@ double get_SN_energy_PopIII(
     return 0;
   }
 }
-#endif
+#endif

src/physics/evolve.c

@@ -177,7 +177,4 @@ void passively_evolve_ghost(galaxy_t* gal, int snapshot)
 
   if (!Flag_IRA)
     delayed_supernova_feedback(gal, snapshot);
-
-  // if (Flag_Metals == true) // You are updating this function to test why probability is decreasing in some cells
-  // calc_metal_bubble(gal, snapshot);
 }

src/physics/supernova_feedback.c

@@ -511,8 +511,6 @@ void calc_metal_bubble(galaxy_t* gal, int snapshot) // result in internal units
       gal->RmetalBubble = gal->PrefactorBubble * pow(IGM_density, -0.2) *
                           pow((gal->TimeBubble - run_globals.LTTime[snapshot] * time_unit), 0.4);
   }
-  if (gal->RmetalBubble > 10.) // Add this to check that there are no huge and unphysical bubbles.
-    mlog("StrangeBubble = %f, Prefactor = %f", MLOG_MESG, gal->RmetalBubble, gal->PrefactorBubble);
 
   if (gal->Type == 0) {
     for (int i_burst = 0; i_burst < n_bursts; i_burst++) {