diff --git a/src/core/BrightnessTemperature.c b/src/core/BrightnessTemperature.c
index d122aa3c..960cb309 100644
--- a/src/core/BrightnessTemperature.c
+++ b/src/core/BrightnessTemperature.c
@@ -18,7 +18,7 @@
  * Written by Bradley Greig.
  */
 
-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) // Added the computation of 21cm Tb without Pop. III (useful to better highlight the contribution of Pop III stars)
 {
 
   int ii, jj, kk, i_real, i_padded, iii;
diff --git a/src/core/PopIII.c b/src/core/PopIII.c
index 5edf768c..079bf692 100644
--- a/src/core/PopIII.c
+++ b/src/core/PopIII.c
@@ -202,7 +202,7 @@ double getIMF_massweighted(double StarMass)
   return StarMass * getIMF(StarMass);
 }
 
-double get_StellarAge(double StarMass) //Star Mass in log10(Msol) to get tstar in log10(yr). Use this so you can do a linear interpolation! 
+double get_StellarAge(double StarMass) //Star Mass in log10(Msol) to get tstar in log10(yr).  
 {
   int PopIIIAgePrescription = run_globals.params.physics.PopIIIAgePrescription; 
   if ((PopIIIAgePrescription == 1) && (run_globals.params.physics.MmaxIMF > 500.0)) {
@@ -432,11 +432,11 @@ double CCSN_PopIII_Number(int i_burst, int curr_snap, int flagMW)
     m_max = MmaxSnII;
     }
   
-  if (m_min > MmaxSnII) //There are no CCSN in this snapshot!
+  if (m_min > MmaxSnII) // There are no CCSN in this snapshot!
     return 0.0;
   
-  else if (m_max < MminSnII) //There are no CCSN in this snapshot!
-    return 0.0; //Maybe put -1 or a break because this condition should stop your while loop!
+  else if (m_max < MminSnII) // There are no CCSN in this snapshot!
+    return 0.0; // Maybe put -1 or a break because this condition should stop your while loop!
     
   else {
     if (m_max > MmaxSnII) // Firstly, you are only interested in stars in the CCSN mass range
@@ -506,7 +506,6 @@ double CCSN_PopIII_Yield(int i_burst, int curr_snap, int yield_type) //0 = Tot,
   double time_unit = run_globals.units.UnitTime_in_Megayears / run_globals.params.Hubble_h * 1e6; //You need result in yrs
   double m_max;
   
-  //double TotalMassSN = run_globals.MassPISN + run_globals.MassSNII;
   double TotalMassSN = run_globals.MassSNII;
   
   double DeltaTime = (LTTime[curr_snap - i_burst] - LTTime[curr_snap]) * time_unit;
@@ -552,15 +551,17 @@ double PISN_PopIII_Yield(int yield_type) //Yield_type = 0 -> Recycling, 1 -> Met
   //Remember that PISN feedback is always contemporaneous and they leave no remnants!
   double MassPISN = run_globals.MassPISN;
   double MassSNII = run_globals.MassSNII;
-  //double PISN_MassFraction = MassPISN / (MassSNII + MassPISN);
-  double PISN_MassFraction = MassPISN / MassPISN; // LINE ABOVE IS PROBABLY A MISTAKE!
+  double PISN_MassFraction = MassPISN / MassPISN; 
   if (yield_type == 0) 
     return PISN_MassFraction; 
   if (yield_type == 1)
     return PISN_MassFraction / 2.0;      
 }
 
-double NLBias(double Dist_Radius, double Halo_Mass, double redshift) //Fitting function written to match the results of Iliev+03 and Dijkstra+08. Input in internal units.
+double NLBias(double Dist_Radius, double Halo_Mass, double redshift) 
+// Fitting function written to match the results of Iliev+03 and Dijkstra+08. 
+// Input in internal units. We don't use this in the small box where we have a grid resolution of about 200ckpc
+// We want to use it for smaller grid resolution (about 1Mpc)
 {
   double Alpha_ind = run_globals.params.physics.AlphaCluster;
   double Beta_ind = run_globals.params.physics.BetaCluster;
diff --git a/src/core/galaxies.c b/src/core/galaxies.c
index e4f0af18..f85d1fa1 100644
--- a/src/core/galaxies.c
+++ b/src/core/galaxies.c
@@ -78,7 +78,7 @@ galaxy_t* new_galaxy(int snapshot, unsigned long halo_ID)
   gal->Metal_Probability = 0.0;
   gal->Metals_IGM = 0.0;
   gal->Gas_IGM = 0.0;
-  gal->Metallicity_IGM = -50.0; // Better to initialize to a very negative value.
+  gal->Metallicity_IGM = -50.0; 
   gal->RmetalBubble = 0.0;
   gal->PrefactorBubble = 0.0;
   gal->TimeBubble = 0.0;
@@ -87,7 +87,7 @@ galaxy_t* new_galaxy(int snapshot, unsigned long halo_ID)
   gal->Flag_ExtMetEnr = 0;
   gal->GalMetal_Probability = gsl_rng_uniform(run_globals.random_generator);
   
-  if (run_globals.params.Flag_IncludeMetalEvo == false) //If you don't have the external metal enrichment you have to initialize this variable (you were doing that in evolve.c)
+  if (run_globals.params.Flag_IncludeMetalEvo == false) // If you don't have the external metal enrichment all galaxies will start as pristine (Pop.III forming)
     gal-> Galaxy_Population == 3;
 #else //If you are not computing PopIII , all the galaxies are PopII. Again you need to initialize the variable otherwise star_formation.c will fail!
   gal->Galaxy_Population == 2;
diff --git a/src/core/magnitudes.c b/src/core/magnitudes.c
index a3a9ec5f..d67e43a9 100644
--- a/src/core/magnitudes.c
+++ b/src/core/magnitudes.c
@@ -457,7 +457,7 @@ void init_magnitudes(void)
 #if USE_MINI_HALOS
     int IMF_Type = run_globals.params.physics.PopIII_IMF;
     if (IMF_Type == 1)
-      strcat(fnameIII, "/Sal500_001.hdf5"); // Expand here once you have more IMFs
+      strcat(fnameIII, "/Sal500_001.hdf5"); 
     else if (IMF_Type == 2)
       strcat(fnameIII, "/Sal500_050.hdf5");
     else if (IMF_Type == 3)
diff --git a/src/core/metal_evo.c b/src/core/metal_evo.c
index 713af215..986af911 100644
--- a/src/core/metal_evo.c
+++ b/src/core/metal_evo.c
@@ -131,135 +131,6 @@ void construct_metal_grids(int snapshot, int local_ngals)
             
               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)
                 buffer_metals[ind] += (4.0 / 3.0 * M_PI * pow((gal->RmetalBubble) * (1 + redshift), 3.0)); // cMpc/h (same units of cell volume)
-                  
-                  //if (gal->RmetalBubble > 0.62 * (box_size / MetalGridDim)) {
-                /*if ((gal->RmetalBubble * (1 + redshift)) > 0.62 * (box_size / MetalGridDim)) { 
-                  double Excess_volume = (4.0 / 3.0 * M_PI * pow((gal->RmetalBubble) * (1 + redshift), 3.0)) - pow((box_size / MetalGridDim), 3.0);
-   
-                  // ixplus and minus needed for your new modification
-          
-                  int ixplus = (int)(pos_to_ngp(gal->Pos[0] + (float)(box_size / MetalGridDim), box_size, MetalGridDim) - slab_ix_start_metals[i_r]);
-                  int ixminus;
-                  
-                  if ((gal->Pos[0] - (float)(box_size / MetalGridDim)) < 0.0)
-                    ixminus = (int)(pos_to_ngp((box_size + gal->Pos[0] - (float)(box_size / MetalGridDim)), box_size, MetalGridDim) - slab_ix_start_metals[i_r]);
-                  else
-                    ixminus = (int)(pos_to_ngp(gal->Pos[0] - (float)(box_size / MetalGridDim), box_size, MetalGridDim) - slab_ix_start_metals[i_r]);
-                  
-                  if (ixplus >= slab_nix_metals[i_r])
-                    ixplus = 0;
-                   
-                  if (ixminus >= slab_nix_metals[i_r])
-                    ixminus = 0;
-                    
-                  if (ixplus < 0)
-                    ixplus = slab_nix_metals[i_r] - 1;
-                   
-                  if (ixminus < 0)
-                    ixminus = slab_nix_metals[i_r] - 1;
-          
-                  assert((ixplus < slab_nix_metals[i_r]) && (ixplus >= 0));
-                  assert((ixminus < slab_nix_metals[i_r]) && (ixminus >= 0));
-                  
-                  int iyplus = pos_to_ngp(gal->Pos[1] + (float)(box_size / MetalGridDim), box_size, MetalGridDim);
-                  int iyminus;
-                  
-                  if ((gal->Pos[1] - (float)(box_size / MetalGridDim)) < 0.0)
-                    iyminus = (int)(pos_to_ngp((box_size + gal->Pos[1] - (float)(box_size / MetalGridDim)), box_size, MetalGridDim));
-                  else
-                    iyminus = (int)(pos_to_ngp(gal->Pos[1] - (float)(box_size / MetalGridDim), box_size, MetalGridDim));
-                    
-                  int izplus = pos_to_ngp(gal->Pos[2] + (float)(box_size / MetalGridDim), box_size, MetalGridDim);
-                  int izminus;
-                  
-                  if ((gal->Pos[2] - (float)(box_size / MetalGridDim)) < 0.0)
-                    izminus = (int)(pos_to_ngp((box_size + gal->Pos[2] - (float)(box_size / MetalGridDim)), box_size, MetalGridDim));
-                  else
-                    izminus = (int)(pos_to_ngp(gal->Pos[2] - (float)(box_size / MetalGridDim), box_size, MetalGridDim));
-                    
-                  assert((iyplus < MetalGridDim) && (iyplus >= 0));
-                  assert((iyminus < MetalGridDim) && (iyminus >= 0));
-                  
-                  assert((izplus < MetalGridDim) && (izplus >= 0));
-                  assert((izminus < MetalGridDim) && (izminus >= 0));
-                  
-                  int indxplus = grid_index(ixplus, iy, iz, MetalGridDim, INDEX_REAL); 
-                  int indxminus = grid_index(ixminus, iy, iz, MetalGridDim, INDEX_REAL); 
-                  int indxypp = grid_index(ixplus, iyplus, iz, MetalGridDim, INDEX_REAL); 
-                  int indxypm = grid_index(ixplus, iyminus, iz, MetalGridDim, INDEX_REAL); 
-                  int indxymp = grid_index(ixminus, iyplus, iz, MetalGridDim, INDEX_REAL); 
-                  int indxymm = grid_index(ixminus, iyminus, iz, MetalGridDim, INDEX_REAL); 
-                  int indxzpp = grid_index(ixplus, iy, izplus, MetalGridDim, INDEX_REAL); 
-                  int indxzpm = grid_index(ixplus, iy, izminus, MetalGridDim, INDEX_REAL);
-                  int indxzmp = grid_index(ixminus, iy, izplus, MetalGridDim, INDEX_REAL); 
-                  int indxzmm = grid_index(ixminus, iy, izminus, MetalGridDim, INDEX_REAL);
-                   
-                  //if (indplus >= slab_nix_metals[i_r] * MetalGridDim * MetalGridDim)
-                  //  indplus = 0;
-                  
-                  //if (indminus < 0)
-                  //  indminus = slab_nix_metals[i_r] * MetalGridDim * MetalGridDim - 1;
-                  
-                  int indyplus = grid_index(ix, iyplus, iz, MetalGridDim, INDEX_REAL); 
-                  int indyminus = grid_index(ix, iyminus, iz, MetalGridDim, INDEX_REAL);
-                  int indzplus = grid_index(ix, iy, izplus, MetalGridDim, INDEX_REAL); 
-                  int indzminus = grid_index(ix, iy, izminus, MetalGridDim, INDEX_REAL); 
-                  int indyzpp = grid_index(ix, iyplus, izplus, MetalGridDim, INDEX_REAL); 
-                  int indyzpm = grid_index(ix, iyplus, izminus, MetalGridDim, INDEX_REAL);
-                  int indyzmp = grid_index(ix, iyminus, izplus, MetalGridDim, INDEX_REAL);
-                  int indyzmm = grid_index(ix, iyminus, izminus, MetalGridDim, INDEX_REAL); 
-                  
-                  assert((indxplus >= 0) && (indxplus < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxminus >= 0) && (indxminus < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxypp >= 0) && (indxypp < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxypm >= 0) && (indxypm < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxymp >= 0) && (indxymp < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxymm >= 0) && (indxymm< slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxzpm >= 0) && (indxzpm < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxzmp >= 0) && (indxzmp < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indxzmm >= 0) && (indxzmm < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  
-                  assert((indyplus >= 0) && (indyplus < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indyminus >= 0) && (indyminus < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indzplus >= 0) && (indzplus < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indzminus >= 0) && (indzminus < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indyzpp >= 0) && (indyzpp < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indyzpm >= 0) && (indyzpm < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indyzmp >= 0) && (indyzmp < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  assert((indyzmm >= 0) && (indyzmm < slab_nix_metals[i_r] * MetalGridDim * MetalGridDim));
-                  
-                  //Adiacent cells in the same axis 
-                  buffer_metals[indzplus] += 0.072 * Excess_volume;
-                  buffer_metals[indzminus] += 0.072 * Excess_volume;
-                  buffer_metals[indyplus] += 0.072 * Excess_volume;
-                  buffer_metals[indyminus] += 0.072 * Excess_volume;
-                  //buffer_metals[ind + (MetalGridDim * MetalGridDim)] += 0.072 * Excess_volume;
-                  //buffer_metals[ind - (MetalGridDim * MetalGridDim)] += 0.072 * Excess_volume;
-                  buffer_metals[indxplus] += 0.072 * Excess_volume; //SLAB NOT JUST ix (still to check!)
-                  buffer_metals[indxminus] += 0.072 * Excess_volume;
-                  
-                  //Adiacent cells obliquos
-                  buffer_metals[indyzpp] += 0.0473 * Excess_volume;
-                  buffer_metals[indyzpm] += 0.0473 * Excess_volume;
-                  buffer_metals[indyzmm] += 0.0473 * Excess_volume;
-                  buffer_metals[indyzmp] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind + (MetalGridDim * MetalGridDim) + 1] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind - (MetalGridDim * MetalGridDim) + 1] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind + (MetalGridDim * MetalGridDim) - 1] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind - (MetalGridDim * MetalGridDim) - 1] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind + (MetalGridDim * MetalGridDim) + MetalGridDim] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind - (MetalGridDim * MetalGridDim) + MetalGridDim] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind + (MetalGridDim * MetalGridDim) - MetalGridDim] += 0.0473 * Excess_volume;
-                  //buffer_metals[ind - (MetalGridDim * MetalGridDim) - MetalGridDim] += 0.0473 * Excess_volume;
-                  buffer_metals[indxzpp] += 0.0473 * Excess_volume;
-                  buffer_metals[indxzmp] += 0.0473 * Excess_volume;
-                  buffer_metals[indxzpm] += 0.0473 * Excess_volume;
-                  buffer_metals[indxzmm] += 0.0473 * Excess_volume;
-                  buffer_metals[indxypp] += 0.0473 * Excess_volume;
-                  buffer_metals[indxymp] += 0.0473 * Excess_volume;
-                  buffer_metals[indxypm] += 0.0473 * Excess_volume;
-                  buffer_metals[indxymm] += 0.0473 * Excess_volume;
-                }*/
               }
 
               break;
