Skip to content

Commit

Permalink
Remove or added comments
Browse files Browse the repository at this point in the history 
Make this for publishing the code
  • Loading branch information
@EMventura @qyx268
EMventura authored and qyx268 committed Jan 20, 2024
1 parent 8752201 commit f25dc03
Show file tree
Hide file tree
Showing 13 changed files with 40 additions and 202 deletions.
2 changes: 1 addition & 1 deletion src/core/BrightnessTemperature.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
17 changes: 9 additions & 8 deletions src/core/PopIII.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)) {
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 2 additions & 2 deletions src/core/galaxies.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
2 changes: 1 addition & 1 deletion src/core/magnitudes.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
136 changes: 1 addition & 135 deletions src/core/metal_evo.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 3 additions & 1 deletion src/core/read_grids.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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];
}
Expand Down
Loading

0 comments on commit f25dc03

Please sign in to comment.