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