00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "thermal.h"
00007 #include "conv.h"
00008 #include "phycon.h"
00009 #include "dense.h"
00010 #include "trace.h"
00011 #include "thirdparty.h"
00012 #include "atoms.h"
00013
00014 void atom_levelN(
00015
00016 long int nLevelCalled,
00017
00018
00019 realnum abund,
00020
00021 const double g[],
00022
00023
00024 const double ex[],
00025
00026 char chExUnits,
00027
00028 double pops[],
00029
00030 double depart[],
00031
00032 double ***AulEscp,
00033
00034 double ***col_str,
00035
00036
00037 double ***AulDest,
00038
00039 double ***AulPump,
00040
00041
00042
00043 double ***CollRate,
00044
00045 const double source[] ,
00046
00047 const double sink[] ,
00048
00049
00050 bool lgCollRateDone,
00051
00052 double *cooltl,
00053 double *coolder,
00054
00055 const char *chLabel,
00056
00057
00058
00059
00060
00061 int *nNegPop,
00062
00063 bool *lgZeroPop ,
00064
00065 bool lgDeBug )
00066 {
00067 bool lgHomogeneous;
00068
00069 long int level,
00070 ihi,
00071 ilo,
00072 j;
00073
00074 int32 ner;
00075
00076 double cool1,
00077 TeInverse,
00078 TeConvFac,
00079 sum;
00080
00081 DEBUG_ENTRY( "atom_levelN()" );
00082
00083 *nNegPop = -1;
00084
00085
00086 if( chExUnits=='K' )
00087 {
00088
00089
00090 TeInverse = 1./phycon.te;
00091
00092 TeConvFac = 1.;
00093 }
00094 else if( chExUnits=='w' )
00095 {
00096
00097 TeInverse = 1./phycon.te_wn;
00098 TeConvFac = T1CM;
00099 }
00100 else
00101 TotalInsanity();
00102
00103 long int nlev = nLevelCalled;
00104
00105 while( ( (dsexp((ex[nlev-1]-ex[0])*TeInverse) + (*AulPump)[0][nlev-1]) < SMALLFLOAT ) &&
00106 source[nlev-1]==0. && nlev > 1)
00107 {
00108 pops[nlev-1] = 0.;
00109 depart[nlev-1] = 0.;
00110 --nlev;
00111 }
00112
00113
00114 ASSERT( abund>= 0. );
00115 if( abund == 0. || nlev==1 )
00116 {
00117 *cooltl = 0.;
00118 *coolder = 0.;
00119
00120 *nNegPop = 0;
00121 *lgZeroPop = true;
00122
00123 pops[0] = abund;
00124 depart[0] = 1.;
00125 for( level=1; level < nlev; level++ )
00126 {
00127 pops[level] = 0.;
00128 depart[level] = 0.;
00129 }
00130
00131
00132
00133 return;
00134 }
00135
00136
00137 auto_vec<int32> ipiv( new int32[nlev] );
00138 auto_vec<double> bvec( new double[nlev] );
00139 multi_arr<double,2,C_TYPE> amat(nlev,nlev);
00140 multi_arr<double,2> excit(nLevelCalled,nLevelCalled);
00141
00142 # ifndef NDEBUG
00143
00144 ASSERT( ex[0] == 0. );
00145
00146 for( ihi=1; ihi < nlev; ihi++ )
00147 {
00148 for( ilo=0; ilo < ihi; ilo++ )
00149 {
00150
00151
00152
00153
00154 ASSERT( (*AulDest)[ilo][ihi] == 0. );
00155 ASSERT( (*AulEscp)[ilo][ihi] == 0 );
00156 ASSERT( (*AulPump)[ihi][ilo] == 0. );
00157
00158 ASSERT( (*AulEscp)[ihi][ilo] >= 0 );
00159 ASSERT( (*AulDest)[ihi][ilo] >= 0 );
00160 ASSERT( (*col_str)[ihi][ilo] >= 0 );
00161 }
00162 }
00163 # endif
00164
00165 if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
00166 {
00167 fprintf( ioQQQ, " atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
00168 fprintf( ioQQQ, " AulDest\n" );
00169
00170 fprintf( ioQQQ, " hi lo" );
00171
00172 for( ilo=0; ilo < nlev-1; ilo++ )
00173 {
00174 fprintf( ioQQQ, "%4ld ", ilo );
00175 }
00176 fprintf( ioQQQ, " \n" );
00177
00178 for( ihi=1; ihi < nlev; ihi++ )
00179 {
00180 fprintf( ioQQQ, "%3ld", ihi );
00181 for( ilo=0; ilo < ihi; ilo++ )
00182 {
00183 fprintf( ioQQQ, "%10.2e", (*AulDest)[ihi][ilo] );
00184 }
00185 fprintf( ioQQQ, "\n" );
00186 }
00187
00188 fprintf( ioQQQ, " A*esc\n" );
00189 fprintf( ioQQQ, " hi lo" );
00190 for( ilo=0; ilo < nlev-1; ilo++ )
00191 {
00192 fprintf( ioQQQ, "%4ld ", ilo );
00193 }
00194 fprintf( ioQQQ, " \n" );
00195
00196 for( ihi=1; ihi < nlev; ihi++ )
00197 {
00198 fprintf( ioQQQ, "%3ld", ihi );
00199 for( ilo=0; ilo < ihi; ilo++ )
00200 {
00201 fprintf( ioQQQ, "%10.2e", (*AulEscp)[ihi][ilo] );
00202 }
00203 fprintf( ioQQQ, "\n" );
00204 }
00205
00206 fprintf( ioQQQ, " AulPump\n" );
00207
00208 fprintf( ioQQQ, " hi lo" );
00209 for( ilo=0; ilo < nlev-1; ilo++ )
00210 {
00211 fprintf( ioQQQ, "%4ld ", ilo );
00212 }
00213 fprintf( ioQQQ, " \n" );
00214
00215 for( ihi=1; ihi < nlev; ihi++ )
00216 {
00217 fprintf( ioQQQ, "%3ld", ihi );
00218 for( ilo=0; ilo < ihi; ilo++ )
00219 {
00220 fprintf( ioQQQ, "%10.2e", (*AulPump)[ilo][ihi] );
00221 }
00222 fprintf( ioQQQ, "\n" );
00223 }
00224
00225 fprintf( ioQQQ, " coll str\n" );
00226 fprintf( ioQQQ, " hi lo" );
00227 for( ilo=0; ilo < nlev-1; ilo++ )
00228 {
00229 fprintf( ioQQQ, "%4ld ", ilo );
00230 }
00231 fprintf( ioQQQ, " \n" );
00232
00233 for( ihi=1; ihi < nlev; ihi++ )
00234 {
00235 fprintf( ioQQQ, "%3ld", ihi );
00236 for( ilo=0; ilo < ihi; ilo++ )
00237 {
00238 fprintf( ioQQQ, "%10.2e", (*col_str)[ihi][ilo] );
00239 }
00240 fprintf( ioQQQ, "\n" );
00241 }
00242
00243 fprintf( ioQQQ, " coll rate\n" );
00244 fprintf( ioQQQ, " hi lo" );
00245 for( ilo=0; ilo < nlev-1; ilo++ )
00246 {
00247 fprintf( ioQQQ, "%4ld ", ilo );
00248 }
00249 fprintf( ioQQQ, " \n" );
00250
00251 if( lgCollRateDone )
00252 {
00253 for( ihi=1; ihi < nlev; ihi++ )
00254 {
00255 fprintf( ioQQQ, "%3ld", ihi );
00256 for( ilo=0; ilo < ihi; ilo++ )
00257 {
00258 fprintf( ioQQQ, "%10.2e", (*CollRate)[ihi][ilo] );
00259 }
00260 fprintf( ioQQQ, "\n" );
00261 }
00262 }
00263 }
00264
00265
00266
00267
00268 for( ilo=0; ilo < (nLevelCalled - 1); ilo++ )
00269 {
00270 for( ihi=ilo + 1; ihi < nLevelCalled; ihi++ )
00271 {
00272
00273 excit[ilo][ihi] = dsexp((ex[ihi]-ex[ilo])*TeInverse);
00274 }
00275 }
00276
00277 if( trace.lgTrace && trace.lgTrLevN )
00278 {
00279 fprintf( ioQQQ, " excit, te=%10.2e\n", phycon.te );
00280 fprintf( ioQQQ, " hi lo" );
00281
00282 for( ilo=0; ilo < (nlev-1); ilo++ )
00283 {
00284 fprintf( ioQQQ, "%4ld ", ilo );
00285 }
00286 fprintf( ioQQQ, " \n" );
00287
00288 for( ihi=1; ihi < nlev; ihi++ )
00289 {
00290 fprintf( ioQQQ, "%3ld", ihi );
00291 for( ilo=0; ilo < ihi; ilo++ )
00292 {
00293 fprintf( ioQQQ, "%10.2e", excit[ilo][ihi] );
00294 }
00295 fprintf( ioQQQ, "\n" );
00296 }
00297 }
00298
00299
00300 *lgZeroPop = false;
00301
00302
00303 for( ilo=0; ilo < (nlev - 1); ilo++ )
00304 {
00305 for( ihi=ilo + 1; ihi < nlev; ihi++ )
00306 {
00307
00308
00309 (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*g[ilo]/g[ihi];
00310 }
00311 }
00312
00313
00314
00315
00316
00317 if( !lgCollRateDone )
00318 {
00319 for( ilo=0; ilo < (nLevelCalled - 1); ilo++ )
00320 {
00321 for( ihi=ilo + 1; ihi < nLevelCalled; ihi++ )
00322 {
00323
00324 ASSERT( (*col_str)[ihi][ilo]>= 0. );
00325
00326 (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/g[ihi]*dense.cdsqte;
00327
00328 (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
00329 excit[ilo][ihi];
00330 }
00331 }
00332 }
00333
00334
00335 amat.zero();
00336
00337
00338
00339
00340
00341 sum = 0.;
00342 lgHomogeneous = false;
00343 for( level=0; level < nlev; level++ )
00344 {
00345 bvec[level] = source[level];
00346 sum += bvec[level];
00347 }
00348 if( sum==0. )
00349 lgHomogeneous = true;
00350
00351
00352
00353
00354 for( level=0; level < nlev; level++ )
00355 {
00356 amat[level][level] = sink[level];
00357
00358
00359 for( ilo=0; ilo < level; ilo++ )
00360 {
00361 amat[level][level] += (*CollRate)[level][ilo] + (*AulEscp)[level][ilo] +
00362 (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
00363
00364
00365 amat[ilo][level] = -(*CollRate)[ilo][level] - (*AulPump)[ilo][level];
00366 }
00367
00368
00369 for( ihi=level + 1; ihi < nlev; ihi++ )
00370 {
00371 amat[level][level] += (*CollRate)[level][ihi] + (*AulPump)[level][ihi];
00372
00373 amat[ihi][level] = -(*CollRate)[ihi][level] - (*AulEscp)[ihi][level] -
00374 (*AulDest)[ihi][level] - (*AulPump)[ihi][level];
00375
00376 }
00377 }
00378
00379
00380
00381
00382 if( lgHomogeneous )
00383 {
00384 for( level=0; level<nlev; ++level )
00385 {
00386 amat[level][0] = 1.0;
00387 }
00388
00389 bvec[0] = abund;
00390 }
00391
00392 if( lgDeBug )
00393 {
00394 fprintf(ioQQQ," amat matrix follows:\n");
00395 for( level=0; level < nlev; level++ )
00396 {
00397 for( j=0; j < nlev; j++ )
00398 {
00399 fprintf(ioQQQ," %.4e" , amat[level][j]);
00400 }
00401 fprintf(ioQQQ,"\n");
00402 }
00403 if( sum==0. )
00404 {
00405 fprintf(ioQQQ," Sum creation zero so pop sum used\n");
00406 }
00407 else
00408 {
00409 fprintf(ioQQQ," Sum creation non-zero (%e), vector follows:\n",sum);
00410 for( j=0; j < nlev; j++ )
00411 {
00412 fprintf(ioQQQ," %.4e" , bvec[j] );
00413 }
00414 fprintf(ioQQQ,"\n");
00415 }
00416 }
00417
00418 ner = 0;
00419 getrf_wrapper(nlev,nlev,amat.data(),nlev,ipiv.data(),&ner);
00420
00421
00422
00423
00424
00425 getrs_wrapper('N',nlev,1,amat.data(),nlev,ipiv.data(),bvec.data(),nlev,&ner);
00426
00427 if( ner != 0 )
00428 {
00429 fprintf( ioQQQ, " atom_levelN: dgetrs finds singular or ill-conditioned matrix\n" );
00430 cdEXIT(EXIT_FAILURE);
00431 }
00432
00433
00434 for( level=0; level < nlev; level++ )
00435 {
00436
00437 pops[level] = bvec[level];
00438 }
00439
00440
00441
00442 *cooltl = 0.;
00443 *coolder = 0.;
00444 for( ihi=1; ihi < nlev; ihi++ )
00445 {
00446 for( ilo=0; ilo < ihi; ilo++ )
00447 {
00448
00449 cool1 = (pops[ilo]*(*CollRate)[ilo][ihi] - pops[ihi]*(*CollRate)[ihi][ilo])*
00450 (ex[ihi] - ex[ilo]);
00451 *cooltl += cool1;
00452
00453
00454
00455 *coolder += cool1*( (ex[ihi] - ex[0])*thermal.tsq1 - thermal.halfte);
00456 }
00457 }
00458
00459
00460 *cooltl *= BOLTZMANN*TeConvFac;
00461 *coolder *= BOLTZMANN*TeConvFac;
00462
00463
00464 if( pops[0] > SMALLFLOAT && excit[0][nlev-1] > SMALLFLOAT )
00465 {
00466
00467 depart[0] = 1.;
00468 for( ihi=1; ihi < nlev; ihi++ )
00469 {
00470
00471 depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit[0][ihi];
00472 }
00473 }
00474
00475 else
00476 {
00477
00478 for( ihi=0; ihi < nlev; ihi++ )
00479 {
00480
00481 depart[ihi] = 0.;
00482 }
00483 depart[0] = 1.;
00484 }
00485
00486
00487 *nNegPop = 0;
00488
00489 double poplimit = 1.0E-10;
00490 for( level=0; level < nlev; level++ )
00491 {
00492 if( pops[level] < 0. )
00493 {
00494 if( fabs(pops[level]/abund) > poplimit )
00495 {
00496
00497 *nNegPop = 1;
00498 }
00499 else
00500 {
00501 pops[level] = SMALLFLOAT;
00502
00503
00504 }
00505 }
00506 }
00507
00508 if( *nNegPop!=0 )
00509 {
00510 ASSERT( *nNegPop==1 );
00511
00512 fprintf( ioQQQ, "\n PROBLEM atom_levelN found negative population, nNegPop=%i3, atom=%s lgSearch=%c\n",
00513 *nNegPop , chLabel , TorF( conv.lgSearch ) );
00514
00515 for( level=0; level < nlev; level++ )
00516 {
00517 fprintf( ioQQQ, "%10.2e", pops[level] );
00518 }
00519
00520 fprintf( ioQQQ, "\n" );
00521 for( level=0; level < nlev; level++ )
00522 {
00523 pops[level] = (double)MAX2(0.,pops[level]);
00524 }
00525 }
00526
00527 if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
00528 {
00529 fprintf( ioQQQ, "\n atom_leveln absolute population " );
00530 for( level=0; level < nlev; level++ )
00531 {
00532 fprintf( ioQQQ, " %10.2e", pops[level] );
00533 }
00534 fprintf( ioQQQ, "\n" );
00535
00536 fprintf( ioQQQ, " departure coefficients" );
00537 for( level=0; level < nlev; level++ )
00538 {
00539 fprintf( ioQQQ, " %10.2e", depart[level] );
00540 }
00541 fprintf( ioQQQ, "\n\n" );
00542 }
00543
00544 # ifndef NDEBUG
00545
00546
00547
00548 for( ihi=1; ihi < nlev; ihi++ )
00549 {
00550 for( ilo=0; ilo < ihi; ilo++ )
00551 {
00552
00553 (*AulDest)[ilo][ihi] = 0.;
00554
00555 (*AulPump)[ihi][ilo] = 0.;
00556 (*AulEscp)[ilo][ihi] = 0;
00557 }
00558 }
00559 # endif
00560
00561
00562
00563 ASSERT( *nNegPop>=0 );
00564
00565 return;
00566 }