17 long int nLevelCalled,
98 else if( chExUnits==
'w' )
101 TeInverse = 1./
phycon.te_wn;
107 long int nlev = nLevelCalled;
110 source[nlev-1]==0. && nlev > 1)
119 if(
abund == 0. || nlev==1 )
129 for( level=1; level < nlev; level++ )
141 for( ihi=1; ihi < nlev; ihi++ )
143 for( ilo=0; ilo < ihi; ilo++ )
160 if( lgDeBug || (
trace.lgTrace &&
trace.lgTrLevN) )
162 fprintf(
ioQQQ,
" atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel,
abund);
163 fprintf(
ioQQQ,
" AulDest\n" );
165 fprintf(
ioQQQ,
" hi lo" );
167 for( ilo=0; ilo < nlev-1; ilo++ )
169 fprintf(
ioQQQ,
"%4ld ", ilo );
171 fprintf(
ioQQQ,
" \n" );
173 for( ihi=1; ihi < nlev; ihi++ )
175 fprintf(
ioQQQ,
"%3ld", ihi );
176 for( ilo=0; ilo < ihi; ilo++ )
180 fprintf(
ioQQQ,
"\n" );
183 fprintf(
ioQQQ,
" A*esc\n" );
184 fprintf(
ioQQQ,
" hi lo" );
185 for( ilo=0; ilo < nlev-1; ilo++ )
187 fprintf(
ioQQQ,
"%4ld ", ilo );
189 fprintf(
ioQQQ,
" \n" );
191 for( ihi=1; ihi < nlev; ihi++ )
193 fprintf(
ioQQQ,
"%3ld", ihi );
194 for( ilo=0; ilo < ihi; ilo++ )
198 fprintf(
ioQQQ,
"\n" );
201 fprintf(
ioQQQ,
" AulPump\n" );
203 fprintf(
ioQQQ,
" hi lo" );
204 for( ilo=0; ilo < nlev-1; ilo++ )
206 fprintf(
ioQQQ,
"%4ld ", ilo );
208 fprintf(
ioQQQ,
" \n" );
210 for( ihi=1; ihi < nlev; ihi++ )
212 fprintf(
ioQQQ,
"%3ld", ihi );
213 for( ilo=0; ilo < ihi; ilo++ )
217 fprintf(
ioQQQ,
"\n" );
220 fprintf(
ioQQQ,
" coll str\n" );
221 fprintf(
ioQQQ,
" hi lo" );
222 for( ilo=0; ilo < nlev-1; ilo++ )
224 fprintf(
ioQQQ,
"%4ld ", ilo );
226 fprintf(
ioQQQ,
" \n" );
228 for( ihi=1; ihi < nlev; ihi++ )
230 fprintf(
ioQQQ,
"%3ld", ihi );
231 for( ilo=0; ilo < ihi; ilo++ )
235 fprintf(
ioQQQ,
"\n" );
238 fprintf(
ioQQQ,
" coll rate\n" );
239 fprintf(
ioQQQ,
" hi lo" );
240 for( ilo=0; ilo < nlev-1; ilo++ )
242 fprintf(
ioQQQ,
"%4ld ", ilo );
244 fprintf(
ioQQQ,
" \n" );
248 for( ihi=1; ihi < nlev; ihi++ )
250 fprintf(
ioQQQ,
"%3ld", ihi );
251 for( ilo=0; ilo < ihi; ilo++ )
255 fprintf(
ioQQQ,
"\n" );
266 for( ilo=0; ilo < (nLevelCalled - 1); ilo++ )
268 for( ihi=ilo + 1; ihi < nLevelCalled; ihi++ )
271 excit[ilo][ihi] =
dsexp((
ex[ihi]-
ex[ilo])*TeInverse);
277 fprintf(
ioQQQ,
" excit, te=%10.2e\n",
phycon.te );
278 fprintf(
ioQQQ,
" hi lo" );
280 for( ilo=0; ilo < (nlev-1); ilo++ )
282 fprintf(
ioQQQ,
"%4ld ", ilo );
284 fprintf(
ioQQQ,
" \n" );
286 for( ihi=1; ihi < nlev; ihi++ )
288 fprintf(
ioQQQ,
"%3ld", ihi );
289 for( ilo=0; ilo < ihi; ilo++ )
291 fprintf(
ioQQQ,
"%10.2e", excit[ilo][ihi] );
293 fprintf(
ioQQQ,
"\n" );
301 for( ilo=0; ilo < (nlev - 1); ilo++ )
303 for( ihi=ilo + 1; ihi < nlev; ihi++ )
307 (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*
g[ilo]/
g[ihi];
315 if( !lgCollRateDone )
317 for( ilo=0; ilo < (nLevelCalled - 1); ilo++ )
319 for( ihi=ilo + 1; ihi < nLevelCalled; ihi++ )
324 (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/
g[ihi]*
dense.cdsqte;
326 (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*
g[ihi]/
g[ilo]*
335 valarray<double> bvec(nlev);
344 for( level=0; level < nlev; level++ )
347 for( ilo=0; ilo < level; ilo++ )
349 double rate = (*CollRate)[level][ilo] + (*AulEscp)[level][ilo] +
350 (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
351 amat[level][level] += rate;
352 amat[level][ilo] = -rate;
356 for( ihi=level + 1; ihi < nlev; ihi++ )
358 double rate = (*CollRate)[level][ihi] + (*AulPump)[level][ihi];
359 amat[level][level] += rate;
360 amat[level][ihi] = -rate;
366 double slowrate = -1.0;
369 for (level=1; level < nlev; ++level)
372 if (rate < slowrate || slowrate < 0.0)
378 if ( slowrate < 3.0 )
380 fprintf(
ioQQQ,
" PROBLEM in atom_levelN: "
381 "non-equilibrium appears important for %s(level=%ld), "
382 "rate * timestep is %g\n",
383 chLabel, slowlevel, slowrate);
388 for( level=0; level < nlev; level++ )
391 bvec[level] =
source[level];
392 if(
amat[level][level] > maxdiag )
393 maxdiag =
amat[level][level];
417 const double CONDITION_SCALE = 1e-10;
418 double conserve_scale = maxdiag*CONDITION_SCALE;
419 ASSERT( conserve_scale > 0.0 );
421 for( level=0; level<nlev; ++level )
423 amat[level][0] += conserve_scale;
425 bvec[0] +=
abund*conserve_scale;
429 fprintf(
ioQQQ,
" amat matrix follows:\n");
430 for( level=0; level < nlev; level++ )
432 for( j=0; j < nlev; j++ )
438 fprintf(
ioQQQ,
" Vector follows:\n");
439 for( j=0; j < nlev; j++ )
441 fprintf(
ioQQQ,
" %.4e" , bvec[j] );
450 fprintf(
ioQQQ,
" atom_levelN: dgetrs finds singular or ill-conditioned matrix\n" );
455 for( level=0; level < nlev; level++ )
458 pops[level] = bvec[level];
465 for( ihi=1; ihi < nlev; ihi++ )
467 for( ilo=0; ilo < ihi; ilo++ )
470 cool1 = (
pops[ilo]*(*CollRate)[ilo][ihi] -
pops[ihi]*(*CollRate)[ihi][ilo])*
474 (*Cool)[ihi][ilo] = cool1;
478 *coolder += dcooldT1;
479 if( dCooldT != NULL )
480 (*dCooldT)[ihi][ilo] = dcooldT1;
489 for( ihi=1; ihi < nlev; ihi++ )
499 for( ihi=0; ihi < nlev; ihi++ )
512 valarray<double> help1(nlev), help2(nlev);
513 double partfun = help1[0] =
g[0];
514 double dpfdT = help2[0] = 0.;
515 for( level=1; level < nlev; level++ )
517 help1[level] =
g[level]*excit[0][level];
518 partfun += help1[level];
519 help2[level] = help1[level]*(
ex[level]-
ex[0])*TeInverse/
phycon.te;
520 dpfdT += help2[level];
525 valarray<double> dndT(nlev);
526 for( level=0; level < nlev; level++ )
528 pops[level] =
abund*help1[level]/partfun;
529 dndT[level] =
abund*(help2[level]*partfun - dpfdT*help1[level])/
pow2(partfun);
536 for( ihi=1; ihi < nlev; ihi++ )
538 for( ilo=0; ilo < ihi; ilo++ )
541 double cool1 = (
pops[ihi]*((*AulEscp)[ihi][ilo] + (*AulPump)[ihi][ilo]) -
542 pops[ilo]*(*AulPump)[ilo][ihi])*deltaE;
545 (*Cool)[ihi][ilo] = cool1;
546 double dcooldT1 = (dndT[ihi]*((*AulEscp)[ihi][ilo] + (*AulPump)[ihi][ilo]) -
547 dndT[ilo]*(*AulPump)[ilo][ihi])*deltaE;
548 *coolder += dcooldT1;
549 if( dCooldT != NULL )
550 (*dCooldT)[ihi][ilo] = dcooldT1;
558 const double poplimit = 1.0e-10;
559 for( level=0; level < nlev; level++ )
561 if(
pops[level] < 0. )
566 *nNegPop = *nNegPop + 1;
581 fprintf(
ioQQQ,
"\n PROBLEM atom_levelN found negative population, nNegPop=%i, atom=%s lgSearch=%c\n",
582 *nNegPop , chLabel ,
TorF(
conv.lgSearch ) );
584 for( level=0; level < nlev; level++ )
589 fprintf(
ioQQQ,
"\n" );
590 for( level=0; level < nlev; level++ )
596 if( lgDeBug || (
trace.lgTrace &&
trace.lgTrLevN) )
598 fprintf(
ioQQQ,
"\n atom_leveln absolute population " );
599 for( level=0; level < nlev; level++ )
601 fprintf(
ioQQQ,
" %10.2e",
pops[level] );
603 fprintf(
ioQQQ,
"\n" );
605 fprintf(
ioQQQ,
" departure coefficients" );
606 for( level=0; level < nlev; level++ )
610 fprintf(
ioQQQ,
"\n\n" );
617 for( ihi=1; ihi < nlev; ihi++ )
619 for( ilo=0; ilo < ihi; ilo++ )
622 (*AulDest)[ilo][ihi] = 0.;
624 (*AulPump)[ihi][ilo] = 0.;
625 (*AulEscp)[ilo][ihi] = 0;
void atom_levelN(long int nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double ***AulEscp, double ***col_str, double ***AulDest, double ***AulPump, double ***CollRate, const double source[], const double sink[], bool lgCollRateDone, double *cooltl, double *coolder, const char *chLabel, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE, multi_arr< double, 2 > *Cool, multi_arr< double, 2 > *dCooldT)