37 CollisIonizatCoolingTotal,
38 CollisIonizatHeatingTotal,
39 dCollisIonizatCoolingTotalDT,
44 CollisIonizatCoolingDT,
47 valarray<double> CollisIonizatCoolingArr,
48 CollisIonizatCoolingDTArr,
52 long int nlo_heat_max , nhi_heat_max;
70 if(
dense.xIonDense[nelem][nelem-ipISO]<=0. ||
71 !
dense.lgElmtOn[nelem] )
87 for(
long ipLo=0; ipLo < ipHi; ++ipLo )
110 CollisIonizatCoolingTotal = 0.;
111 CollisIonizatHeatingTotal = 0.;
112 dCollisIonizatCoolingTotalDT = 0.;
118 CollisIonizatCooling =
121 CollisIonizatCoolingTotal += CollisIonizatCooling;
122 CollisIonizatHeating =
125 dense.xIonDense[nelem][nelem+1-ipISO];
126 CollisIonizatHeatingTotal += CollisIonizatHeating;
128 double CollisIonizatDiff = CollisIonizatCooling-CollisIonizatHeating;
130 CollisIonizatCoolingDT = CollisIonizatDiff*
134 dCollisIonizatCoolingTotalDT += CollisIonizatCoolingDT;
138 CollisIonizatCoolingArr[n] = CollisIonizatDiff;
139 CollisIonizatCoolingDTArr[n] = CollisIonizatCoolingDT;
143 double CollisIonizatNetCooling =
144 CollisIonizatCoolingTotal-CollisIonizatHeatingTotal;
147 sp->
coll_ion = CollisIonizatNetCooling;
150 thermal.dCooldT += dCollisIonizatCoolingTotalDT;
153 sprintf(chLabel ,
"IScion%2s%2s" ,
elementnames.chElementSym[ipISO] ,
169 thermal.heating[0][3] += CollisIonizatHeatingTotal;
175 if (CollisIonizatNetCooling < 0)
176 thermal.heating[0][3] -= CollisIonizatNetCooling;
182 fprintf(
ioQQQ,
"DEBUG coll ioniz cool contributors:");
185 if( CollisIonizatCoolingArr[n] /
SDIV( CollisIonizatNetCooling ) > 0.1 )
186 fprintf(
ioQQQ,
" %2li %.1e",
188 CollisIonizatCoolingArr[n]/
SDIV( CollisIonizatNetCooling ) );
191 fprintf(
ioQQQ,
"DEBUG coll ioniz derivcontributors:");
194 if( CollisIonizatCoolingDTArr[n] /
SDIV( dCollisIonizatCoolingTotalDT ) > 0.1 )
195 fprintf(
ioQQQ,
" %2li %.1e",
197 CollisIonizatCoolingDTArr[n]/
SDIV( dCollisIonizatCoolingTotalDT ) );
209 edenIonAbund =
dense.eden*
dense.xIonDense[nelem][nelem+1-ipISO];
212 if(
dense.gas_phase[nelem] > 1e-3 *
dense.xNucleiTotal )
234 enum {DEBUG_LOC=
false};
261 ASSERT( sp->
fb[n].PhotoHeat >= 0. );
264 SavePhotoHeat[n] = sp->
st[n].Pop()*
266 HeatExcited += SavePhotoHeat[n];
267 if( SavePhotoHeat[n] > biggest )
269 biggest = SavePhotoHeat[n];
275 enum {DEBUG_LOC=
false};
276 if( DEBUG_LOC && ipISO==0 && nelem==0 &&
nzone > 700)
281 fprintf(
ioQQQ,
"ipISO = %li nelem=%li ipbig=%li biggest=%g HeatExcited=%.2e ctot=%.2e\n",
290 fprintf(
ioQQQ,
"DEBUG phot heat%2li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
292 SavePhotoHeat[n]/HeatExcited,
293 dense.xIonDense[nelem][nelem+1-ipISO],
316 sprintf(chLabel ,
"ISrcol%2s%2s" ,
elementnames.chElementSym[ipISO] ,
334 SaveInducCool[n] = sp->
fb[n].RecomInducCool_Coef*sp->
fb[n].PopLTE*edenIonAbund;
336 thermal.dCooldT += SaveInducCool[n]*
342 enum {DEBUG_LOC=
false};
343 if( DEBUG_LOC && ipISO==0 && nelem==5 )
345 fprintf(
ioQQQ,
" ipISO=%li nelem=%li ctot = %.2e\n",
349 fprintf(
ioQQQ,
"sum\t%.2e\t%.2e\t%.2e\n",
353 fprintf(
ioQQQ,
"sum\tp ht\tr cl\ti cl\n");
358 fprintf(
ioQQQ,
"%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e \n",
361 sp->
fb[n].RadRecCoolCoef,
363 sp->
fb[n].RecomInducCool_Coef,
365 sp->
fb[n].RecomInducRate);
367 fprintf(
ioQQQ,
" \n");
371 sprintf(chLabel ,
"ISicol%2s%2s" ,
elementnames.chElementSym[ipISO] ,
386 fprintf(
ioQQQ,
"%ld %ld %ld %g\n",nelem,ipISO,ipHi,
387 sp->
st[ipHi].Pop()/sp->
st[ipHi].g());
390 fprintf(
ioQQQ,
"%ld %ld %ld %g\n",nelem,ipISO,ipHi,
391 sp->
st[ipHi].Pop()/(sp->
st[ipHi].Boltzmann()*
400 for (
unsigned lev = 0; lev < Boltzmann.size(); ++lev)
401 Boltzmann[lev] = 1.0;
405 double arg = (sp->
st[ipHi].energy().K()-sp->
st[ipHi-1].energy().K()) /
408 double bstep =
sexp( arg );
409 for (
long ipLo = 0; ipLo < ipHi; ++ipLo )
410 Boltzmann[ipLo] *= bstep;
412 for(
long ipLo=0; ipLo < ipHi; ++ipLo )
419 sp->
st[ipHi].g()/sp->
st[ipLo].g() -
444 if( hlone < heat_max )
450 dCdT_all -= hlone*
thermal.halfte;
456 enum {DEBUG_LOC=
false};
460 fprintf(
ioQQQ,
"%li %li %.2e\n", nlo_heat_max, nhi_heat_max, heat_max );
474 sp->
st[ipHi].Boltzmann()*
475 sp->
st[ipHi].g()/sp->
st[0].g() -
479 if( ipHi ==
iso_ctrl.nLyaLevel[ipISO] )
505 sexp( (sp->
st[ipHi].energy().K()-sp->
st[ipHi-1].energy().K()) /
516 sp->
st[ipHi].g()/sp->
st[ipLo].g() -
519 ASSERT( sp->
st[ipHi].energy().Erg() >
520 sp->
st[ipLo].energy().Erg() );
527 sp->
st[ipHi].g()/sp->
st[ipLo].g() -
530 ASSERT( sp->
st[ipHi].energy().Erg() >
531 sp->
st[ipLo].energy().Erg() );
543 for (
unsigned lev = 0; lev < Boltzmann.size(); ++lev)
549 sexp( (sp->
st[ipHi].energy().K()-sp->
st[ipHi-1].energy().K()) /
551 for (
long ipLo = 0; ipLo < ipHi; ++ipLo )
552 Boltzmann[ipLo] *= bstep;
554 for(
long ipLo=3; ipLo < ipHi; ++ipLo )
560 sp->
st[ipHi].g()/sp->
st[ipLo].g() -
573 sprintf(chLabel ,
"ISclin%2s%2s" ,
elementnames.chElementSym[ipISO] ,
587 sp->
dLTot = -dCdT_all;
592 enum {DEBUG_LOC=
false};
597 fprintf(
ioQQQ,
"%.2e la %.2f restly %.2f barest %.2f hrest %.2f\n",
609 enum {DEBUG_LOC=
false};
612 if( DEBUG_LOC && (nelem==1 || nelem==0) )
620 fprintf(
ioQQQ,
"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
621 sp->
fb[n].gamnc *sp->
fb[n].PopLTE,
626 sp->
fb[n].RecomInducRate*sp->
fb[n].PopLTE ,
628 sp->
fb[n].RecomInducRate*sp->
fb[n].PopLTE ,
632 sp->
fb[n].PhotoHeat*sp->
fb[n].PopLTE,
633 sp->
fb[n].RecomInducCool_Coef*sp->
fb[n].PopLTE+sp->
fb[n].RadRecCoolCoef,
635 sp->
fb[n].RecomInducCool_Coef*sp->
fb[n].PopLTE ,
637 sp->
fb[n].RadRecCoolCoef );
645 fprintf(
ioQQQ,
"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
646 sp->
fb[n].gamnc*sp->
st[n].Pop(),
649 sp->
fb[n].RecomInducRate*sp->
fb[n].PopLTE *edenIonAbund ,
651 sp->
fb[n].RecomInducRate*sp->
fb[n].PopLTE *edenIonAbund ,
654 SaveInducCool[n]+sp->
fb[n].RadRecCoolCoef*edenIonAbund ,
658 sp->
fb[n].RadRecCoolCoef*edenIonAbund );
665 for(
int i = coolnum; i <
thermal.ncltot ; i++ )