Commits

Michael Kuhlen committed 04dc362

A few bug fixes. UV background interpolation appears to work.

Comments (0)

Files changed (6)

input/UVB_rates_FG2011.hdf5

Binary file modified.

input/UVB_rates_HM2005.hdf5

Binary file modified.

input/UVB_rates_HM2012.hdf5

Binary file modified.

src/clib/freefall.C

   gr_float gravitational_constant = 4.0 * 3.1415926 * 6.6726e-8 * 
     my_units.density_units * POW(my_units.time_units, 2);
 
-  gr_float a_value = 1.0;
+  // gr_float a_value = 1.0;
+  gr_float a_value = 1.0/(1.0+10.0);
 
   if (initialize_chemistry_data(my_chemistry, my_units, a_value) == FAIL) {
     fprintf(stderr, "Error in initialize_chemistry_data.\n");
     return FAIL;
   }
 
+  my_chemistry.UVbackground = 1;
   my_chemistry.UVbackground_file = (char*) "UVB_rates_HM2012.hdf5";
   if (initialize_UVbackground_data(my_chemistry, my_units, a_value) == FAIL) {
     fprintf(stderr, "Error in initialize_UVbackground_data.\n");
     return FAIL;
   }
 
+  update_UVbackground_rates(my_chemistry, my_units, a_value);
+
   gr_float *density, *energy, *x_velocity, *y_velocity, *z_velocity;
   gr_float *HI_density, *HII_density, *HM_density,
     *HeI_density, *HeII_density, *HeIII_density,

src/clib/grackle.h

 
 int initialize_UVbackground_data(chemistry_data &my_chemistry,
 				 code_units &my_units, gr_float a_value);
+int update_UVbackground_rates(chemistry_data &my_chemistry,
+			      code_units &my_units, float a_value);
 
 int solve_chemistry(chemistry_data &my_chemistry,
 		    code_units &my_units,

src/clib/update_UVbackground_rates.C

   // find interpolation index
   gr_float *zvec = my_chemistry.UVbackground_table.z;
   gr_int index=0;
-  while (Redshift >= zvec[index++]);
+  while (Redshift > zvec[index])
+    index++;
+  if(index == 0) index=1;
   if(index == my_chemistry.UVbackground_table.Nz) index--;
 
+  // printf("index = %d, %.3f <= %.3f <= %.3f\n",index,zvec[index-1],Redshift,zvec[index]);
+
   // *** k24 ***
   gr_float slope = (my_chemistry.UVbackground_table.k24[index] - my_chemistry.UVbackground_table.k24[index-1]) / (zvec[index] - zvec[index-1]);
   my_chemistry.k24 = (Redshift - zvec[index-1]) * slope + my_chemistry.UVbackground_table.k24[index-1];
   my_chemistry.piHeI = (Redshift - zvec[index-1]) * slope + my_chemistry.UVbackground_table.piHeI[index-1];
 
 
+  // printf("%e %e %e\n",my_chemistry.k24,my_chemistry.k25,my_chemistry.k26);
+
 
   // Now apply the Ramp factor
   my_chemistry.k24 *= Ramp;