diff --git a/source/direct_ion.c b/source/direct_ion.c index 0db39a4d8..58ed1970f 100644 --- 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 index dd2596142..174fcc75d 100644 --- 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 index efdc25b17..156986dd1 100644 --- 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