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