39 double edensave =
dense.eden;
43 if (
dense.lgElmtOn[nelem])
46 dense.SetGasPhaseDensity( nelem,
dense.gas_phase[nelem] * factor);
54 fprintf(
ioQQQ,
" EDEN change PressureChange from to %10.3e %10.3e %10.3e\n",
58 hmi.H2_total *= factor;
59 h2.ortho_density *= factor;
60 h2.para_density *= factor;
61 for(
long mol=0; mol <
mole_global.num_total; mol++ )
63 mole.species[mol].den *= factor;
65 dense.updateXMolecules();
73 for(
long ion=0; ion < (nelem + 2); ion++ )
75 dense.xIonDense[nelem][ion] *= factor;
76 if (nelem-ion >= 0 && nelem-ion <
NISO)
104 for( ChemAtomList::iterator atom =
atom_list.begin(); atom !=
atom_list.end(); ++atom )
106 long nelem = (*atom)->el->Z-1;
108 if(
dense.lgElmtOn[nelem] )
112 double sum_monatomic = 0.;
113 for(
long ion=0; ion<nelem+2; ++ion )
115 sum_monatomic +=
dense.xIonDense[nelem][ion] * (*atom)->frac;
117 realnum sum_molecules =
dense.xMolecules(nelem) * (*atom)->frac;
118 realnum sum_gas_phase =
dense.gas_phase[nelem] * (*atom)->frac;
119 if( sum_monatomic + sum_molecules <=
SMALLFLOAT &&
122 fprintf(
ioQQQ,
"PROBLEM non-conservation of nuclei %s\tions %g moles %g error %g of %g\n",
123 (*atom)->label().c_str(),
126 sum_monatomic + sum_molecules-sum_gas_phase,
135 if( fabs( sum_monatomic + sum_molecules - sum_gas_phase ) >
136 conv.GasPhaseAbundErrorAllowed * sum_gas_phase )
138 fprintf(
ioQQQ,
"PROBLEM non-conservation of nuclei %s\t nzone %li atoms %.12e moles %.12e sum %.12e tot gas %.12e rel err %.3e\n",
139 (*atom)->label().c_str(),
143 sum_monatomic + sum_molecules,
145 (sum_monatomic + sum_molecules - sum_gas_phase)/sum_gas_phase );
162 if( 0 &&
conv.lgConvIoniz() &&
163 dense.lgElmtOn[nelem] &&
164 dense.xIonDense[nelem][ionStage]/
dense.eden >
conv.EdenErrorAllowed/10. &&
168 for (
long ipLev=0; ipLev<numStates; ipLev++)
170 abund += states[ipLev].Pop();
173 double rel_err = fabs(
abund-
dense.xIonDense[nelem][ionStage])/
178 if( rel_err > err_tol )
181 if( 0 &&
nzone > 0 && loop_ion > 0 )
183 fprintf(
ioQQQ,
"PROBLEM Inconsistent states/stage pops nzone %3ld loop_ion %2ld nelem %2ld ion %2ld states = %e stage = %e error = %e\n",
189 dense.xIonDense[nelem][ionStage],
193 sprintf( chConvIoniz ,
"States!=stage pops nelem %ld ion %ld ", nelem, ionStage );
194 conv.setConvIonizFail( chConvIoniz,
abund,
195 dense.xIonDense[nelem][ionStage] );
209 if(
dense.lgElmtOn[nelem] )
211 for(
long ion=0; ion<=nelem+1; ++ion )
212 DenseAtomsIons +=
dense.xIonDense[nelem][ion];
219 ASSERT( DenseAtomsIons > 0. );
227 dense.xNucleiTotal = den_mole + DenseAtomsIons;
230 fprintf(
ioQQQ,
"PROBLEM DISASTER SumDensities has found "
231 "dense.xNucleiTotal with an insane density.\n");
232 fprintf(
ioQQQ,
"The density was %.2e\n",
244 for(
long i=0; i <
LIMELM; i++ )
263 if(
dense.xMassDensity0 < 0.0 )
284 bool lgChange =
false;
291 for(
long nelem=1; nelem <
LIMELM; nelem++ )
293 if(
abund.lgAbunTabl[nelem] )
298 double hold = abun/
dense.gas_phase[nelem];
301 for(
long ion=0; ion < (nelem + 2); ion++ )
313 if( !
dense.lgDenFlucOn )
315 static double FacAbunSav;
316 double OldAbun = 0.0;
321 OldAbun = FacAbunSav;
325 if(
dense.lgDenFlucRadius )
340 FacAbun = FacAbunSav/OldAbun;
351 if (
dense.lgElmtOn[nelem])
358 for(
long mol=0; mol <
mole_global.num_total; mol++ )
375 const int SCALE_NEW = 1;
390 return struc.hden[i];
double AbundancesTable(double r0, double depth, long int iel)
const int INPUT_LINE_LENGTH
NORETURN void TotalInsanity(void)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
void EdenChange(double EdenNew)
UNUSED const realnum BIGFLOAT
void ScaleAllDensities(realnum factor)
bool lgElemsConserved(void)
realnum scalingDensity(void)
void ScaleIonDensities(const long nelem, const realnum factor)
void lgStatesConserved(long nelem, long ionStage, qList states, long numStates, realnum err_tol, long loop_ion)
realnum scalingZoneDensity(long i)
void SetGasPhaseDeuterium(const realnum &Hdensity)
void ScaleDensitiesDeuterium(const realnum &factor)
t_elementnames elementnames
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
void iso_renorm(long nelem, long ipISO, double &renorm)
t_mole_global mole_global
realnum total_molecules_gasphase(void)
void mole_print_species_reactions(molecule *speciesToPrint)
molecule * findspecies(const char buf[])
void total_molecule_elems(realnum total[LIMELM])
UNUSED const double ATOMIC_MASS_UNIT
realnum m_xMolecules[LIMELM]
realnum gas_phase[LIMELM]
void SetGasPhaseDensity(const long nelem, const realnum density)
void TempChange(double TempNew, bool lgForceUpdate)