@@ -288,13 +159,12 @@ void construct_metal_grids(int snapshot, int local_ngals)
             case prop_mass_ej_metals:
             
               if (gal->RmetalBubble >= 3 * gal->Rvir) // Add this condition to be consistent with above
-                buffer_metals[ind] += gal->MetalsEjectedGas; //Internal units (same of gas_cell)
+                buffer_metals[ind] += gal->MetalsEjectedGas; // Internal units (same of gas_cell)
                 
               break;
               
             case prop_mass_ej_gas:
             
-              //buffer_metals[ind] += (gal->EjectedGas - gal->HotGas - gal->ColdGas);
               buffer_metals[ind] -= (gal->HotGas + gal->ColdGas);
               if (gal->RmetalBubble >= 3 * gal->Rvir)
                 buffer_metals[ind] += gal->EjectedGas; // Add this condition to be consistent with above
@@ -544,10 +414,6 @@ void malloc_metal_grids()
 
   metal_grids_t* grids = &(run_globals.metal_grids);
 
-  // run_globals.NStoreSnapshots is set in `initialize_halo_storage` Leave it commented because you might use that in the future when looking for scaling relations
-  //run_globals.SnapshotDeltax = (float**)calloc((size_t)run_globals.NStoreSnapshots, sizeof(float*)); //?
-  //run_globals.SnapshotVel = (float**)calloc((size_t)run_globals.NStoreSnapshots, sizeof(float*));  //?
-
   grids->galaxy_to_slab_map_metals = NULL;
 
   grids->mass_metals = NULL;
