corrected lack of T^3/2 term in q_recomb

jhmatthews · Nov 26, 2015 · 3c316af · 3c316af
1 parent 4480157
commit 3c316af
Show file tree

Hide file tree

Showing 3 changed files with 13 additions and 5 deletions.
diff --git a/source/direct_ion.c b/source/direct_ion.c
@@ -337,6 +337,7 @@ q_recomb_dere (cont_ptr, electron_temperature)
   int nion;
   double u0;
   double gaunt, coeff;
+  double root_etemp;
 
   gaunt = 0.1;      //for now - from Mihalas for hydrogen
   u0 = cont_ptr->freq[0] * H_OVER_K / electron_temperature;
@@ -345,10 +346,13 @@ q_recomb_dere (cont_ptr, electron_temperature)
      only do this if it is the ground state */
   if (ion[nion].dere_di_flag == 1 && config[cont_ptr->nlev].ilv == 1)
     {
-      coeff = 2.07e-16 * config[cont_ptr->nlev].g / config[cont_ptr->uplev].g;
+      root_etemp = sqrt(electron_temperature);     
+      coeff = 2.07e-16 / (root_etemp * root_etemp * root_etemp);
+      coeff *= config[cont_ptr->nlev].g / config[cont_ptr->uplev].g;
       coeff *= exp(u0);
       coeff *= q_ioniz_dere(nion, electron_temperature);
 
+
       /* do a sane check here, as there's an exponential which could blow up */
       if (sane_check(coeff))
         Error("q_recomb is %8.4e for ion %i at temperature %8.4e\n",

diff --git a/source/para_update.c b/source/para_update.c
@@ -50,9 +50,9 @@ communicate_estimators_para ()
      for each cell, plus three arrays, each of length NXBANDS */
   plasma_double_helpers = (15 + 3 * NXBANDS) * NPLASMA;
 
-  /* The size of the helper array for integers. We transmit 6 numbers 
+  /* The size of the helper array for integers. We transmit 7 numbers 
      for each cell, plus one array of length NXBANDS */
-  plasma_int_helpers = (6 + NXBANDS) * NPLASMA;
+  plasma_int_helpers = (7 + NXBANDS) * NPLASMA;
 
 
   maxfreqhelper = calloc (sizeof (double), NPLASMA);
@@ -163,9 +163,11 @@ communicate_estimators_para ()
       iredhelper[mpi_i + 3 * NPLASMA] = plasmamain[mpi_i].ntot_disk;
       iredhelper[mpi_i + 4 * NPLASMA] = plasmamain[mpi_i].ntot_wind;
       iredhelper[mpi_i + 5 * NPLASMA] = plasmamain[mpi_i].ntot_agn;
+      iredhelper[mpi_i + 6 * NPLASMA] = plasmamain[mpi_i].nioniz;
+
       for (mpi_j = 0; mpi_j < NXBANDS; mpi_j++)
 	{
-	  iredhelper[mpi_i + (6 + mpi_j) * NPLASMA] = plasmamain[mpi_i].nxtot[mpi_j];
+	  iredhelper[mpi_i + (7 + mpi_j) * NPLASMA] = plasmamain[mpi_i].nxtot[mpi_j];
 	}
     }
   MPI_Reduce (iredhelper, iredhelper2, plasma_int_helpers, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
@@ -183,9 +185,10 @@ communicate_estimators_para ()
       plasmamain[mpi_i].ntot_disk = iredhelper2[mpi_i + 3 * NPLASMA];
       plasmamain[mpi_i].ntot_wind = iredhelper2[mpi_i + 4 * NPLASMA];
       plasmamain[mpi_i].ntot_agn = iredhelper2[mpi_i + 5 * NPLASMA];
+	  plasmamain[mpi_i].nioniz = iredhelper2[mpi_i + 6 * NPLASMA];
       for (mpi_j = 0; mpi_j < NXBANDS; mpi_j++)
 	{
-	  plasmamain[mpi_i].nxtot[mpi_j] = iredhelper2[mpi_i + (6 + mpi_j) * NPLASMA];
+	  plasmamain[mpi_i].nxtot[mpi_j] = iredhelper2[mpi_i + (7 + mpi_j) * NPLASMA];
 	}
     }
 

diff --git a/source/wind_updates2d.c b/source/wind_updates2d.c
@@ -180,6 +180,7 @@ WindPtr (w);
 
       nwind = plasmamain[n].nwind;
       volume = w[nwind].vol;
+
       //volume = plasmamain[n].vol; 	// 1411 JM -- we want volume to be the filled volume
       //1504 -- JM No we don't! The mean intensity shouldn't change if we *just* clump