20 void dgeev_(
const char*,
const char*,
int*,
double*,
int*,
double*,
22 double*,
int*,
double*,
int*,
double*,
int*,
int*,
int,
int);
37 for(
long ihi=1; ihi < nlev; ++ihi )
39 double fac =
dsexp((ex[ihi]-ex[ihi-1])*TeInverse);
40 for(
long ilo=0; ilo < ihi-1; ++ilo )
44 excit[ihi][ihi-1] = fac;
51 for(
long ilo=0; ilo < nlev-1; ilo++ )
57 for(
long ihi=1; ihi < nlev; ihi++ )
60 for(
long ilo=0; ilo < ihi; ilo++ )
70 for(
long ilo=0; ilo < (nlev-1); ilo++ )
76 for(
long ihi=1; ihi < nlev; ihi++ )
79 for(
long ilo=0; ilo < ihi; ilo++ )
87 for(
long ihi=1; ihi < nlev; ihi++ )
89 for(
long ilo=0; ilo < ihi; ilo++ )
92 ASSERT( (col_str)[ihi][ilo]>= 0. );
94 (CollRate)[ihi][ilo] = (col_str)[ihi][ilo]/g[ihi]*
dense.
cdsqte;
96 (CollRate)[ilo][ihi] = (CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
103 long int nLevelCalled,
130 const double source[],
139 const bool lgPrtMatrix,
173 if (grnd_excit != NULL) *grnd_excit = 0.0;
184 else if( chExUnits==
'w' )
194 for( ihi=1; ihi < nLevelCalled; ++ihi )
195 arg[ihi] = -(ex[ihi]-ex[0])*TeInverse;
196 vexp( arg.ptr0(), excit_gnd.
ptr0(), 1, nLevelCalled );
198 long int nlev = nLevelCalled;
200 while( nlev > 1 && excit_gnd[nlev-1]+AulPump[0][nlev-1] <
SMALLFLOAT &&
201 source[nlev-1] == 0. )
210 if( abund == 0. || nlev==1 )
220 for( level=1; level < nlev; level++ )
232 for( ihi=0; ihi < nlev; ihi++ )
234 for( ilo=0; ilo < nlev; ilo++ )
238 ASSERT( ihi == ilo || (AulDest)[ihi][ilo] >= 0. );
239 ASSERT( ihi == ilo || (AulEscp)[ihi][ilo] >= 0 );
243 for( ihi=1; ihi < nlev; ihi++ )
245 for( ilo=0; ilo < ihi; ilo++ )
247 ASSERT( (AulPump)[ihi][ilo] >= 0. );
254 fprintf(
ioQQQ,
" atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
259 for( ilo=0; ilo < nlev-1; ilo++ )
265 for( ihi=1; ihi < nlev; ihi++ )
268 for( ilo=0; ilo < ihi; ilo++ )
277 for( ilo=0; ilo < nlev-1; ilo++ )
283 for( ihi=1; ihi < nlev; ihi++ )
286 for( ilo=0; ilo < ihi; ilo++ )
296 for( ilo=0; ilo < nlev-1; ilo++ )
302 for( ihi=1; ihi < nlev; ihi++ )
305 for( ilo=0; ilo < ihi; ilo++ )
314 for( ilo=0; ilo < nlev-1; ilo++ )
320 for( ihi=1; ihi < nlev; ihi++ )
323 for( ilo=0; ilo < ihi; ilo++ )
345 for( level=0; level < nlev; level++ )
348 for( ilo=0; ilo < level; ilo++ )
350 double rate = (CollRate)[level][ilo] + (AulEscp)[level][ilo] +
351 (AulDest)[level][ilo] + (AulPump)[level][ilo];
352 amat[level][level] += rate;
353 amat[level][ilo] = -rate;
357 for( ihi=level + 1; ihi < nlev; ihi++ )
359 double rate = (CollRate)[level][ihi] + (AulPump)[level][ihi];
360 amat[level][level] += rate;
361 amat[level][ihi] = -rate;
364 if (grnd_excit != NULL) *grnd_excit =
amat[0][0];
369 double slowrate = -1.0;
372 for (level=1; level < nlev; ++level)
375 if (rate < slowrate || slowrate < 0.0)
381 if ( slowrate < 3.0 )
383 static map<string,double> failures;
385 oss << chLabel <<
"[level=" << slowlevel <<
']';
386 string str=oss.str();
387 if (failures.find(str) == failures.end())
389 failures[str] = slowrate;
393 "non-equilibrium appears important for %s, "
394 "rate * timestep is %g\n",
395 str.c_str(), slowrate);
399 if ( slowrate < failures[str])
400 failures[str] = slowrate;
406 for( level=0; level < nlev; level++ )
409 bvec[level] = source[level];
410 if(
amat[level][level] > maxdiag )
411 maxdiag =
amat[level][level];
412 amat[level][level] += sink[level];
435 const double CONDITION_SCALE = 1e-10;
436 double conserve_scale = maxdiag*CONDITION_SCALE;
437 ASSERT( conserve_scale > 0.0 );
439 for( level=0; level<nlev; ++level )
441 amat[level][0] += conserve_scale;
443 bvec[0] += abund*conserve_scale;
448 for( level=0; level < nlev; level++ )
450 for( j=0; j < nlev; j++ )
457 for( j=0; j < nlev; j++ )
468 vl(nlev,nlev),vr(nlev,nlev);
472 for( level=0; level < nlev; level++ )
475 for( ilo=0; ilo < level; ilo++ )
477 double rate = (*CollRate)[level][ilo];
478 mcol[level][level] += rate;
479 mcol[level][ilo] = -rate;
480 rate = (*AulEscp)[level][ilo] +
481 (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
482 mrad[level][level] += rate;
483 mrad[level][ilo] = -rate;
486 for( ihi=level + 1; ihi < nlev; ihi++ )
488 double rate = (*CollRate)[level][ihi];
489 mcol[level][level] += rate;
490 mcol[level][ihi] = -rate;
491 rate = (*AulPump)[level][ihi];
492 mrad[level][level] += rate;
493 mrad[level][ihi] = -rate;
497 for( ilo=0; ilo < nlev; ilo++ )
499 for( ihi=0; ihi < nlev; ++ihi)
501 mtst[ilo][ihi] = mcol[ilo][ihi];
504 const char job[2] =
"V";
505 int lwork = 4*nlev,info,inlev=nlev;
506 vector<double> wr(nlev),wi(nlev),work(lwork);
511 for( level=0; level < nlev; level++ )
513 fprintf(
ioQQQ,
"%ld %15.8g+%15.8gi\n",level,wr[level],wi[level]);
522 for( level=0; level < nlev; level++ )
525 for( ilo=0; ilo < nlev; ilo++ )
527 tot += vl[level][ilo]*vr[level][ilo];
530 for( ilo=0; ilo < nlev; ilo++ )
532 vl[level][ilo] *= tot;
537 for( ilo=0; ilo < nlev; ilo++ )
539 for( ihi=0; ihi < nlev; ihi++ )
541 for( level=0; level < nlev; level++ )
543 mtst[ilo][ihi] += mrad[level][ihi] * vr[ilo][level];
549 for( ilo=0; ilo < nlev; ilo++ )
551 for( ihi=0; ihi < nlev; ihi++ )
553 for( level=0; level < nlev; level++ )
555 mtst1[ilo][ihi] += vl[ihi][level]*mtst[ilo][level];
559 for( level=0; level < nlev; level++ )
561 fprintf(
ioQQQ,
"%ld %15.8g %15.8g\n",level,wr[level],mtst1[level][level]);
575 fprintf(
ioQQQ,
" PROBLEM atom_levelN: dgetrs finds singular or ill-conditioned matrix for %s\n",
581 for( level=0; level < nlev; level++ )
584 pops[level] =
bvec[level];
591 for( ihi=1; ihi < nlev; ihi++ )
593 for( ilo=0; ilo < ihi; ilo++ )
596 cool1 = (pops[ilo]*(CollRate)[ilo][ihi] - pops[ihi]*(CollRate)[ihi][ilo])*
597 (ex[ihi] - ex[ilo])*BOLTZMANN*TeConvFac;
600 (*Cool)[ihi][ilo] = cool1;
604 *coolder += dcooldT1;
605 if( dCooldT != NULL )
606 (*dCooldT)[ihi][ilo] = dcooldT1;
615 for( ihi=1; ihi < nlev; ihi++ )
618 depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit_gnd[ihi];
625 for( ihi=0; ihi < nlev; ihi++ )
638 valarray<double> help1(nlev), help2(nlev);
639 double partfun = help1[0] = g[0];
640 double dpfdT = help2[0] = 0.;
641 for( level=1; level < nlev; level++ )
643 help1[level] = g[level]*excit_gnd[level];
644 partfun += help1[level];
645 help2[level] = help1[level]*(ex[level]-ex[0])*TeInverse/
phycon.
te;
646 dpfdT += help2[level];
651 valarray<double> dndT(nlev);
652 for( level=0; level < nlev; level++ )
654 pops[level] = abund*help1[level]/partfun;
655 dndT[level] = abund*(help2[level]*partfun - dpfdT*help1[level])/
pow2(partfun);
662 for( ihi=1; ihi < nlev; ihi++ )
664 for( ilo=0; ilo < ihi; ilo++ )
666 double deltaE = (ex[ihi] - ex[ilo])*BOLTZMANN*TeConvFac;
667 double cool1 = (pops[ihi]*((AulEscp)[ihi][ilo] + (AulPump)[ihi][ilo]) -
668 pops[ilo]*(AulPump)[ilo][ihi])*deltaE;
671 (*Cool)[ihi][ilo] = cool1;
672 double dcooldT1 = (dndT[ihi]*((AulEscp)[ihi][ilo] + (AulPump)[ihi][ilo]) -
673 dndT[ilo]*(AulPump)[ilo][ihi])*deltaE;
674 *coolder += dcooldT1;
675 if( dCooldT != NULL )
676 (*dCooldT)[ihi][ilo] = dcooldT1;
684 const double poplimit = 1.0e-10;
685 for( level=0; level < nlev; level++ )
687 if( pops[level] < 0. )
689 if( fabs(pops[level]/abund) > poplimit )
692 *nNegPop = *nNegPop + 1;
707 fprintf(
ioQQQ,
"\n PROBLEM atom_levelN found negative population, nNegPop=%i, atom=%s lgSearch=%c\n",
710 for( level=0; level < nlev; level++ )
716 for( level=0; level < nlev; level++ )
718 pops[level] = (double)
MAX2(0.,pops[level]);
724 fprintf(
ioQQQ,
"\n atom_leveln absolute population " );
725 for( level=0; level < nlev; level++ )
732 for( level=0; level < nlev; level++ )
743 for( ilo=0; ilo < nlev; ilo++ )
745 for( ihi=ilo+1; ihi < nlev; ihi++ )
748 AulDest[ilo][ihi] = 0.;
751 for( ihi=1; ihi < nlev; ihi++ )
753 for( ilo=0; ilo < ihi; ilo++ )
756 AulPump[ihi][ilo] = 0.;
759 for( ilo=0; ilo < nlev; ilo++ )
761 for( ihi=ilo+1; ihi < nlev; ihi++ )
763 AulEscp[ilo][ihi] = 0;
NORETURN void TotalInsanity(void)
double depart(const genericState &gs)
void resize(long int nlev)
void prtRates(const long nlevels_local, const multi_arr< double, 2, C_TYPE > &a, valarray< double > &b)
bool lgTimeDependentStatic
void vexp(const double x[], double y[], long nlo, long nhi)
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
void operator()(long nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double **AulEscp, double **AulDest, double **AulPump, double **CollRate, const double create[], const double destroy[], double *cooltl, double *coolder, const char *chLabel, bool lgPrtMatrix, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE=false, multi_arr< double, 2 > *Cool=NULL, multi_arr< double, 2 > *dCooldT=NULL, double *grnd_excit=NULL)
multi_arr< double, 2, C_TYPE > amat
#define DEBUG_ENTRY(funcname)
int fprintf(const Output &stream, const char *format,...)
void operator()(long nlev, double TeInverse, double **col_str, const double ex[], const double g[], double **CollRate)