diff --git a/src/core/read_grids.c b/src/core/read_grids.c
index b558dfff..dfd6dfbc 100644
--- a/src/core/read_grids.c
+++ b/src/core/read_grids.c
@@ -54,7 +54,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)
+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)
 {
   mlog("Smoothing the overdensity of the reionization grid into the metal grid...", MLOG_MESG);
   reion_grids_t* reiogrids = &(run_globals.reion_grids);
@@ -108,8 +108,10 @@ void smooth_Densitygrid_real(int snapshot) //Need this to put the overdensity in
   for (int ii = 0; ii < slab_n_real_metals; ii++) { 
     metalgrids->mass_IGM[ii] /= (resample_factorReal * resample_factorReal * resample_factorReal); // Average
     metalgrids->mass_IGM[ii] += 1.0; // overdensity is 1 + deltax
+    
     //To compute the mass we need to multiply by total baryon content 
     metalgrids->mass_IGM[ii] *= run_globals.rhocrit[snapshot] * run_globals.params.OmegaM * run_globals.params.BaryonFrac * pow((pixel_length_metals / (1.0 + run_globals.ZZ[snapshot])), 3.0);
+    
     //Add the contribution coming from galaxies (+ Ejected - Hot - Cold)
     metalgrids->mass_IGM[ii] += metalgrids->mass_gas[ii];
     }
diff --git a/src/core/save.c b/src/core/save.c
index be59c55b..808b7840 100644
--- a/src/core/save.c
+++ b/src/core/save.c
@@ -745,7 +745,7 @@ void create_master_file()
   }
   
 #if USE_MINI_HALOS
