diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index ec1b01b45167bc7d6cdb777191ca04f5e501640b..b4af701343ae309bee2d9858f8620e1d4d79a44d 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -11,7 +11,7 @@ Cosmology:
Omega_lambda: 0.693
Omega_b: 0.0455
h: 0.6777
- a_begin: 0.0078125
+ a_begin: 0.25
a_end: 1.0
# Parameters governing the time integration
@@ -66,21 +66,34 @@ GrackleCooling:
ConvergenceLimit: 1e-2
EagleCooling:
- filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
+ filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
+ reionisation_redshift: 8.898
EAGLEChemistry:
- InitMetallicity: 0.
- InitAbundance_Hydrogen: 0.752
- InitAbundance_Helium: 0.248
- InitAbundance_Carbon: 0.000
- InitAbundance_Nitrogen: 0.000
- InitAbundance_Oxygen: 0.000
- InitAbundance_Neon: 0.000
- InitAbundance_Magnesium: 0.000
- InitAbundance_Silicon: 0.000
- InitAbundance_Iron: 0.000
+ InitMetallicity: 0.014
+ InitAbundance_Hydrogen: 0.706
+ InitAbundance_Helium: 0.280
+ InitAbundance_Carbon: 0.00207
+ InitAbundance_Nitrogen: 8.35e-4
+ InitAbundance_Oxygen: 0.00549
+ InitAbundance_Neon: 0.00141
+ InitAbundance_Magnesium: 5.91e-4
+ InitAbundance_Silicon: 6.83e-4
+ InitAbundance_Iron: 0.0011
CalciumOverSilicon: 0.0941736
SulphurOverSilicon: 0.6054160
+ #InitMetallicity: 0.
+ #InitAbundance_Hydrogen: 0.752
+ #InitAbundance_Helium: 0.248
+ #InitAbundance_Carbon: 0.000
+ #InitAbundance_Nitrogen: 0.000
+ #InitAbundance_Oxygen: 0.000
+ #InitAbundance_Neon: 0.000
+ #InitAbundance_Magnesium: 0.000
+ #InitAbundance_Silicon: 0.000
+ #InitAbundance_Iron: 0.000
+ #CalciumOverSilicon: 0.0941736
+ #SulphurOverSilicon: 0.6054160
GearChemistry:
InitialMetallicity: 0.01295
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index b6d4e565990e832c46f4902491ad0bc8c4b41699..f330d2660da70b288dcb5c83d4e15e022cc9828f 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -40,10 +40,12 @@
/**
* @brief Return a string containing the name of a given #chemistry_element.
*/
-__attribute__((always_inline)) INLINE static const char*
+//__attribute__((always_inline)) INLINE static const char*
+__attribute__((always_inline)) INLINE const char*
chemistry_get_element_name(enum chemistry_element elem) {
- static const char* chemistry_element_names[chemistry_element_count] = {
+ //static const char* chemistry_element_names[chemistry_element_count] = {
+ const char* chemistry_element_names[chemistry_element_count] = {
"Hydrogen", "Helium", "Carbon", "Nitrogen", "Oxygen",
"Neon", "Magnesium", "Silicon", "Iron"};
@@ -172,22 +174,22 @@ static INLINE void chemistry_print_backend(
* @param p particle struct
* @param elem enum value of element
*/
-__attribute__((always_inline)) INLINE static double chemistry_get_number_density(const struct part* restrict p, enum chemistry_element elem, const struct phys_const* restrict internal_const) {
+__attribute__((always_inline)) INLINE static double chemistry_get_number_density(const struct part* restrict p, const struct cosmology* restrict cosmo, enum chemistry_element elem, const struct phys_const* restrict internal_const) {
double number_density;
int atomic_number;
switch(elem){
- case chemistry_element_H : atomic_number = 1;
- case chemistry_element_He: atomic_number = 4;
- case chemistry_element_C : atomic_number = 12;
- case chemistry_element_N : atomic_number = 14;
- case chemistry_element_O : atomic_number = 16;
- case chemistry_element_Ne: atomic_number = 20;
- case chemistry_element_Mg: atomic_number = 24;
- case chemistry_element_Si: atomic_number = 28;
- case chemistry_element_Fe: atomic_number = 56;
+ case chemistry_element_H : atomic_number = 1; break;
+ case chemistry_element_He: atomic_number = 4; break;
+ case chemistry_element_C : atomic_number = 12; break;
+ case chemistry_element_N : atomic_number = 14; break;
+ case chemistry_element_O : atomic_number = 16; break;
+ case chemistry_element_Ne: atomic_number = 20; break;
+ case chemistry_element_Mg: atomic_number = 24; break;
+ case chemistry_element_Si: atomic_number = 28; break;
+ case chemistry_element_Fe: atomic_number = 56; break;
}
double element_mass = internal_const->const_proton_mass*atomic_number;
- number_density = p->chemistry_data.metal_mass_fraction[elem]*hydro_get_comoving_density(p)/element_mass;
+ number_density = p->chemistry_data.metal_mass_fraction[elem]*hydro_get_physical_density(p,cosmo)/element_mass;
return number_density;
}
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index c5e2c95479f0484c21f14e84f820215c8ac05612..ef3f104a6ca9aeaad92c4578afee6ac68f97d85e 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -33,7 +33,7 @@
#include <hdf5.h>
/* Local includes. */
-// #include "cooling_struct.h"
+#include "cooling_struct.h"
#include "error.h"
#include "hydro.h"
#include "chemistry.h"
@@ -129,7 +129,7 @@ __attribute__((always_inline)) INLINE void get_index_1d_therm(float *table, int
*dx = 1;
} else {
*i = (int)floor(((float)x - table[0]) * dxm1);
- //printf("i, x, dxm1: %d, %f, %f, %f\n", *i, (float) x, table[0], dxm1);
+ printf("i, x, dxm1: %d, %f, %f, %f\n", *i, (float) x, table[0], dxm1);
*dx = ((float)x - table[*i]) * dxm1;
}
}
@@ -315,6 +315,141 @@ __attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int
return result;
}
+
+//int static Silicon_SPH_Index = -1;
+//int static Calcium_SPH_Index = -1;
+//int static Sulphur_SPH_Index = -1;
+//
+//int static Silicon_CoolHeat_Index = -1;
+//int static Calcium_CoolHeat_Index = -1;
+//int static Sulphur_CoolHeat_Index = -1;
+
+static int first_call = 0;
+
+inline int set_cooling_SolarAbundances(const float *element_abundance,
+ double *cooling_element_abundance,
+ const struct cooling_function_data* restrict cooling) {
+ int i, index;
+ int Silicon_SPH_Index = -1;
+ int Calcium_SPH_Index = -1;
+ int Sulphur_SPH_Index = -1;
+
+ int Silicon_CoolHeat_Index = -1;
+ int Calcium_CoolHeat_Index = -1;
+ int Sulphur_CoolHeat_Index = -1;
+ for (int k = 0; k < cooling->N_Elements; k++){
+ printf("element, abundances %d, %.5e\n", k, element_abundance[k]);
+ }
+
+ float *cooling_ElementAbundance_SOLARM1 = malloc(cooling->N_SolarAbundances*sizeof(float));
+
+ //if (first_call == 0) {
+ /* determine (inverse of) solar abundance of these elements */
+ for (i = 0; i < cooling->N_Elements; i++) {
+ index =
+ get_element_index(cooling->SolarAbundanceNames,
+ cooling->N_SolarAbundances, cooling->ElementNames[i]);
+
+ if (index < 0) error("Eagle cooling.h index out of bounds");
+
+ //for (int j = 0; j < cooling->N_SolarAbundances; j++) {
+ // if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) == 0)
+ // index = j;
+ //}
+
+ index = cooling->SolarAbundanceNamePointers[i];
+
+ if(cooling->SolarAbundances[index] != 0) cooling_ElementAbundance_SOLARM1[i] = 1. / cooling->SolarAbundances[index];
+ else cooling_ElementAbundance_SOLARM1[i] = 0.0;
+
+ //index = ElementNamePointers[i];
+
+ //if (index < 0 && ThisTask == 0)
+ // printf("[bg_cooling] element not found %s\n", cooling_ElementNames[i]);
+ }
+
+ /* Sulphur tracks Silicon: may choose not to follow Sulphur as SPH element
+ */
+ /* Same is true for Calcium */
+ /* We will assume the code tracks Silicon, and may need to scale Calcium and
+ * Sulphur accordingly */
+
+ Silicon_SPH_Index = element_index("Silicon",cooling);
+ Calcium_SPH_Index = element_index("Calcium",cooling);
+ Sulphur_SPH_Index = element_index("Sulphur",cooling);
+
+ Silicon_CoolHeat_Index =
+ get_element_index(cooling->ElementNames, cooling->N_Elements, "Silicon");
+ Calcium_CoolHeat_Index =
+ get_element_index(cooling->ElementNames, cooling->N_Elements, "Calcium");
+ Sulphur_CoolHeat_Index =
+ get_element_index(cooling->ElementNames, cooling->N_Elements, "Sulphur");
+
+ if (Silicon_CoolHeat_Index == -1 || Calcium_CoolHeat_Index == -1 ||
+ Sulphur_CoolHeat_Index == -1) {
+ //if (ThisTask == 0)
+ // printf("[bg_cooling] error: did not find Si or Ca or S??\n");
+ //endrun(-1233);
+ error("Si, Ca, or S index out of bound\n");
+ }
+
+ first_call = 1;
+ //}
+
+ int sili_index;
+ for (i = 0; i < cooling->N_Elements; i++) {
+ if (strcmp(chemistry_get_element_name((enum chemistry_element) i), "Silicon") == 0) sili_index = i;
+ }
+
+ printf("Eagle cooling.h calcium cool index, sph index %d, %d\n", Calcium_CoolHeat_Index, Calcium_SPH_Index);
+ printf("Eagle cooling.h sulphur cool index, sph index %d, %d\n", Sulphur_CoolHeat_Index, Sulphur_SPH_Index);
+ for (i = 0; i < cooling->N_Elements; i++) {
+ if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1)
+ /* SPH does not track Calcium: use Si abundance */
+ if (Silicon_SPH_Index == -1)
+ cooling_element_abundance[i] = 0.0;
+ else{
+ cooling_element_abundance[i] =
+ element_abundance[Silicon_SPH_Index] *
+ cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
+ printf("Eagle cooling.h index %d calcium abundance %.5e\n", i, cooling_element_abundance[i]);
+ }
+ else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index == -1)
+ /* SPH does not track Sulphur: use Si abundance */
+ if (Silicon_SPH_Index == -1)
+ cooling_element_abundance[i] = 0.0;
+ else{
+ cooling_element_abundance[i] =
+ element_abundance[Silicon_SPH_Index] *
+ cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
+ printf("Eagle cooling.h index %d sulphur abundance %.5e\n", i, cooling_element_abundance[i]);
+ }
+ else{
+ /* introduce index to replace ElementNamePointers array, see function MakeNamePointers in eagle_init_cool.c */
+ //int element_index;
+ //if (strcmp(cooling->ElementNames[i], "Sulphur") == 0 ||
+ // strcmp(cooling->ElementNames[i], "Calcium") ==
+ // 0) /* These elements are tracked! */
+ // element_index = -1 * sili_index;
+ //else {
+ // for (int j = 0; j < cooling->N_Elements; j++) {
+ // if (strcmp(cooling->ElementNames[i], chemistry_get_element_name((enum chemistry_element) j)) == 0)
+ // element_index = j;
+ // //printf("Eagle cooling.h n_elements, element name, get element name, element index %d, %s, %s, %d\n", cooling->N_Elements, cooling->ElementNames[i], chemistry_get_element_name((enum chemistry_element) j), element_index);
+ // }
+ //}
+ //cooling_element_abundance[i] = element_abundance[element_index] *
+ // cooling_ElementAbundance_SOLARM1[i];
+ //printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",element_index,cooling->ElementNames[i], element_abundance[element_index],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
+ cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
+ cooling_ElementAbundance_SOLARM1[i];
+ printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
+ }
+ }
+
+ return 0;
+}
+
/*
* ----------------------------------------------------------------------
@@ -375,9 +510,8 @@ __attribute__((always_inline)) INLINE void get_redshift_index(float z, int *z_in
*
*/
__attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(const struct part* restrict p, const struct cooling_function_data* restrict cooling, const struct cosmology* restrict cosmo, const struct phys_const *internal_const) {
- float d_z;
- double inn_h = chemistry_get_number_density(p,chemistry_element_H,internal_const)*cooling->number_density_scale;
- double inHe = chemistry_get_number_density(p,chemistry_element_He,internal_const)*cooling->number_density_scale; // chemistry data
+ double inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const)*cooling->number_density_scale;
+ double inHe = chemistry_get_number_density(p,cosmo,chemistry_element_He,internal_const)*cooling->number_density_scale; // chemistry data
double u = hydro_get_physical_internal_energy(p,cosmo)*cooling->internal_energy_scale;
int u_i, He_i, n_h_i;
@@ -386,30 +520,34 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(cons
float z = cosmo->z,dz;
int z_index;
- printf("Eagle cooling.h number density non-dim, dim, scale factor %.5e, %.5e, %.5e\n", chemistry_get_number_density(p,chemistry_element_H,internal_const), inn_h, cooling->number_density_scale);
+ //printf("Eagle cooling.h number density non-dim, dim, scale factor %.5e, %.5e, %.5e\n", chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const), inn_h, cooling->number_density_scale);
- get_redshift_index(z,&z_index,&dz,cooling); // ADD OFFSET (d_z) CALCULATION INTO THIS FUNCTION
- get_index_1d(cooling->Redshifts, cooling->N_Redshifts, z, &z_index, &d_z);
+ get_redshift_index(z,&z_index,&dz,cooling);
get_index_1d(cooling->HeFrac, cooling->N_He, inHe, &He_i, &d_He);
get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
- get_index_1d_therm(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
-
- //printf("Eagle cooling.h actual z, HeFrac, nH, internal energy %.5e, %.5e, %.5e, %.5e \n", z, inHe, log10(inn_h), log10(u));
- //printf("Eagle cooling.h interp z, HeFrac, nH, internal energy %.5e, %.5e, %.5e, %.5e \n", cooling->Redshifts[z_index], cooling->HeFrac[He_i], cooling->nH[n_h_i], cooling->Therm[u_i]);
-
- if (z_index == cooling->N_Redshifts+1){
- logT = interpol_4d(cooling->table.photoionisation_cooling.temperature, 0, He_i, n_h_i,
- u_i, d_z, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- } else if (z_index > cooling->N_Redshifts+1){
- logT = interpol_4d(cooling->table.collisional_cooling.temperature, 0, He_i, n_h_i,
- u_i, d_z, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- } else {
+ get_index_1d(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
+
+ printf("Eagle cooling.h actual z, HeFrac, nH, internal energy %.5e, %.5e, %.5e, %.5e \n", z, inHe, log10(inn_h), log10(u));
+ printf("Eagle cooling.h interp z, HeFrac, nH, internal energy %.5e, %.5e, %.5e, %.5e, %.5e, %d \n", cooling->Redshifts[z_index], cooling->HeFrac[He_i], cooling->nH[n_h_i], cooling->Therm[u_i], cooling->Redshifts[z_index+1], z_index);
+
+ //if (z_index == cooling->N_Redshifts+1){
+ //logT = interpol_4d(cooling->table.photoionisation_cooling.temperature, 0, He_i, n_h_i,
+ // u_i, d_z, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+ //} else if (z_index > cooling->N_Redshifts+1){
+ //logT = interpol_4d(cooling->table.collisional_cooling.temperature, 0, He_i, n_h_i,
+ // u_i, d_z, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+ //} else {
+ //logT = interpol_4d(cooling->table.element_cooling.temperature, 0, He_i, n_h_i,
+ // u_i, d_z, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+ //}
+
+ //logT = interpol_4d(cooling->table.element_cooling.temperature, z_index, He_i, n_h_i,
+ // u_i, dz, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
logT = interpol_4d(cooling->table.element_cooling.temperature, 0, He_i, n_h_i,
- u_i, d_z, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- }
+ u_i, dz, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+
T = pow(10.0, logT);
- //T = 1.0e6;
if (u_i == 0 && d_u == 0) T *= u / pow(10.0, cooling->Therm[0]);
@@ -422,56 +560,75 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(cons
*/
__attribute__((always_inline)) INLINE static double eagle_cooling_rate(const struct part* restrict p, const struct cooling_function_data* restrict cooling, const struct cosmology* restrict cosmo, const struct phys_const *internal_const) {
double T_gam, LambdaNet = 0.0, solar_electron_abundance;
- double n_h = chemistry_get_number_density(p,chemistry_element_H,internal_const)*cooling->number_density_scale; // chemistry_data
+ double n_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const)*cooling->number_density_scale; // chemistry_data
double z = cosmo->z;
- float d_z,dz; // ARE THESE ACTUALLY MEANT TO BE THE SAME VALUE???
+ float dz; // ARE THESE ACTUALLY MEANT TO BE THE SAME VALUE???
int z_index;
float h_plus_he_electron_abundance;
- float *cooling_solar_abundances = malloc(cooling->N_Elements*sizeof(float));
+ double *cooling_solar_abundances = malloc(cooling->N_Elements*sizeof(double));
int i;
double temp;
int n_h_i, He_i, temp_i;
float d_n_h, d_He, d_temp;
float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He]/(p->chemistry_data.metal_mass_fraction[chemistry_element_H]+p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- //float HeFrac = 0.248;
- //n_h = pow(10,-4.4);
- //z = 0.0;
- get_redshift_index(z,&z_index,&dz,cooling); // ADD OFFSET (d_z) CALCULATION INTO THIS FUNCTION
- get_index_1d(cooling->Redshifts, cooling->N_Redshifts, z, &z_index, &d_z);
+ get_redshift_index(z,&z_index,&dz,cooling);
temp = eagle_convert_u_to_temp(p,cooling,cosmo,internal_const);
+ printf("Eagle cooling.h temperature cgs, internal energy %.5e %.5e\n",temp,hydro_get_physical_internal_energy(p,cosmo)*cooling->internal_energy_scale);
get_index_1d(cooling->Temp, cooling->N_Temp, log10(temp), &temp_i, &d_temp);
get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
get_index_1d(cooling->nH, cooling->N_nH, log10(n_h), &n_h_i, &d_n_h);
- printf("Eagle cooling.h actual z, HeFrac, log nH, log temperature %.5e, %.5e, %.5e, %.5e \n", z, HeFrac, log10(n_h), log10(temp));
- printf("Eagle cooling.h interp z, HeFrac, log nH, log temperature %.5e, %.5e, %.5e, %.5e \n", cooling->Redshifts[z_index], cooling->HeFrac[He_i], cooling->nH[n_h_i], cooling->Temp[temp_i]);
-
+ //printf("Eagle cooling.h actual z, HeFrac, log nH, log temperature %.5e, %.5e, %.5e, %.5e \n", z, HeFrac, log10(n_h), log10(temp));
+ //printf("Eagle cooling.h interp z, HeFrac, log nH, log temperature %.5e, %.5e, %.5e, %.5e \n", cooling->Redshifts[z_index], cooling->HeFrac[He_i], cooling->nH[n_h_i], cooling->Temp[temp_i]);
+
+ char output_filename[21];
+ FILE** output_file = malloc(9*sizeof(FILE*));
+ for (int element = 0; element < cooling->N_Elements; element++){
+ sprintf(output_filename, "%s%d%s", "cooling_output_",element,".dat");
+ output_file[element] = fopen(output_filename, "a");
+ if (output_file == NULL)
+ {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ }
+
+ char *output_filename3 = "iron_output.dat";
+ FILE *output_file3 = fopen(output_filename3, "a");
+ if (output_file3 == NULL)
+ {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+
/* ------------------ */
/* Metal-free cooling */
/* ------------------ */
- if (z_index == cooling->N_Redshifts+1){
- LambdaNet =
- interpol_4d(cooling->table.photoionisation_cooling.H_plus_He_heating, 0, He_i, n_h_i,
- temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- printf("Eagle cooling.h 1 LambdaNet = %.5e\n", LambdaNet);
- h_plus_he_electron_abundance =
- interpol_4d(cooling->table.photoionisation_cooling.H_plus_He_electron_abundance, 0, He_i, n_h_i,
- temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- } else if (z_index > cooling->N_Redshifts+1){
+ //LambdaNet =
+ // interpol_4d(cooling->table.collisional_cooling.H_plus_He_heating, 0, He_i, n_h_i,
+ // temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+ //h_plus_he_electron_abundance =
+ // interpol_4d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, 0, He_i, n_h_i,
+ // temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
LambdaNet =
interpol_4d(cooling->table.collisional_cooling.H_plus_He_heating, 0, He_i, n_h_i,
- temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- printf("Eagle cooling.h 2 LambdaNet = %.5e\n", LambdaNet);
+ temp_i, dz, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
h_plus_he_electron_abundance =
interpol_4d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, 0, He_i, n_h_i,
- temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
- }
+ temp_i, dz, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+
+ //LambdaNet =
+ // interpol_4d(cooling->table.photoionisation_cooling.H_plus_He_heating, 0, He_i, n_h_i,
+ // temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+ //h_plus_he_electron_abundance =
+ // interpol_4d(cooling->table.photoionisation_cooling.H_plus_He_electron_abundance, 0, He_i, n_h_i,
+ // temp_i, 0, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
/* ------------------ */
/* Compton cooling */
@@ -487,52 +644,39 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
LambdaNet -= eagle_compton_rate * (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
h_plus_he_electron_abundance / n_h; // WATCH OUT WHERE h_plus_he_electron_abundance GETS DEFINED!!!
- printf("Eagle cooling.h 3 LambdaNet = %.5e, %.5e, %.5e, %.5e, %.5e, %.5e, %.5e\n", LambdaNet, eagle_compton_rate, temp, eagle_cmb_temperature*(1+z), pow((1+z),4), h_plus_he_electron_abundance, n_h);
+ //printf("Eagle cooling.h 3 LambdaNet = %.5e, %.5e, %.5e, %.5e, %.5e, %.5e, %.5e\n", LambdaNet, eagle_compton_rate, temp, eagle_cmb_temperature*(1+z), pow((1+z),4), h_plus_he_electron_abundance, n_h);
}
/* ------------- */
/* Metal cooling */
/* ------------- */
- if (eagle_metal_cooling_on) {
/* for each element, find the abundance and multiply it
by the interpolated heating-cooling */
set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling_solar_abundances, cooling);
solar_electron_abundance =
- interpol_3d(cooling->table.element_cooling.electron_abundance, 0, n_h_i, temp_i, d_z,
+ interpol_3d(cooling->table.element_cooling.electron_abundance, 0, n_h_i, temp_i, dz,
d_n_h, d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); /* ne/n_h */
- //printf("Eagle cooling.h solar electron abundance %.5e, %.5e\n",solar_electron_abundance, cooling->table.element_cooling.electron_abundance[n_h_i*cooling->N_Temp + temp_i]);
- //for(int ii = 0; ii < cooling->N_Redshifts; ii++){
- // for(int ie = 0; ie < cooling->N_Elements; ie++){
- // for(int ij = 0; ij < cooling->N_nH; ij++){
- // for(int ik = 0; ik < cooling->N_Temp; ik++){
- // double table_metal_heating = cooling->table.element_cooling.metal_heating[ii*cooling->N_Elements*cooling->N_nH*cooling->N_Temp + ie*cooling->N_nH*cooling->N_Temp + ij*cooling->N_Temp + ik];
- // if (table_metal_heating != 0) printf("Eagle cooling.h cooling->table.element_cooling.metal_heating(z,elem,nh,temp),iz,elem,i_h,i_temp, %.5e, %.5e, %d, %.5e, %.5e\n",table_metal_heating,cooling->Redshifts[ii],ie,cooling->nH[ij],cooling->Temp[ik]);
- // //printf("Eagle cooling.h cooling->table.element_cooling.metal_heating(iz,i_h,i_temp),%.5e\n",table_metal_heating);
- // }
- // }
- // }
- //}
+ double element_cooling_lambda;
for (i = 0; i < cooling->N_Elements; i++){
- if (cooling_solar_abundances[i] > 0 && h_plus_he_electron_abundance != 0 && solar_electron_abundance != 0){
- LambdaNet +=
- interpol_4d(cooling->table.element_cooling.metal_heating, 0, i, n_h_i,
- temp_i, d_z, 0.0, d_n_h, d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp) *
+ element_cooling_lambda = interpol_4d(cooling->table.element_cooling.metal_heating, 0, i, n_h_i,
+ temp_i, dz, 0.0, d_n_h, d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp) *
(h_plus_he_electron_abundance / solar_electron_abundance) *
cooling_solar_abundances[i];
- printf("Eagle cooling.h 4 LambdaNet = %.5e\n", LambdaNet);
- } else if (cooling_solar_abundances[i] > 0 && (h_plus_he_electron_abundance == 0 || solar_electron_abundance == 0)){
- LambdaNet +=
- interpol_4d(cooling->table.element_cooling.metal_heating, 0, i, n_h_i,
- temp_i, d_z, 0.0, d_n_h, d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp) *
- cooling_solar_abundances[i];
- printf("Eagle cooling.h 5 LambdaNet = %.5e\n", LambdaNet);
- }
+ LambdaNet += element_cooling_lambda;
+ fprintf(output_file[i],"%.5e\n",element_cooling_lambda);
}
- }
+
+
+ for (i = 0; i < cooling->N_Elements; i++) fclose(output_file[i]);
+
+ float iron_cooling_lambda = interpol_4d(cooling->table.element_cooling.metal_heating, 0, i, n_h_i,
+ temp_i, dz, 0.0, d_n_h, d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp);
+ fprintf(output_file3,"Eagle cooling.h iron interpol, electron abundance ratio, solar abundance %.5e %.5e %.5e\n",iron_cooling_lambda,(h_plus_he_electron_abundance / solar_electron_abundance),
+ cooling_solar_abundances[i]);
return LambdaNet;
}
@@ -555,8 +699,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
const struct cooling_function_data* restrict cooling,
struct part* restrict p, struct xpart* restrict xp, float dt) {
- double u_old = hydro_get_comoving_internal_energy(p)*cooling->internal_energy_scale;
- //double rho = hydro_get_comoving_density(p);
+ double u_old = hydro_get_physical_internal_energy(p,cosmo)*cooling->internal_energy_scale;
float dz;
int z_index;
get_redshift_index(cosmo->z,&z_index,&dz,cooling);
@@ -579,7 +722,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* convert Hydrogen mass fraction in Hydrogen number density */
//inn_h = rho * XH / PROTONMASS;
- inn_h = chemistry_get_number_density(p,chemistry_element_H,phys_const)*cooling->number_density_scale;
+ inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,phys_const)*cooling->number_density_scale;
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
ratefact = inn_h * (XH / eagle_proton_mass_cgs);
@@ -680,11 +823,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* Update the internal energy time derivative */
hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
- //hydro_set_internal_energy_dt(p, hydro_du_dt);
/* Store the radiated energy */
xp->cooling_data.radiated_energy += -hydro_get_mass(p) * cooling_du_dt * dt;
- //xp->cooling_data.radiated_energy += 0.0;
}
@@ -766,9 +907,11 @@ static INLINE void cooling_init_backend(
char fname[200];
parser_get_param_string(parameter_file, "EagleCooling:filename",cooling->cooling_table_path);
+ cooling->reionisation_redshift = parser_get_param_float(parameter_file, "EagleCooling:reionisation_redshift");
GetCoolingRedshifts(cooling);
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
ReadCoolingHeader(fname,cooling);
+ MakeNamePointers(cooling);
cooling->table = eagle_readtable(cooling->cooling_table_path,cooling);
printf("Eagle cooling.h read table \n");
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index a793cec2255881ed69ddce3f88cac2759c1746b1..e9f488549bb0c3c45c5775fdebbc6f8ad9c3f506 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -99,6 +99,8 @@ struct cooling_function_data {
float *SolarElectronAbundance;
char **ElementNames;
char **SolarAbundanceNames;
+ int *ElementNamePointers;
+ int *SolarAbundanceNamePointers;
/*! Constant multiplication factor for time-step criterion */
float cooling_tstep_mult;
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index 68f32c52ffcce411aff4b50bf0820c19965844c4..1af3621875ae9cf2c9bbe777f97b6eb34edcca0b 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/src/cooling/EAGLE/eagle_cool_tables.h
@@ -604,6 +604,17 @@ inline struct eagle_cooling_table eagle_readtable(char *cooling_table_path, cons
return table;
}
+inline int element_index(char *element_name, const struct cooling_function_data* restrict cooling) {
+ int i;
+
+ for (i = 0; i < cooling->N_Elements; i++)
+ if (strcmp(cooling->ElementNames[i], element_name) == 0) return i;
+
+ /* element not found */
+ return -1;
+}
+
+
inline int get_element_index(char *table[20], int size, char *element_name) {
int i;
@@ -614,19 +625,153 @@ inline int get_element_index(char *table[20], int size, char *element_name) {
return -1;
}
-inline int set_cooling_SolarAbundances(const float *element_abundance,
- float *cooling_element_abundance,
- const struct cooling_function_data* restrict cooling) {
- int i, index;
-
- for(i = 0; i < cooling->N_Elements; i++){
- index = get_element_index(cooling->SolarAbundanceNames,
- cooling->N_SolarAbundances, cooling->ElementNames[i]);
- if (cooling->SolarAbundances[index] != 0.0) cooling_element_abundance[i] = element_abundance[i]/cooling->SolarAbundances[index];
- else cooling_element_abundance[i] = 0.0;
+inline void MakeNamePointers(struct cooling_function_data* cooling) {
+ int i, j, sili_index = 0;
+ //char *ElementNames;
+ char ElementNames[cooling->N_Elements][eagle_element_name_length];
+
+ /* This is ridiculous, way too many element name arrays. Needs to be changed */
+ //ElementNames = malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
+ strcpy(ElementNames[0], "Hydrogen");
+ strcpy(ElementNames[1], "Helium");
+ strcpy(ElementNames[2], "Carbon");
+ strcpy(ElementNames[3], "Nitrogen");
+ strcpy(ElementNames[4], "Oxygen");
+ strcpy(ElementNames[5], "Neon");
+ strcpy(ElementNames[6], "Magnesium");
+ strcpy(ElementNames[7], "Silicon");
+ strcpy(ElementNames[8], "Iron");
+
+ cooling->ElementNamePointers = malloc(cooling->N_Elements * sizeof(int));
+ cooling->SolarAbundanceNamePointers = malloc(cooling->N_Elements * sizeof(int));
+
+ for (i = 0; i < cooling->N_Elements; i++) {
+ if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
}
- return 0;
+ for (i = 0; i < cooling->N_Elements; i++) {
+ cooling->SolarAbundanceNamePointers[i] = -999;
+ cooling->ElementNamePointers[i] = -999;
+
+ for (j = 0; j < cooling->N_SolarAbundances; j++) {
+ if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) == 0)
+ cooling->SolarAbundanceNamePointers[i] = j;
+ }
+
+ if (strcmp(cooling->ElementNames[i], "Sulphur") == 0 ||
+ strcmp(cooling->ElementNames[i], "Calcium") ==
+ 0) /* These elements are tracked! */
+ cooling->ElementNamePointers[i] = -1 * sili_index;
+ else {
+ for (j = 0; j < cooling->N_Elements; j++) {
+ if (strcmp(cooling->ElementNames[i], ElementNames[j]) == 0)
+ cooling->ElementNamePointers[i] = j;
+ }
+ }
+ }
}
+
+//inline int set_cooling_SolarAbundances(const float *element_abundance,
+// float *cooling_element_abundance,
+// const struct cooling_function_data* restrict cooling) {
+// int i, index;
+//
+// for(i = 0; i < cooling->N_Elements; i++){
+// index = get_element_index(cooling->SolarAbundanceNames,
+// cooling->N_SolarAbundances, cooling->ElementNames[i]);
+// if (cooling->SolarAbundances[index] != 0.0) cooling_element_abundance[i] = element_abundance[i]/cooling->SolarAbundances[index];
+// else cooling_element_abundance[i] = 0.0;
+// printf ("eagle_cool_tables.h element, name, abundance, solar abundance, cooling abundance %d %s %.5e %.5e %.5e\n",index,cooling->ElementNames[i], element_abundance[i],cooling->SolarAbundances[index], cooling_element_abundance[i]);
+// }
+//
+// return 0;
+//}
+
+//inline int set_cooling_SolarAbundances(const float *element_abundance,
+// float *cooling_element_abundance,
+// const struct cooling_function_data* restrict cooling) {
+// int i, index;
+//
+// int static Silicon_SPH_Index = -1;
+// int static Calcium_SPH_Index = -1;
+// int static Sulphur_SPH_Index = -1;
+//
+// int static Silicon_CoolHeat_Index = -1;
+// int static Calcium_CoolHeat_Index = -1;
+// int static Sulphur_CoolHeat_Index = -1;
+//
+// static int first_call = 0;
+//
+// if (first_call == 0) {
+// /* determine (inverse of) solar abundance of these elements */
+// for (i = 0; i < cooling_N_Elements; i++) {
+// index =
+// get_element_index(cooling_SolarAbundanceNames,
+// cooling_N_SolarAbundances, cooling_ElementNames[i]);
+//
+// if (index < 0) endrun(-12345);
+//
+// index = SolarAbundanceNamePointers[i];
+//
+// cooling_ElementAbundance_SOLARM1[i] = 1. / cooling_SolarAbundances[index];
+//
+// index = ElementNamePointers[i];
+//
+// if (index < 0 && ThisTask == 0)
+// printf("[bg_cooling] element not found %s\n", cooling_ElementNames[i]);
+// }
+//
+// /* Sulphur tracks Silicon: may choose not to follow Sulphur as SPH element
+// */
+// /* Same is true for Calcium */
+// /* We will assume the code tracks Silicon, and may need to scale Calcium and
+// * Sulphur accordingly */
+//
+// Silicon_SPH_Index = element_index("Silicon");
+// Calcium_SPH_Index = element_index("Calcium");
+// Sulphur_SPH_Index = element_index("Sulphur");
+//
+// Silicon_CoolHeat_Index =
+// get_element_index(cooling_ElementNames, cooling_N_Elements, "Silicon");
+// Calcium_CoolHeat_Index =
+// get_element_index(cooling_ElementNames, cooling_N_Elements, "Calcium");
+// Sulphur_CoolHeat_Index =
+// get_element_index(cooling_ElementNames, cooling_N_Elements, "Sulphur");
+//
+// if (Silicon_CoolHeat_Index == -1 || Calcium_CoolHeat_Index == -1 ||
+// Sulphur_CoolHeat_Index == -1) {
+// if (ThisTask == 0)
+// printf("[bg_cooling] error: did not find Si or Ca or S??\n");
+// endrun(-1233);
+// }
+//
+// first_call = 1;
+// }
+// for (i = 0; i < cooling_N_Elements; i++) {
+// if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1)
+// /* SPH does not track Calcium: use Si abundance */
+// if (Silicon_SPH_Index == -1)
+// cooling_element_abundance[i] = 0.0;
+// else
+// cooling_element_abundance[i] =
+// element_abundance[Silicon_SPH_Index] *
+// cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
+// else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index == -1)
+// /* SPH does not track Sulphur: use Si abundance */
+// if (Silicon_SPH_Index == -1)
+// cooling_element_abundance[i] = 0.0;
+// else
+// cooling_element_abundance[i] =
+// element_abundance[Silicon_SPH_Index] *
+// cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
+// else
+// cooling_element_abundance[i] = element_abundance[ElementNamePointers[i]] *
+// cooling_ElementAbundance_SOLARM1[i];
+// }
+//
+// return 0;
+//}
+
+
#endif
diff --git a/tests/plots.py b/tests/plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..844dea6d1745b99650b96d7a3a8d4de1a4936665
--- /dev/null
+++ b/tests/plots.py
@@ -0,0 +1,38 @@
+import matplotlib.pyplot as plt
+
+elements = 9
+temperature = []
+cooling_rate = [[] for i in range(elements+1)]
+length = 0
+file_in = open('cooling_output.dat', 'r')
+for line in file_in:
+ data = line.split()
+ temperature.append(float(data[0]))
+ cooling_rate[0].append(-float(data[1]))
+
+file_in.close()
+
+for elem in range(elements):
+ file_in = open('cooling_output_'+str(elem)+'.dat','r')
+ for line in file_in:
+ data = line.split()
+ cooling_rate[elem+1].append(-float(data[0]))
+ file_in.close()
+
+
+p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'b', label = 'Carbon')
+p2, = plt.loglog(temperature, cooling_rate[2], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+p3, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'r', label = 'Oxygen')
+p4, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'c', label = 'Neon')
+p5, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'm', label = 'Magnesium')
+p6, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'y', label = 'Silicon')
+p7, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+p8, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'olive', label = 'Calcium')
+p9, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
+plt.ylim([1e-24,1e-21])
+plt.xlim([1e4,1e8])
+plt.xlabel('Temperature (K)')
+plt.ylabel('Cooling rate (eagle units...)')
+plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9])
+plt.show()
diff --git a/tests/testCooling.c b/tests/testCooling.c
index 420a5ad5a5c2602bec4189ab5e05279713d54590..33f6989d563bcb2dfe480986ba14bc33d61c9038 100644
--- a/tests/testCooling.c
+++ b/tests/testCooling.c
@@ -30,6 +30,7 @@ int main() {
struct unit_system us;
struct chemistry_data chemistry_data;
struct part p;
+ struct xpart xp;
struct phys_const internal_const;
struct cooling_function_data cooling;
struct cosmology cosmo;
@@ -50,43 +51,60 @@ int main() {
units_init(&us, params, "InternalUnitSystem");
phys_const_init(&us, params, &internal_const);
- double number_density_cgs = 0.1;
- double temperature_cgs = 1.0e6;
- //double power_scale = units_cgs_conversion_factor(&us,UNIT_CONV_POWER)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS);
- double power_per_num_density_factor = units_cgs_conversion_factor(&us,UNIT_CONV_POWER)*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3)/number_density_cgs;
-
- double gamma = 5.0/3.0;
-
- double number_density = number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
- p.rho = number_density*(1.0/0.6*internal_const.const_proton_mass);
-
- for (int i = 0; i < 9; i++) p.chemistry_data.metal_mass_fraction[i] = 0.0;
- p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen] = 0.752;
- p.chemistry_data.metal_mass_fraction[EAGLE_Helium] = 0.248;
-
- double temperature = temperature_cgs/units_cgs_conversion_factor(&us,UNIT_CONV_TEMPERATURE);
- double pressure = number_density*internal_const.const_boltzmann_k*temperature;
- printf("non-dim number density, code number density, number density scale %.5e, %.5e, %.5e\n",number_density, chemistry_get_number_density(&p,chemistry_element_H,&internal_const), pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3));
- printf("number density, boltzmann constant, temperature: %.5e, %.5e, %.5e\n",number_density_cgs, internal_const.const_boltzmann_k, temperature);
- printf("proton mass, boltzmann constant, %.5e, %.5e\n", internal_const.const_proton_mass, internal_const.const_boltzmann_k);
- double internal_energy = pressure/(p.rho*(gamma - 1.0));
- //p.entropy = internal_energy/((gamma - 1.0)*pow(p.rho,gamma - 1.0));
- p.entropy = pressure/(pow(p.rho,gamma));
- printf("double check pressure, actual pressure: %.5e, %.5e\n", hydro_get_comoving_pressure(&p), pressure);
- printf("temperature, pressure, internal energy, entropy: %.5e, %.5e, %.5e, %.5e, %.5e\n", temperature,pressure,internal_energy,p.entropy,internal_const.const_boltzmann_k);
+ double hydrogen_number_density_cgs = 1e-4;
+ double u, hydrogen_number_density, pressure, gamma, cooling_du_dt, temperature_cgs;
+ int n_t_i = 2000;
+
+ char *output_filename = "cooling_output.dat";
+ FILE *output_file = fopen(output_filename, "w");
+ if (output_file == NULL)
+ {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ char *output_filename2 = "temperature_output.dat";
+ FILE *output_file2 = fopen(output_filename2, "w");
+ if (output_file2 == NULL)
+ {
+ printf("Error opening file!\n");
+ exit(1);
+ }
cosmology_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo);
- printf("testCooling.c redshift %.5e\n", cosmo.z);
cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&cooling);
+ //for(int j = 0; j < cooling.N_Redshifts; j++) printf("redshift %.5e\n",cooling.Redshifts[j]);
chemistry_init(params, &us, &internal_const, &chemistry_data);
+ chemistry_first_init_part(&p,&xp,&chemistry_data);
chemistry_print(&chemistry_data);
-
- double cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const)*power_per_num_density_factor;
- printf("cooling rate: %.5e\n",cooling_du_dt);
+ for (int k = 0; k < cooling.N_Elements; k++){
+ printf("element, abundances %d, %.5e\n", k, p.chemistry_data.metal_mass_fraction[k]);
+ }
+
+ for(int t_i = 0; t_i < n_t_i; t_i++){
+ //for (int element = 0; element < cooling.N_Elements; element++){
+ // printf("element %d abundance %.5e particle abundance %.5e\n",element,chemistry_data.initial_metal_mass_fraction[element],p.chemistry_data.metal_mass_fraction[element]);
+ //}
+
+ u = 1.0*pow(10.0,11 + t_i*6.0/n_t_i)/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
+ pressure = u*p.rho*(gamma -1.0);
+ hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
+ p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ p.entropy = pressure/(pow(p.rho,gamma));
+
+ gamma = 5.0/3.0;
+
+ cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const);
+ temperature_cgs = eagle_convert_u_to_temp(&p,&cooling,&cosmo,&internal_const);
+ //fprintf(output_file,"%.5e %.5e\n",u*units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS),cooling_du_dt);
+ fprintf(output_file,"%.5e %.5e\n",temperature_cgs,cooling_du_dt);
+ fprintf(output_file2,"%.5e %.5e\n",u*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)), temperature_cgs);
+ }
+ fclose(output_file);
+ fclose(output_file2);
free(params);