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 float 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 long int level,
00066 ihi,
00067 ilo,
00068 i,
00069 j;
00070
00071 bool lgHomogeneous;
00072
00073 double cool1,
00074 TeInverse,
00075 TeConvFac;
00076 double
00077 *bvec,
00078 *work,
00079 *amat,
00080 **excit,
00081 sum;
00082
00083 int32 ner;
00084 int32 *ipiv;
00085
00086 DEBUG_ENTRY( "atom_levelN()" );
00087
00088
00089 if( abund <= 0. )
00090 {
00091 *cooltl = 0.;
00092 *coolder = 0.;
00093
00094 *nNegPop = false;
00095 *lgZeroPop = true;
00096
00097 for( level=0; level < nlev; level++ )
00098 {
00099 pops[level] = 0.;
00100 depart[level] = 0.;
00101 }
00102
00103 depart[0] = 1.;
00104
00105
00106
00107 DEBUG_EXIT( "atom_levelN()" );
00108 return;
00109 }
00110
00111
00112 bvec = (double *)MALLOC( sizeof(double)*(size_t)(nlev+1) );
00113 work = (double *)MALLOC( sizeof(double)*(size_t)(nlev) );
00114 ipiv = (int32 *)MALLOC( sizeof(int32)*(size_t)(nlev) );
00115
00116 # ifdef AMAT
00117 # undef AMAT
00118 # endif
00119 # define AMAT(I_,J_) (*(amat+(I_)*(nlev)+(J_)))
00120
00121 amat=(double*)MALLOC(sizeof(double)*(size_t)(nlev*nlev));
00122
00123
00124 excit = ((double **)MALLOC((size_t)(nlev+1)*sizeof(double *)));
00125
00126 for( i=0; i<(nlev+1); ++i )
00127 {
00128 excit[i] = ((double *)MALLOC((size_t)(nlev+1)*sizeof(double )));
00129 }
00130
00131 # ifndef NDEBUG
00132
00133 ASSERT( ex[0] == 0. );
00134
00135 for( ihi=1; ihi < nlev; ihi++ )
00136 {
00137 for( ilo=0; ilo < ihi; ilo++ )
00138 {
00139
00140 ASSERT( (*AulDest)[ilo][ihi] == 0. );
00141 ASSERT( (*AulDest)[ihi][ilo] >= 0 );
00142
00143 ASSERT( (*AulPump)[ihi][ilo] == 0. );
00144 ASSERT( (*AulEscp)[ihi][ilo] >= 0 );
00145 ASSERT( (*AulEscp)[ilo][ihi] == 0 );
00146 ASSERT( (*col_str)[ihi][ilo] >= 0 );
00147 }
00148 }
00149 # endif
00150
00151 if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
00152 {
00153 fprintf( ioQQQ, " atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
00154 fprintf( ioQQQ, " AulDest\n" );
00155
00156 fprintf( ioQQQ, " hi lo" );
00157
00158 for( ilo=0; ilo < nlev-1; ilo++ )
00159 {
00160 fprintf( ioQQQ, "%4ld ", ilo );
00161 }
00162 fprintf( ioQQQ, " \n" );
00163
00164 for( ihi=1; ihi < nlev; ihi++ )
00165 {
00166 fprintf( ioQQQ, "%3ld", ihi );
00167 for( ilo=0; ilo < ihi; ilo++ )
00168 {
00169 fprintf( ioQQQ, "%10.2e", (*AulDest)[ihi][ilo] );
00170 }
00171 fprintf( ioQQQ, "\n" );
00172 }
00173
00174 fprintf( ioQQQ, " A*esc\n" );
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", (*AulEscp)[ihi][ilo] );
00188 }
00189 fprintf( ioQQQ, "\n" );
00190 }
00191
00192 fprintf( ioQQQ, " AulPump\n" );
00193
00194 fprintf( ioQQQ, " hi lo" );
00195 for( ilo=0; ilo < nlev-1; ilo++ )
00196 {
00197 fprintf( ioQQQ, "%4ld ", ilo );
00198 }
00199 fprintf( ioQQQ, " \n" );
00200
00201 for( ihi=1; ihi < nlev; ihi++ )
00202 {
00203 fprintf( ioQQQ, "%3ld", ihi );
00204 for( ilo=0; ilo < ihi; ilo++ )
00205 {
00206 fprintf( ioQQQ, "%10.2e", (*AulPump)[ilo][ihi] );
00207 }
00208 fprintf( ioQQQ, "\n" );
00209 }
00210
00211 fprintf( ioQQQ, " coll str\n" );
00212 fprintf( ioQQQ, " hi lo" );
00213 for( ilo=0; ilo < nlev-1; ilo++ )
00214 {
00215 fprintf( ioQQQ, "%4ld ", ilo );
00216 }
00217 fprintf( ioQQQ, " \n" );
00218
00219 for( ihi=1; ihi < nlev; ihi++ )
00220 {
00221 fprintf( ioQQQ, "%3ld", ihi );
00222 for( ilo=0; ilo < ihi; ilo++ )
00223 {
00224 fprintf( ioQQQ, "%10.2e", (*col_str)[ihi][ilo] );
00225 }
00226 fprintf( ioQQQ, "\n" );
00227 }
00228
00229 fprintf( ioQQQ, " coll rate\n" );
00230 fprintf( ioQQQ, " hi lo" );
00231 for( ilo=0; ilo < nlev-1; ilo++ )
00232 {
00233 fprintf( ioQQQ, "%4ld ", ilo );
00234 }
00235 fprintf( ioQQQ, " \n" );
00236
00237 for( ihi=1; ihi < nlev; ihi++ )
00238 {
00239 fprintf( ioQQQ, "%3ld", ihi );
00240 for( ilo=0; ilo < ihi; ilo++ )
00241 {
00242 fprintf( ioQQQ, "%10.2e", (*CollRate)[ihi][ilo] );
00243 }
00244 fprintf( ioQQQ, "\n" );
00245 }
00246 }
00247
00248
00249 if( chExUnits=='K' )
00250 {
00251
00252
00253 TeInverse = 1./phycon.te;
00254
00255 TeConvFac = 1.;
00256 }
00257 else if( chExUnits=='w' )
00258 {
00259
00260 TeInverse = 1./phycon.te_wn;
00261 TeConvFac = T1CM;
00262 }
00263 else
00264 TotalInsanity();
00265
00266
00267 for( ilo=0; ilo < (nlev - 1); ilo++ )
00268 {
00269 for( ihi=ilo + 1; ihi < nlev; ihi++ )
00270 {
00271
00272
00273 excit[ilo][ihi] = sexp((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=1; ilo <= nlev; ilo++ )
00283 {
00284 fprintf( ioQQQ, "%4ld ", ilo );
00285 }
00286 fprintf( ioQQQ, " \n" );
00287
00288 for( ihi=0; ihi < nlev; ihi++ )
00289 {
00290 fprintf( ioQQQ, "%3ld", ihi );
00291 for( ilo=0; ilo < nlev; ilo++ )
00292 {
00293 fprintf( ioQQQ, "%10.2e", excit[ilo][ihi] );
00294 }
00295 fprintf( ioQQQ, "\n" );
00296 }
00297 }
00298
00299
00300
00301 if( (excit[0][nlev-1] + (*AulPump)[0][nlev-1] < 1e-13 ) && (source[nlev-1]==0.) )
00302 {
00303 *cooltl = 0.;
00304 *coolder = 0.;
00305
00306
00307 *nNegPop = -1;
00308 *lgZeroPop = true;
00309
00310 for( level=1; level < nlev; level++ )
00311 {
00312 pops[level] = 0.;
00313
00314 depart[level] = 0.;
00315 }
00316
00317
00318 pops[0] = abund;
00319 depart[0] = 1.;
00320
00321 free( bvec );
00322 free( work );
00323 free( ipiv );
00324 free( amat );
00325 for( i=0; i < (nlev+1); ++i )
00326 {
00327 free( excit[i] );
00328 }
00329 free( excit );
00330
00331
00332
00333 DEBUG_EXIT( "atom_levelN()" );
00334 return;
00335 }
00336
00337
00338 *lgZeroPop = false;
00339
00340
00341 for( ilo=0; ilo < (nlev - 1); ilo++ )
00342 {
00343 for( ihi=ilo + 1; ihi < nlev; ihi++ )
00344 {
00345
00346
00347 (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*g[ilo]/g[ihi];
00348 }
00349 }
00350
00351
00352
00353 if( !lgCollRateDone )
00354 {
00355
00356 for( ilo=0; ilo < (nlev - 1); ilo++ )
00357 {
00358 for( ihi=ilo + 1; ihi < nlev; ihi++ )
00359 {
00360
00361 ASSERT( (*col_str)[ihi][ilo]>= 0. );
00362
00363
00364
00365 (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/g[ihi]*dense.cdsqte;
00366
00367 (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
00368 excit[ilo][ihi];
00369 }
00370 }
00371 }
00372
00373
00374 for( level=0; level < nlev; level++ )
00375 {
00376 for( j=0; j < nlev; j++ )
00377 {
00378 AMAT(level,j) = 0.;
00379 }
00380 }
00381
00382
00383
00384
00385
00386 sum = 0.;
00387 lgHomogeneous = false;
00388 for( level=0; level < nlev; level++ )
00389 {
00390 bvec[level] = source[level];
00391 sum += bvec[level];
00392 }
00393 if( sum==0. )
00394 lgHomogeneous = true;
00395
00396
00397
00398
00399 for( level=0; level < nlev; level++ )
00400 {
00401 AMAT(level,level) = sink[level];
00402
00403
00404 for( ilo=0; ilo < level; ilo++ )
00405 {
00406
00407 AMAT(level,level) += (*CollRate)[level][ilo] + (*AulEscp)[level][ilo] + (*AulDest)[level][ilo] +
00408 (*AulPump)[level][ilo];
00409
00410
00411 AMAT(ilo,level) = -(*CollRate)[ilo][level] - (*AulPump)[ilo][level];
00412 }
00413
00414
00415 for( ihi=level + 1; ihi < nlev; ihi++ )
00416 {
00417 AMAT(level,level) += (*CollRate)[level][ihi] + (*AulPump)[level][ihi];
00418
00419 AMAT(ihi,level) = -(*CollRate)[ihi][level] - (*AulEscp)[ihi][level] - (*AulDest)[ihi][level] -
00420 (*AulPump)[ihi][level];
00421
00422 }
00423 }
00424
00425
00426
00427
00428 if( lgHomogeneous )
00429 {
00430 for( level=0; level<nlev; ++level )
00431 {
00432 AMAT(level,0) = 1.0;
00433 }
00434
00435 bvec[0] = abund;
00436 }
00437
00438 if( lgDeBug )
00439 {
00440 fprintf(ioQQQ," AMAT matrix follows:\n");
00441 for( level=0; level < nlev; level++ )
00442 {
00443 for( j=0; j < nlev; j++ )
00444 {
00445 fprintf(ioQQQ," %.4e" , AMAT(level,j));
00446 }
00447 fprintf(ioQQQ,"\n");
00448 }
00449 if( sum==0. )
00450 {
00451 fprintf(ioQQQ," Sum creation zero so pop sum used\n");
00452 }
00453 else
00454 {
00455 fprintf(ioQQQ," Sum creation non-zero (%e), vector follows:\n",sum);
00456 for( j=0; j < nlev; j++ )
00457 {
00458 fprintf(ioQQQ," %.4e" , bvec[j] );
00459 }
00460 fprintf(ioQQQ,"\n");
00461 }
00462 }
00463
00464
00465
00466 ner = 0;
00467 getrf_wrapper(nlev,nlev,amat,nlev,ipiv,&ner);
00468
00469
00470
00471
00472
00473
00474
00475
00476 getrs_wrapper('N',nlev,1,amat,nlev,ipiv,bvec,nlev,&ner);
00477
00478 if( ner != 0 )
00479 {
00480 fprintf( ioQQQ, " atom_levelN: dgetrs finds singular or ill-conditioned matrix\n" );
00481 puts( "[Stop in leveln]" );
00482 cdEXIT(EXIT_FAILURE);
00483 }
00484
00485
00486 for( level=0; level < nlev; level++ )
00487 {
00488
00489 pops[level] = bvec[level];
00490 }
00491
00492
00493
00494 *cooltl = 0.;
00495 *coolder = 0.;
00496 for( ihi=1; ihi < nlev; ihi++ )
00497 {
00498 for( ilo=0; ilo < ihi; ilo++ )
00499 {
00500
00501 cool1 = (pops[ilo]*(*CollRate)[ilo][ihi] - pops[ihi]*(*CollRate)[ihi][ilo])*
00502 (ex[ihi] - ex[ilo]);
00503 *cooltl += cool1;
00504
00505
00506
00507
00508 *coolder += cool1*( (ex[ihi] - ex[0])*thermal.tsq1 - thermal.halfte);
00509 }
00510 }
00511
00512
00513 *cooltl *= BOLTZMANN*TeConvFac;
00514 *coolder *= BOLTZMANN*TeConvFac;
00515
00516
00517 if( pops[0] > 1e-20 && excit[0][nlev-1] > 1e-20 )
00518 {
00519
00520 depart[0] = 1.;
00521 for( ihi=1; ihi < nlev; ihi++ )
00522 {
00523
00524 depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/
00525 excit[0][ihi];
00526 }
00527 }
00528
00529 else
00530 {
00531
00532 for( ihi=0; ihi < nlev; ihi++ )
00533 {
00534
00535 depart[ihi] = 0.;
00536 }
00537 depart[0] = 1.;
00538 }
00539
00540
00541 *nNegPop = false;
00542 for( level=0; level < nlev; level++ )
00543 {
00544 if( pops[level] < 0. )
00545 {
00546 *nNegPop = true;
00547 }
00548 }
00549
00550 if( *nNegPop )
00551 {
00552 fprintf( ioQQQ, "\n PROBLEM atom_levelN found negative population, atom=%s\n",
00553 chLabel );
00554 fprintf( ioQQQ, " absolute population=" );
00555
00556 for( level=0; level < nlev; level++ )
00557 {
00558 fprintf( ioQQQ, "%10.2e", pops[level] );
00559 }
00560
00561 fprintf( ioQQQ, "\n" );
00562 for( level=0; level < nlev; level++ )
00563 {
00564 pops[level] = (double)MAX2(0.,pops[level]);
00565 }
00566 }
00567
00568 if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
00569 {
00570 fprintf( ioQQQ, "\n atom_leveln absolute population " );
00571 for( level=0; level < nlev; level++ )
00572 {
00573 fprintf( ioQQQ, " %10.2e", pops[level] );
00574 }
00575 fprintf( ioQQQ, "\n" );
00576
00577 fprintf( ioQQQ, " departure coefficients" );
00578 for( level=0; level < nlev; level++ )
00579 {
00580 fprintf( ioQQQ, " %10.2e", depart[level] );
00581 }
00582 fprintf( ioQQQ, "\n\n" );
00583 }
00584
00585 free( bvec );
00586 free( work );
00587 free( ipiv );
00588 free( amat );
00589 for( i=0; i < (nlev+1); ++i )
00590 {
00591 free( excit[i] );
00592 }
00593 free( excit );
00594
00595 # ifndef NDEBUG
00596
00597
00598
00599 for( ihi=1; ihi < nlev; ihi++ )
00600 {
00601 for( ilo=0; ilo < ihi; ilo++ )
00602 {
00603
00604 (*AulDest)[ilo][ihi] = 0.;
00605
00606 (*AulPump)[ihi][ilo] = 0.;
00607 (*AulEscp)[ilo][ihi] = 0;
00608 }
00609 }
00610 # endif
00611
00612 DEBUG_EXIT( "atom_levelN()" );
00613 return;
00614
00615 }
00616