-  if (run_globals.params.Flag_IncludeMetalEvo) { // Here there might be some issue with little h!
+  if (run_globals.params.Flag_IncludeMetalEvo) { 
     const char* group_name = { "Units/MetalGrids" };
     group_id = H5Gcreate(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     H5LTset_attribute_string(file_id, group_name, "mass_metals", "solMass");
@@ -753,7 +753,7 @@ void create_master_file()
     H5LTset_attribute_string(file_id, group_name, "Zigm_box", "Zsol");
     H5LTset_attribute_string(file_id, group_name, "Probability_metals", "none");
     H5LTset_attribute_string(file_id, group_name, "N_bubbles", "none");
-    H5LTset_attribute_string(file_id, group_name, "R_ave", "cMpc"); // I believe are cMpc / h! But there might be something wrong! You might want to add Conversion grids!!
+    H5LTset_attribute_string(file_id, group_name, "R_ave", "cMpc"); 
     H5LTset_attribute_string(file_id, group_name, "R_max", "cMpc");
     H5LTset_attribute_string(file_id, group_name, "mass_IGM", "solMass");
     
@@ -928,20 +928,6 @@ void create_master_file()
           H5Lcreate_external(relative_source_file, "/", snap_group_id, "Grids", H5P_DEFAULT, H5P_DEFAULT);
           H5Fclose(source_file_id);
         }
-
-        // sprintf(source_ds, "Snap%03d/PowerSpectrum", run_globals.ListOutputSnaps[i_out]);
-        // if ((H5LTpath_valid(source_file_id, source_ds, FALSE)))
-        // {
-        //   sprintf(target_ds, "PowerSpectrum");
-        //   H5Lcreate_external(relative_source_file, source_ds, snap_group_id, target_ds, H5P_DEFAULT, H5P_DEFAULT);
-        // }
-
-        // sprintf(source_ds, "Snap%03d/RegionSizeDist", run_globals.ListOutputSnaps[i_out]);
-        // if ((H5LTpath_valid(source_file_id, source_ds, FALSE)))
-        // {
-        //   sprintf(target_ds, "RegionSizeDist");
-        //   H5Lcreate_external(relative_source_file, source_ds, snap_group_id, target_ds, H5P_DEFAULT, H5P_DEFAULT);
-        // }
       }
       
 #if USE_MINI_HALOS
@@ -954,20 +940,6 @@ void create_master_file()
           H5Lcreate_external(relative_source_file, "/", snap_group_id, "MetalGrids", H5P_DEFAULT, H5P_DEFAULT);
           H5Fclose(source_file_id);
         }
