cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*atom_levelN compute an arbitrary N level atom */ 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 /* nlev is the number of levels to compute*/ 00015 long int nlev, 00016 /* ABUND is total abundance of species, used for nth equation 00017 * if balance equations are homogeneous */ 00018 realnum abund, 00019 /* G(nlev) is stat weight of levels */ 00020 const double g[], 00021 /* EX(nlev) is excitation potential of levels, deg K or wavenumbers 00022 * 0 for lowest level, all are energy rel to ground NOT d(ENER) */ 00023 const double ex[], 00024 /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */ 00025 char chExUnits, 00026 /* populations [cm-3] of each level as deduced here */ 00027 double pops[], 00028 /* departure coefficient, derived below */ 00029 double depart[], 00030 /* net transition rate, A * esc prob, s-1 */ 00031 double ***AulEscp, 00032 /* col str from up to low */ 00033 double ***col_str, 00034 /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob, 00035 * asserts confirm that [ihi][ilo] is zero */ 00036 double ***AulDest, 00037 /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */ 00038 double ***AulPump, 00039 /* collision rates (s^-1), evaluated here and returned for cooling by calling function, 00040 * unless following flag is true. If true then calling function has already filled 00041 * in these rates. CollRate[i][j] is rate from i to j */ 00042 double ***CollRate, 00043 /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */ 00044 const double source[] , 00045 /* this is an additional destruction rate to continuum, normally zero, units s-1 */ 00046 const double sink[] , 00047 /* flag saying whether CollRate already done, or we need to do it here, 00048 * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/ 00049 bool lgCollRateDone, 00050 /* total cooling and its derivative, set here but nothing done with it*/ 00051 double *cooltl, 00052 double *coolder, 00053 /* string used to identify calling program in case of error */ 00054 const char *chLabel, 00055 /* nNegPop flag indicating what we have done 00056 * positive if negative populations occurred 00057 * zero if normal calculation done 00058 * negative if too cold (for some atoms other routine will be called in this case) */ 00059 int *nNegPop, 00060 /* true if populations are zero, either due to zero abundance of very low temperature */ 00061 bool *lgZeroPop , 00062 /* option to print debug information */ 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 /* check for zero abundance and exit if so */ 00082 if( abund <= 0. ) 00083 { 00084 *cooltl = 0.; 00085 *coolder = 0.; 00086 /* says calc was ok */ 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 /* there are TWO abort returns in this sub, 00099 * this one is for zero abundance */ 00100 return; 00101 } 00102 00103 // these are all automatically deallocated when they go out of scope 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 /* excitation temperature of lowest level must be zero */ 00111 ASSERT( ex[0] == 0. ); 00112 00113 for( ihi=1; ihi < nlev; ihi++ ) 00114 { 00115 for( ilo=0; ilo < ihi; ilo++ ) 00116 { 00117 /* following must be zero: 00118 * AulDest[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low 00119 * AulEscp[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low 00120 * AulPump[ihi][ilo] - so that pumping only proceeds from low energy to high */ 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 /* >>chng 05 dec 14, units of ex[] can be Kelvin (old default) or wavenumbers */ 00233 if( chExUnits=='K' ) 00234 { 00235 /* ex[] is in temperature units - this will multiply ex[] to 00236 * obtain Boltzmann factor */ 00237 TeInverse = 1./phycon.te; 00238 /* this multiplies ex[] to obtain energy difference between levels */ 00239 TeConvFac = 1.; 00240 } 00241 else if( chExUnits=='w' ) 00242 { 00243 /* ex[] is in wavenumber units */ 00244 TeInverse = 1./phycon.te_wn; 00245 TeConvFac = T1CM; 00246 } 00247 else 00248 TotalInsanity(); 00249 00250 /* find set of Boltzmann factors */ 00251 for( ilo=0; ilo < (nlev - 1); ilo++ ) 00252 { 00253 for( ihi=ilo + 1; ihi < nlev; ihi++ ) 00254 { 00255 /* >>chng 05 dec 14, option to have ex be either Kelvin or wavenumbers */ 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 /* punt if total excitation rate is zero */ 00283 /* >>chng 04 sep 16, add test on non-zero source */ 00284 if( (excit[0][nlev-1] + (*AulPump)[0][nlev-1] < 1e-13 ) && (source[nlev-1]==0.) ) 00285 { 00286 *cooltl = 0.; 00287 *coolder = 0.; 00288 /* special flag saying too cool for highest level to be computed. 00289 * some routines will call another routine for lower levels in this case */ 00290 *nNegPop = -1; 00291 *lgZeroPop = true; 00292 00293 for( level=1; level < nlev; level++ ) 00294 { 00295 pops[level] = 0.; 00296 /* these are off by one - lowest index is zero */ 00297 depart[level] = 0.; 00298 } 00299 00300 /* everything in ground */ 00301 pops[0] = abund; 00302 depart[0] = 1.; 00303 00304 /* there are two error exits from this routine, 00305 * previous one for zero abundance, and this one for zero excitation */ 00306 return; 00307 } 00308 00309 /* we will predict populations */ 00310 *lgZeroPop = false; 00311 00312 /* already have excitation pumping, now get deexcitation */ 00313 for( ilo=0; ilo < (nlev - 1); ilo++ ) 00314 { 00315 for( ihi=ilo + 1; ihi < nlev; ihi++ ) 00316 { 00317 /* (*AulPump)[low][ihi] is excitation rate, low -> ihi, due to external continuum, 00318 * so derive rate from upper to lower */ 00319 (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*g[ilo]/g[ihi]; 00320 } 00321 } 00322 00323 /* evaluate collision rates from collision strengths, but only if calling 00324 * routine has not already done this */ 00325 if( !lgCollRateDone ) 00326 { 00327 for( ilo=0; ilo < (nlev - 1); ilo++ ) 00328 { 00329 for( ihi=ilo + 1; ihi < nlev; ihi++ ) 00330 { 00331 /* this should be a collision strength */ 00332 ASSERT( (*col_str)[ihi][ilo]>= 0. ); 00333 /* this is deexcitation rate */ 00334 (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/g[ihi]*dense.cdsqte; 00335 /* this is excitation rate */ 00336 (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*g[ihi]/g[ilo]* 00337 excit[ilo][ihi]; 00338 } 00339 } 00340 } 00341 00342 /* rate equations equal zero */ 00343 amat.zero(); 00344 00345 /* following is column of vector - represents source terms from elsewhere, 00346 * if this is zero then matrix is singular and must replace one row with 00347 * population sum equation - if sum is non-zero then get total abundance 00348 * from source and sink terms */ 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 /* eqns for destruction of level 00360 * AulEscp[iho][ilo] are A*esc p, CollRate is coll excit in direction 00361 * AulPump[low][high] is excitation rate from low to high */ 00362 for( level=0; level < nlev; level++ ) 00363 { 00364 amat[level][level] = sink[level]; 00365 00366 /* leaving level to below */ 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 /* >>chng 97 jul 31, added pumping down 00372 * coming to level I from below */ 00373 amat[ilo][level] = -(*CollRate)[ilo][level] - (*AulPump)[ilo][level]; 00374 } 00375 00376 /* leaving level to above */ 00377 for( ihi=level + 1; ihi < nlev; ihi++ ) 00378 { 00379 amat[level][level] += (*CollRate)[level][ihi] + (*AulPump)[level][ihi]; 00380 /* coming to level from above */ 00381 amat[ihi][level] = -(*CollRate)[ihi][level] - (*AulEscp)[ihi][level] - 00382 (*AulDest)[ihi][level] - (*AulPump)[ihi][level]; 00383 /* >>chng 97 jul 31, added pumping down */ 00384 } 00385 } 00386 00387 /* homogeneous case, all source terms add up to zero, so use the population, 00388 * system has no total population information, 00389 * equation to replace redundant equation */ 00390 if( lgHomogeneous ) 00391 { 00392 for( level=0; level<nlev; ++level ) 00393 { 00394 amat[level][0] = 1.0; 00395 } 00396 /* these add up to total abundance */ 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 /* usage DGETRS, 'N' = no transpose 00429 * order of matrix, 00430 * number of cols in bvec, =1 00431 * array 00432 * leading dim of array */ 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 /* set populations */ 00442 for( level=0; level < nlev; level++ ) 00443 { 00444 /* save bvec into populations */ 00445 pops[level] = bvec[level]; 00446 } 00447 00448 /* now find total energy exchange rate, normally net cooling and its 00449 * temperature derivative */ 00450 *cooltl = 0.; 00451 *coolder = 0.; 00452 for( ihi=1; ihi < nlev; ihi++ ) 00453 { 00454 for( ilo=0; ilo < ihi; ilo++ ) 00455 { 00456 /* this is now net cooling rate [K cm-3 s-1] */ 00457 cool1 = (pops[ilo]*(*CollRate)[ilo][ihi] - pops[ihi]*(*CollRate)[ihi][ilo])* 00458 (ex[ihi] - ex[ilo]); 00459 *cooltl += cool1; 00460 00461 /* derivative wrt temperature - use Boltzmann factor relative to ground */ 00462 /* >>chng 03 aug 28, use real cool1 */ 00463 *coolder += cool1*( (ex[ihi] - ex[0])*thermal.tsq1 - thermal.halfte); 00464 } 00465 } 00466 /* convert from units of ex[] into ergs */ 00467 /* >>chng 05 dec 14, ex[] may be K or wn, TeConvFac will take care of either case */ 00468 *cooltl *= BOLTZMANN*TeConvFac; 00469 *coolder *= BOLTZMANN*TeConvFac; 00470 00471 /* fill in departure coefficients */ 00472 if( pops[0] > 1e-20 && excit[0][nlev-1] > 1e-20 ) 00473 { 00474 /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */ 00475 depart[0] = 1.; 00476 for( ihi=1; ihi < nlev; ihi++ ) 00477 { 00478 /* these are off by one - lowest index is zero */ 00479 depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit[0][ihi]; 00480 } 00481 } 00482 00483 else 00484 { 00485 /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */ 00486 for( ihi=0; ihi < nlev; ihi++ ) 00487 { 00488 /* Boltzmann factor or abundance too small, set departure coefficient to zero */ 00489 depart[ihi] = 0.; 00490 } 00491 depart[0] = 1.; 00492 } 00493 00494 /* sanity check for valid solution - non negative populations */ 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 /* these were reset to non zero values by the solver, but we will 00541 * assert that they are zero (for safety) when routine reenters so must 00542 * set to zero here, but only if asserts are in place */ 00543 for( ihi=1; ihi < nlev; ihi++ ) 00544 { 00545 for( ilo=0; ilo < ihi; ilo++ ) 00546 { 00547 /* zero ots destruction rate */ 00548 (*AulDest)[ilo][ihi] = 0.; 00549 /* both AulDest and AulPump (low, hi) are not used, should be zero */ 00550 (*AulPump)[ihi][ilo] = 0.; 00551 (*AulEscp)[ilo][ihi] = 0; 00552 } 00553 } 00554 # endif 00555 return; 00556 }