44 const double source[] ,
90 for( level=0; level < nlev; level++ )
113 for( ihi=1; ihi < nlev; ihi++ )
115 for( ilo=0; ilo < ihi; ilo++ )
121 ASSERT( (*AulDest)[ilo][ihi] == 0. );
122 ASSERT( (*AulEscp)[ilo][ihi] == 0 );
123 ASSERT( (*AulPump)[ihi][ilo] == 0. );
125 ASSERT( (*AulEscp)[ihi][ilo] >= 0 );
126 ASSERT( (*AulDest)[ihi][ilo] >= 0 );
127 ASSERT( (*col_str)[ihi][ilo] >= 0 );
134 fprintf(
ioQQQ,
" atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
135 fprintf(
ioQQQ,
" AulDest\n" );
137 fprintf(
ioQQQ,
" hi lo" );
139 for( ilo=0; ilo < nlev-1; ilo++ )
141 fprintf(
ioQQQ,
"%4ld ", ilo );
143 fprintf(
ioQQQ,
" \n" );
145 for( ihi=1; ihi < nlev; ihi++ )
147 fprintf(
ioQQQ,
"%3ld", ihi );
148 for( ilo=0; ilo < ihi; ilo++ )
150 fprintf(
ioQQQ,
"%10.2e", (*AulDest)[ihi][ilo] );
152 fprintf(
ioQQQ,
"\n" );
155 fprintf(
ioQQQ,
" A*esc\n" );
156 fprintf(
ioQQQ,
" hi lo" );
157 for( ilo=0; ilo < nlev-1; ilo++ )
159 fprintf(
ioQQQ,
"%4ld ", ilo );
161 fprintf(
ioQQQ,
" \n" );
163 for( ihi=1; ihi < nlev; ihi++ )
165 fprintf(
ioQQQ,
"%3ld", ihi );
166 for( ilo=0; ilo < ihi; ilo++ )
168 fprintf(
ioQQQ,
"%10.2e", (*AulEscp)[ihi][ilo] );
170 fprintf(
ioQQQ,
"\n" );
173 fprintf(
ioQQQ,
" AulPump\n" );
175 fprintf(
ioQQQ,
" hi lo" );
176 for( ilo=0; ilo < nlev-1; ilo++ )
178 fprintf(
ioQQQ,
"%4ld ", ilo );
180 fprintf(
ioQQQ,
" \n" );
182 for( ihi=1; ihi < nlev; ihi++ )
184 fprintf(
ioQQQ,
"%3ld", ihi );
185 for( ilo=0; ilo < ihi; ilo++ )
187 fprintf(
ioQQQ,
"%10.2e", (*AulPump)[ilo][ihi] );
189 fprintf(
ioQQQ,
"\n" );
192 fprintf(
ioQQQ,
" coll str\n" );
193 fprintf(
ioQQQ,
" hi lo" );
194 for( ilo=0; ilo < nlev-1; ilo++ )
196 fprintf(
ioQQQ,
"%4ld ", ilo );
198 fprintf(
ioQQQ,
" \n" );
200 for( ihi=1; ihi < nlev; ihi++ )
202 fprintf(
ioQQQ,
"%3ld", ihi );
203 for( ilo=0; ilo < ihi; ilo++ )
205 fprintf(
ioQQQ,
"%10.2e", (*col_str)[ihi][ilo] );
207 fprintf(
ioQQQ,
"\n" );
210 fprintf(
ioQQQ,
" coll rate\n" );
211 fprintf(
ioQQQ,
" hi lo" );
212 for( ilo=0; ilo < nlev-1; ilo++ )
214 fprintf(
ioQQQ,
"%4ld ", ilo );
216 fprintf(
ioQQQ,
" \n" );
220 for( ihi=1; ihi < nlev; ihi++ )
222 fprintf(
ioQQQ,
"%3ld", ihi );
223 for( ilo=0; ilo < ihi; ilo++ )
225 fprintf(
ioQQQ,
"%10.2e", (*CollRate)[ihi][ilo] );
227 fprintf(
ioQQQ,
"\n" );
241 else if( chExUnits==
'w' )
251 for( ilo=0; ilo < (nlev - 1); ilo++ )
253 for( ihi=ilo + 1; ihi < nlev; ihi++ )
256 excit[ilo][ihi] =
sexp((ex[ihi]-ex[ilo])*TeInverse);
263 fprintf(
ioQQQ,
" hi lo" );
265 for( ilo=0; ilo < (nlev-1); ilo++ )
267 fprintf(
ioQQQ,
"%4ld ", ilo );
269 fprintf(
ioQQQ,
" \n" );
271 for( ihi=1; ihi < nlev; ihi++ )
273 fprintf(
ioQQQ,
"%3ld", ihi );
274 for( ilo=0; ilo < ihi; ilo++ )
276 fprintf(
ioQQQ,
"%10.2e", excit[ilo][ihi] );
278 fprintf(
ioQQQ,
"\n" );
284 if( (excit[0][nlev-1] + (*AulPump)[0][nlev-1] < 1e-13 ) && (source[nlev-1]==0.) )
293 for( level=1; level < nlev; level++ )
313 for( ilo=0; ilo < (nlev - 1); ilo++ )
315 for( ihi=ilo + 1; ihi < nlev; ihi++ )
319 (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*g[ilo]/g[ihi];
325 if( !lgCollRateDone )
327 for( ilo=0; ilo < (nlev - 1); ilo++ )
329 for( ihi=ilo + 1; ihi < nlev; ihi++ )
332 ASSERT( (*col_str)[ihi][ilo]>= 0. );
334 (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/g[ihi]*
dense.
cdsqte;
336 (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
350 lgHomogeneous =
false;
351 for( level=0; level < nlev; level++ )
353 bvec[level] = source[level];
357 lgHomogeneous =
true;
362 for( level=0; level < nlev; level++ )
364 amat[level][level] = sink[level];
367 for( ilo=0; ilo < level; ilo++ )
369 amat[level][level] += (*CollRate)[level][ilo] + (*AulEscp)[level][ilo] +
370 (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
373 amat[ilo][level] = -(*CollRate)[ilo][level] - (*AulPump)[ilo][level];
377 for( ihi=level + 1; ihi < nlev; ihi++ )
379 amat[level][level] += (*CollRate)[level][ihi] + (*AulPump)[level][ihi];
381 amat[ihi][level] = -(*CollRate)[ihi][level] - (*AulEscp)[ihi][level] -
382 (*AulDest)[ihi][level] - (*AulPump)[ihi][level];
392 for( level=0; level<nlev; ++level )
394 amat[level][0] = 1.0;
402 fprintf(
ioQQQ,
" amat matrix follows:\n");
403 for( level=0; level < nlev; level++ )
405 for( j=0; j < nlev; j++ )
407 fprintf(
ioQQQ,
" %.4e" , amat[level][j]);
413 fprintf(
ioQQQ,
" Sum creation zero so pop sum used\n");
417 fprintf(
ioQQQ,
" Sum creation non-zero (%e), vector follows:\n",sum);
418 for( j=0; j < nlev; j++ )
420 fprintf(
ioQQQ,
" %.4e" , bvec[j] );
437 fprintf(
ioQQQ,
" atom_levelN: dgetrs finds singular or ill-conditioned matrix\n" );
442 for( level=0; level < nlev; level++ )
445 pops[level] = bvec[level];
452 for( ihi=1; ihi < nlev; ihi++ )
454 for( ilo=0; ilo < ihi; ilo++ )
457 cool1 = (pops[ilo]*(*CollRate)[ilo][ihi] - pops[ihi]*(*CollRate)[ihi][ilo])*
472 if( pops[0] > 1e-20 && excit[0][nlev-1] > 1e-20 )
476 for( ihi=1; ihi < nlev; ihi++ )
479 depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit[0][ihi];
486 for( ihi=0; ihi < nlev; ihi++ )
496 for( level=0; level < nlev; level++ )
498 if( pops[level] < 0. )
506 fprintf(
ioQQQ,
"\n PROBLEM atom_levelN found negative population, atom=%s\n",
508 fprintf(
ioQQQ,
" absolute population=" );
510 for( level=0; level < nlev; level++ )
512 fprintf(
ioQQQ,
"%10.2e", pops[level] );
515 fprintf(
ioQQQ,
"\n" );
516 for( level=0; level < nlev; level++ )
518 pops[level] = (double)
MAX2(0.,pops[level]);
524 fprintf(
ioQQQ,
"\n atom_leveln absolute population " );
525 for( level=0; level < nlev; level++ )
527 fprintf(
ioQQQ,
" %10.2e", pops[level] );
529 fprintf(
ioQQQ,
"\n" );
531 fprintf(
ioQQQ,
" departure coefficients" );
532 for( level=0; level < nlev; level++ )
534 fprintf(
ioQQQ,
" %10.2e", depart[level] );
536 fprintf(
ioQQQ,
"\n\n" );
543 for( ihi=1; ihi < nlev; ihi++ )
545 for( ilo=0; ilo < ihi; ilo++ )
548 (*AulDest)[ilo][ihi] = 0.;
550 (*AulPump)[ihi][ilo] = 0.;
551 (*AulEscp)[ilo][ihi] = 0;