-
-        // sprintf(source_ds, "Snap%03d/PowerSpectrum", run_globals.ListOutputSnaps[i_out]);
-        // if ((H5LTpath_valid(source_file_id, source_ds, FALSE)))
-        // {
-        //   sprintf(target_ds, "PowerSpectrum");
-        //   H5Lcreate_external(relative_source_file, source_ds, snap_group_id, target_ds, H5P_DEFAULT, H5P_DEFAULT);
-        // }
-
-        // sprintf(source_ds, "Snap%03d/RegionSizeDist", run_globals.ListOutputSnaps[i_out]);
-        // if ((H5LTpath_valid(source_file_id, source_ds, FALSE)))
-        // {
-        //   sprintf(target_ds, "RegionSizeDist");
-        //   H5Lcreate_external(relative_source_file, source_ds, snap_group_id, target_ds, H5P_DEFAULT, H5P_DEFAULT);
-        // }
       }
 #endif
     }
diff --git a/src/core/stellar_feedback.c b/src/core/stellar_feedback.c
index ec621eaa..b90babf2 100644
--- a/src/core/stellar_feedback.c
+++ b/src/core/stellar_feedback.c
@@ -178,9 +178,9 @@ double get_total_SN_energy(void)
 }
 
 #if USE_MINI_HALOS
-// Adding stuff for Pop III feedback 
+// Stuff for Pop. III feedback 
 
-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!)
+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!)
 {
   double NumberPISN = run_globals.NumberPISN;
   double NumberSNII = run_globals.NumberSNII;
@@ -189,9 +189,9 @@ double get_SN_energy_PopIII(int i_burst, int snapshot, int SN_type) //SN_type =
   if (SN_type == 0) {
     Enova = ENOVA_CC; 
     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
+    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.
   }
-  //PISN (feedback here is contemporaneous), this part is probably useless
+  //PISN (feedback here is contemporaneous).
   if (SN_type == 1) {
     if (i_burst != 0) {
       mlog_error("PISN feedback is instantaneous!");
diff --git a/src/physics/cooling.c b/src/physics/cooling.c
index d8e20b7c..45dbd3c0 100644
--- a/src/physics/cooling.c
+++ b/src/physics/cooling.c
@@ -75,15 +75,15 @@ double gas_cooling(galaxy_t* gal)
 
     if (halo_type != 0) {
 
-      x = PROTONMASS * BOLTZMANN * Tvir / lambda;              // now this has units sec g/cm^3
+      x = PROTONMASS * BOLTZMANN * Tvir / lambda; // now this has units sec g/cm^3
       x /= (units->UnitDensity_in_cgs * units->UnitTime_in_s); // now in internal units
       
       if (halo_type == 1)
-        rho_r_cool = x / t_cool * 0.885;                         // 0.885 = 3/2 * mu, mu=0.59 for a fully ionized gas
+        rho_r_cool = x / t_cool * 0.885; // 0.885 = 3/2 * mu, mu=0.59 for a fully ionized gas
 
 #if USE_MINI_HALOS      
       else if (halo_type == 2)
-        rho_r_cool = x / t_cool * 1.83;                        // 1.83 = 3/2 * mu, mu = 1.22 for a fully neutral gas  
+        rho_r_cool = x / t_cool * 1.83;  // 1.83 = 3/2 * mu, mu = 1.22 for a fully neutral gas  
 #endif
 
       assert(rho_r_cool > 0);
diff --git a/src/physics/evolve.c b/src/physics/evolve.c
index 33e6938f..f55231f3 100644
--- a/src/physics/evolve.c
+++ b/src/physics/evolve.c
@@ -6,7 +6,6 @@
 #include "core/PopIII.h"
 #include "core/virial_properties.h"
 #include "core/misc_tools.h"
-//#include <gsl/gsl_rng.h>
 #endif
 #include "infall.h"
 #include "meraxes.h"
diff --git a/src/physics/infall.c b/src/physics/infall.c
index 91cb2059..3b2b5b76 100644
--- a/src/physics/infall.c
+++ b/src/physics/infall.c
@@ -96,7 +96,7 @@ void add_infall_to_hot(galaxy_t* central, double infall_mass)
     central->HotGas += infall_mass;
 #if USE_MINI_HALOS
     if (Flag_Metals == true) {
-      if (central->Flag_ExtMetEnr == 1)
+      if (central->Flag_ExtMetEnr == 1) // If the halo is externally enriched, it will accrete polluted gas (metals).
         central->MetalsHotGas += infall_mass * central->Metallicity_IGM;
     }
 #endif
diff --git a/src/physics/mergers.c b/src/physics/mergers.c
index ff7ff8f8..46447324 100644
--- a/src/physics/mergers.c
+++ b/src/physics/mergers.c
@@ -208,9 +208,6 @@ void merge_with_target(galaxy_t* gal, int* dead_gals, int snapshot)
   // If I have a Merger between Pop III and Pop II the result is a Pop. II. Actually I should compute metallicity
   // TODO: this could be improved in the future!
   
-  //if (gal->Galaxy_Population == 2) 
-  //  parent->Galaxy_Population = 2; // Since now you updated the internal enrichment probably this condition is useless
-  
   if (gal->RmetalBubble > parent->RmetalBubble) {
     
     parent->RmetalBubble = gal->RmetalBubble; // This is to account the evolution of metal bubbles after a merger event
diff --git a/src/physics/supernova_feedback.c b/src/physics/supernova_feedback.c
index 9c46bb74..1a042cc6 100644
--- a/src/physics/supernova_feedback.c
+++ b/src/physics/supernova_feedback.c
@@ -296,13 +296,13 @@ void delayed_supernova_feedback(galaxy_t* gal, int snapshot)
       new_metals += m_stars_II * get_metal_yield(i_burst, metallicity);
 #if USE_MINI_HALOS
       new_metals += m_stars_III * CCSN_PopIII_Yield(i_burst, snapshot, 1) * MassSNII;
-      m_remnant += m_stars_III * CCSN_PopIII_Yield(i_burst, snapshot, 2) * MassSNII; // Pop. III will have remnants
+      m_remnant += m_stars_III * CCSN_PopIII_Yield(i_burst, snapshot, 2) * MassSNII; // Remnants from Pop. III 
 #endif
       // Calculate SNII energy
       
       sn_energy_II += get_SN_energy(i_burst, metallicity) * m_stars_II; 
 #if USE_MINI_HALOS
-      sn_energy_III += get_SN_energy_PopIII(i_burst, snapshot, 0) * m_stars_III; // Tis is DeltaM reheat (eq.16 Mutch+16) * ENOVA
+      sn_energy_III += get_SN_energy_PopIII(i_burst, snapshot, 0) * m_stars_III; // This is DeltaM reheat (eq.16 Mutch+16) * ENOVA
 #endif
     }
   }
@@ -311,7 +311,9 @@ void delayed_supernova_feedback(galaxy_t* gal, int snapshot)
   sn_energy_II *= calc_sn_ejection_eff(gal, snapshot, 2);
 #if USE_MINI_HALOS
   m_reheat_III = calc_sn_reheat_eff(gal, snapshot, 3) * sn_energy_III / ENOVA_CC; 
-  sn_energy_III *= (calc_sn_ejection_eff(gal, snapshot, 3) * NumberSNII * 1e10 / run_globals.params.Hubble_h / energy_unit); //Maybe for the SN ejection efficiency is more important to distinguish between PISN/CC rather than Pop.III/II 
+  sn_energy_III *= (calc_sn_ejection_eff(gal, snapshot, 3) * NumberSNII * 1e10 / run_globals.params.Hubble_h / energy_unit); 
+  
+  // Maybe for the SN ejection efficiency is more important to distinguish between PISN/CC rather than Pop.III/II 
 #endif
   m_recycled = m_recycled_II + m_recycled_III;
   m_reheat = m_reheat_II + m_reheat_III;
@@ -406,7 +408,6 @@ void contemporaneous_supernova_feedback(galaxy_t* gal,
       }
     else if (gal->Galaxy_Population == 3){
       *m_recycled = *m_stars * (CCSN_PopIII_Yield(0, snapshot, 0)) * MassSNII;
-      //*new_metals = *m_stars * CCSN_PopIII_Yield(0, snapshot, 1); // Shouldn't have also MassSNII here?
       *new_metals = *m_stars * CCSN_PopIII_Yield(0, snapshot, 1) * MassSNII;
       *m_remnant = *m_stars * (MassBHs + CCSN_PopIII_Yield(0, snapshot, 2) * MassSNII);
       if (MassPISN > 0) { // Account for PISN
@@ -442,8 +443,8 @@ void contemporaneous_supernova_feedback(galaxy_t* gal,
 #if USE_MINI_HALOS
   }
   else if (gal->Galaxy_Population == 3){ 
-    sn_energy = (*m_stars) * (get_SN_energy_PopIII(0, snapshot, 0) + get_SN_energy_PopIII(0, snapshot, 1)); //erg
-    sn_energy /= energy_unit; //Convert this because you need in internal units it for m_ejected
+    sn_energy = (*m_stars) * (get_SN_energy_PopIII(0, snapshot, 0) + get_SN_energy_PopIII(0, snapshot, 1)); // erg
+    sn_energy /= energy_unit; // Convert this because you need in internal units it for m_ejected
     *m_reheat = calc_sn_reheat_eff(gal, snapshot, 3) * (*m_stars) * (get_SN_energy_PopIII(0, snapshot, 1) / ENOVA_PISN + get_SN_energy_PopIII(0, snapshot, 0) / ENOVA_CC );
     sn_energy *= calc_sn_ejection_eff(gal, snapshot, 3);  
     }
@@ -455,8 +456,6 @@ void contemporaneous_supernova_feedback(galaxy_t* gal,
   if (*m_reheat > gal->ColdGas)
     *m_reheat = gal->ColdGas;
     
-  // You might add here a condition for m_remnant!
-    
   // attenuate the star formation if necessary, so that we are being consistent
   // if (*m_reheat + *m_stars - *m_recycled > gal->ColdGas)
   if (*m_reheat + *m_stars > gal->ColdGas) {
@@ -510,7 +509,7 @@ void calc_metal_bubble(galaxy_t* gal, int snapshot) // result in internal units
     else
       gal->RmetalBubble = gal->PrefactorBubble * pow(IGM_density, -0.2) * pow((gal->TimeBubble - run_globals.LTTime[snapshot] * time_unit), 0.4);
   }  
-  if (gal->RmetalBubble > 10.)
+  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) {
@@ -519,10 +518,10 @@ void calc_metal_bubble(galaxy_t* gal, int snapshot) // result in internal units
       double m_stars_III = gal->NewStars_III[i_burst];
       double m_stars = m_stars_II + m_stars_III;
     
-      // Compute the SN energy that drives the metal bubble. Should you include the ejection efficiency?  
+      // Compute the SN energy that drives the metal bubble.
       if (m_stars_II > 1e-10) {
         double metallicity = calc_metallicity(m_stars, gal->NewMetals[i_burst]);
-        sn_energy += m_stars_II * get_SN_energy(0, metallicity) * energy_unit * calc_sn_ejection_eff(gal, snapshot, 2); // Are you sure of the unit? (You want this in erg)
+        sn_energy += m_stars_II * get_SN_energy(0, metallicity) * energy_unit * calc_sn_ejection_eff(gal, snapshot, 2); 
         }
       else if (m_stars_III > 1e-10) {
         if (i_burst == 0) // You have both CC and PISN
@@ -541,7 +540,9 @@ void calc_metal_bubble(galaxy_t* gal, int snapshot) // result in internal units
         }
         else
           gal->Radii[n_bursts - i_burst] = 0.0;
-        if (gal->Radii[n_bursts - i_burst] > gal->RmetalBubble) { //Look if one of the new bubbles is bigger than RmetalBubble
+        if (gal->Radii[n_bursts - i_burst] > gal->RmetalBubble) { 
+          //Look if one of the new bubbles is bigger than RmetalBubble 
+          //and in this case this will be the new metal bubble associated to the galaxy.
           gal->RmetalBubble = gal->Radii[n_bursts - i_burst];
           gal->PrefactorBubble = gal->Prefactor[n_bursts - i_burst];
           gal->TimeBubble = gal->Times[n_bursts - i_burst];