00001
00002
00003
00004 #include "cddefines.h"
00005 #include "phycon.h"
00006 #include "physconst.h"
00007 #include "dense.h"
00008 #include "transition.h"
00009 #include "thermal.h"
00010 #include "cooling.h"
00011
00012 #include "atoms.h"
00013 #include "rfield.h"
00014
00015 void atom_level3(const TransitionProxy& t10,
00016 const TransitionProxy& t21,
00017 const TransitionProxy& t20)
00018 {
00019 char chLab[5],
00020 chLab10[11];
00021 double AbunxIon,
00022 a,
00023 a10,
00024 a20,
00025 a21,
00026 b,
00027 beta,
00028 bolt01,
00029 bolt02,
00030 bolt12,
00031 c,
00032 ener10,
00033 ener20,
00034 ener21,
00035 g0,
00036 g010,
00037 g020,
00038 g1,
00039 g110,
00040 g121,
00041 g2,
00042 g220,
00043 g221,
00044 o10,
00045 o20,
00046 o21,
00047 p0,
00048 p1,
00049 p2,
00050 pump01,
00051 pump02,
00052 pump12;
00053
00054 int nLev3Fail;
00055
00056 double TotCool,
00057 TotHeat,
00058 TotInten,
00059 alpha,
00060 alpha1,
00061 alpha2,
00062 c01,
00063 c02,
00064 c10,
00065 c12,
00066 c20,
00067 c21,
00068 cnet01,
00069 cnet02,
00070 cnet12,
00071 cool01,
00072 cool02,
00073 cool12,
00074 heat10,
00075 heat20,
00076 heat21,
00077 hnet01,
00078 hnet02,
00079 hnet12,
00080 pump10,
00081 pump20,
00082 pump21,
00083 r01,
00084 r02,
00085 r10,
00086 r12,
00087 r20,
00088 r21,
00089 temp01,
00090 temp02,
00091 temp12;
00092
00093 DEBUG_ENTRY( "atom_level3()" );
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 g010 = (*t10.Lo()).g();
00106 g110 = (*t10.Hi()).g();
00107
00108 g121 = (*t21.Lo()).g();
00109 g221 = (*t21.Hi()).g();
00110
00111 g020 = (*t20.Lo()).g();
00112 g220 = (*t20.Hi()).g();
00113
00114
00115 if( g010 > 0. )
00116 {
00117 g0 = g010;
00118 }
00119
00120 else if( g020 > 0. )
00121 {
00122 g0 = g020;
00123 }
00124
00125 else
00126 {
00127 g0 = -1.;
00128 strcpy( chLab10, chLineLbl(t10) );
00129 fprintf( ioQQQ, " PROBLEM atom_level3: insane stat weights g0 :%10.2e%10.2e %10.10s\n",
00130 g010, g020, chLab10 );
00131 TotalInsanity();
00132 }
00133
00134 if( g110 > 0. )
00135 {
00136 g1 = g110;
00137 }
00138
00139 else if( g121 > 0. )
00140 {
00141 g1 = g121;
00142 }
00143
00144 else
00145 {
00146 g1 = -1.;
00147 strcpy( chLab10, chLineLbl(t10) );
00148 fprintf( ioQQQ, " PROBLEM atom_level3: insane stat weights g1 :%10.2e%10.2e %10.10s\n",
00149 g010, g020, chLab );
00150 TotalInsanity();
00151 }
00152
00153 if( g220 > 0. )
00154 {
00155 g2 = g220;
00156 }
00157 else if( g221 > 0. )
00158 {
00159 g2 = g221;
00160 }
00161
00162 else
00163 {
00164 g2 = -1.;
00165 strcpy( chLab10, chLineLbl(t20) );
00166 fprintf( ioQQQ, " PROBLEM atom_level3: insane stat weights g2 :%10.2e%10.2e %10.10s\n",
00167 g010, g020, chLab10 );
00168 TotalInsanity();
00169 }
00170
00171
00172
00173 if( g010 > 0. )
00174 {
00175
00176 chIonLbl(chLab, t10);
00177 AbunxIon = dense.xIonDense[ (*t10.Hi()).nelem() -1][(*t10.Hi()).IonStg()-1];
00178 }
00179
00180 else if( g121 > 0. )
00181 {
00182
00183 chIonLbl(chLab, t21);
00184 AbunxIon = dense.xIonDense[(*t21.Hi()).nelem() -1][(*t21.Hi()).IonStg()-1];
00185 }
00186
00187 else
00188
00189 {
00190 chLab[0] = ' ';
00191 AbunxIon = 0.;
00192 fprintf( ioQQQ, " PROBLEM atom_level3: insanity at g010 g121 branch \n" );
00193 TotalInsanity();
00194 }
00195
00196 a = t10.EnergyK()*phycon.teinv;
00197 b = t21.EnergyK()*phycon.teinv;
00198 c = t20.EnergyK()*phycon.teinv;
00199
00200 if( c == 0. )
00201 {
00202 c = a + b;
00203 }
00204
00205
00206
00207 nLev3Fail = -1;
00208
00209
00210 if( ( t10.EnergyErg()/EN1RYD < rfield.plsfrq && t20.EnergyErg()/EN1RYD < rfield.plsfrq ) )
00211 {
00212 atom_level2( t21 );
00213 return;
00214 }
00215 else if( ( t10.EnergyErg()/EN1RYD < rfield.plsfrq && t21.EnergyErg()/EN1RYD < rfield.plsfrq ) )
00216 {
00217 atom_level2( t20 );
00218 return;
00219 }
00220 else if( ( t20.EnergyErg()/EN1RYD < rfield.plsfrq && t21.EnergyErg()/EN1RYD < rfield.plsfrq ) )
00221 {
00222 atom_level2( t10 );
00223 return;
00224 }
00225
00228 if( AbunxIon <= 1e-30 || c > 60. )
00229 {
00230 nLev3Fail = 0;
00231
00232
00233 atoms.PopLevels[0] = AbunxIon;
00234 atoms.PopLevels[1] = 0.;
00235 atoms.PopLevels[2] = 0.;
00236
00238 atoms.DepLTELevels[0] = 1.;
00239 atoms.DepLTELevels[1] = 0.;
00240 atoms.DepLTELevels[2] = 0.;
00241
00242
00243 t21.Emis().PopOpc() = 0.;
00244 t10.Emis().PopOpc() = AbunxIon;
00245 t20.Emis().PopOpc() = AbunxIon;
00246 (*t21.Lo()).Pop() = 0.;
00247 (*t10.Lo()).Pop() = AbunxIon;
00248 (*t20.Lo()).Pop() = AbunxIon;
00249 (*t21.Hi()).Pop() = 0.;
00250 (*t10.Hi()).Pop() = 0.;
00251 (*t20.Hi()).Pop() = 0.;
00252
00253
00254 t20.Coll().heat() = 0.;
00255 t21.Coll().heat() = 0.;
00256 t10.Coll().heat() = 0.;
00257
00258
00259 t21.Emis().xIntensity() = 0.;
00260 t10.Emis().xIntensity() = 0.;
00261 t20.Emis().xIntensity() = 0.;
00262
00263
00264 t20.Coll().cool() = 0.;
00265 t21.Coll().cool() = 0.;
00266 t10.Coll().cool() = 0.;
00267
00268
00269 # if 0
00270
00271 t20.ots() = 0.;
00272 t21.ots() = 0.;
00273 t10.ots() = 0.;
00274 # endif
00275
00276
00277 t21.Emis().phots() = 0.;
00278 t10.Emis().phots() = 0.;
00279 t20.Emis().phots() = 0.;
00280
00281
00282 t21.Emis().ColOvTot() = 0.;
00283 t10.Emis().ColOvTot() = 0.;
00284 t20.Emis().ColOvTot() = 0.;
00285
00286
00287 CoolAdd(chLab, t21.WLAng() ,0.);
00288 CoolAdd(chLab, t10.WLAng() ,0.);
00289 CoolAdd(chLab, t20.WLAng() ,0.);
00290 return;
00291 }
00292
00293
00294 o10 = t10.Coll().col_str();
00295 o21 = t21.Coll().col_str();
00296 o20 = t20.Coll().col_str();
00297
00298
00299
00300 ASSERT( g010*g020 == 0. || fp_equal( g010, g020 ) );
00301
00302 ASSERT( g110*g121 == 0. || fp_equal( g110, g121 ) );
00303
00304 ASSERT( g221*g220 == 0. || fp_equal( g221, g220 ) );
00305
00306
00307
00308 ASSERT( ((*t10.Hi()).IonStg()*(*t21.Hi()).IonStg() == 0) || ((*t10.Hi()).IonStg() == (*t21.Hi()).IonStg()));
00309
00310 ASSERT( ((*t20.Hi()).IonStg()*(*t21.Hi()).IonStg() == 0) || ((*t20.Hi()).IonStg() == (*t21.Hi()).IonStg() ) );
00311
00312 ASSERT( ((*t10.Hi()).nelem() * (*t21.Hi()).nelem() == 0) || ((*t10.Hi()).nelem() == (*t21.Hi()).nelem()) );
00313
00314 ASSERT( ((*t20.Hi()).nelem() * (*t21.Hi()).nelem() == 0) || ((*t20.Hi()).nelem() == (*t21.Hi()).nelem()) );
00315
00316 ASSERT( o10 > 0. && o21 > 0. && o20 > 0. );
00317
00318
00319
00320
00321 a21 = t21.Emis().Aul() * (t21.Emis().Pesc()+ t21.Emis().Pelec_esc() + t21.Emis().Pdest());
00322 a10 = t10.Emis().Aul() * (t10.Emis().Pesc()+ t10.Emis().Pelec_esc() + t10.Emis().Pdest());
00323 a20 = t20.Emis().Aul() * (t20.Emis().Pesc()+ t20.Emis().Pelec_esc() + t20.Emis().Pdest());
00324
00325
00326
00327 if( t10.Emis().Aul() == 0. )
00328 {
00329 ener20 = t20.EnergyErg();
00330 ener21 = t21.EnergyErg();
00331 ener10 = ener20 - ener21;
00332 bolt12 = exp(-t21.EnergyK()/phycon.te);
00333 bolt02 = exp(-t20.EnergyK()/phycon.te);
00334 bolt01 = bolt02/bolt12;
00335 temp12 = t21.EnergyK();
00336 temp02 = t20.EnergyK();
00337 temp01 = temp02 - temp12;
00338 }
00339
00340 else if( t21.Emis().Aul() == 0. )
00341 {
00342 ener10 = t10.EnergyErg();
00343 ener20 = t20.EnergyErg();
00344 ener21 = ener20 - ener10;
00345 bolt01 = exp(-t10.EnergyK()/phycon.te);
00346 bolt02 = exp(-t20.EnergyK()/phycon.te);
00347 bolt12 = bolt02/bolt01;
00348 temp02 = t20.EnergyK();
00349 temp01 = t10.EnergyK();
00350 temp12 = temp02 - temp01;
00351 }
00352
00353 else if( t20.Emis().Aul() == 0. )
00354 {
00355 ener10 = t10.EnergyErg();
00356 ener21 = t21.EnergyErg();
00357 ener20 = ener21 + ener10;
00358 bolt01 = exp(-t10.EnergyK()/phycon.te);
00359 bolt12 = exp(-t21.EnergyK()/phycon.te);
00360 bolt02 = bolt01*bolt12;
00361 temp01 = t10.EnergyK();
00362 temp12 = t21.EnergyK();
00363 temp02 = temp01 + temp12;
00364 }
00365
00366 else
00367 {
00368
00369 ener10 = t10.EnergyErg();
00370 ener20 = t20.EnergyErg();
00371 ener21 = t21.EnergyErg();
00372 bolt01 = exp(-t10.EnergyK()/phycon.te);
00373 bolt12 = exp(-t21.EnergyK()/phycon.te);
00374 bolt02 = bolt01*bolt12;
00375 temp02 = t20.EnergyK();
00376 temp01 = t10.EnergyK();
00377 temp12 = t21.EnergyK();
00378 }
00379
00380
00381 ASSERT( ener10 > 0. && ener20 > 0. && ener21 > 0. );
00382
00383
00384 ASSERT( ener10 < ener20 && ener21 < ener20 );
00385
00386
00387 ASSERT( fabs((ener10+ener21)/ener20-1.) < 1e-4 );
00388
00389 pump01 = t10.Emis().pump();
00390 pump10 = pump01*g0/g1;
00391 pump12 = t21.Emis().pump();
00392 pump21 = pump12*g1/g2;
00393 pump02 = t20.Emis().pump();
00394 pump20 = pump02*g0/g2;
00395
00396
00397 c01 = o10*bolt01*dense.cdsqte/g0;
00398 r01 = c01 + pump01;
00399 c10 = o10*dense.cdsqte/g1;
00400 r10 = c10 + a10 + pump10;
00401 c20 = o20*dense.cdsqte/g2;
00402 r20 = c20 + a20 + pump20;
00403 c02 = o20*bolt02*dense.cdsqte/g0;
00404 r02 = c02 + pump02;
00405 c12 = o21*bolt12*dense.cdsqte/g1;
00406 r12 = c12 + pump12;
00407 c21 = o21*dense.cdsqte/g2;
00408 r21 = c21 + a21 + pump21;
00409
00410 alpha1 = (double)(AbunxIon)*(r01+r02)/(r10+r01+r02);
00411 alpha2 = (double)(AbunxIon)*(r01)/(r10+r12+r01);
00412 alpha = alpha1 - alpha2;
00413
00414
00415
00416 beta = (r21 - r01)/(r10 + r12 + r01) + (r20 + r01 + r02)/(r10 +
00417 r01 + r02);
00418
00419 if( alpha/MAX2(alpha1,alpha2) < 1e-11 )
00420 {
00421
00422 p2 = 0.;
00423 alpha = 0.;
00424 nLev3Fail = 1;
00425 }
00426
00427 else
00428 {
00429 p2 = alpha/beta;
00430 }
00431 atoms.PopLevels[2] = p2;
00432
00433 if( alpha < 0. || beta < 0. )
00434 {
00435 fprintf( ioQQQ, " atom_level3: insane n2 pop alf, bet, p2=%10.2e%10.2e%10.2e %10.10s t=%10.2e\n",
00436 alpha, beta, p2, chLab, phycon.te );
00437 fprintf( ioQQQ, " gs are%5.1f%5.1f%5.1f\n", g0, g1,
00438 g2 );
00439 fprintf( ioQQQ, " Bolts are%10.2e%10.2e%10.2e\n",
00440 bolt01, bolt12, bolt02 );
00441 fprintf( ioQQQ, " As are%10.2e%10.2e%10.2e\n", a10,
00442 a21, a20 );
00443 fprintf( ioQQQ, " Energies are%10.2e%10.2e%10.2e\n",
00444 ener10, ener21, ener20 );
00445 fprintf( ioQQQ, " 2 terms, dif of alpha are%15.6e%15.6e\n",
00446 (r01 + r02)/(r10 + r01 + r02), r01/(r10 + r12 + r01) );
00447 ShowMe();
00448 cdEXIT(EXIT_FAILURE);
00449 }
00450
00451 alpha = (double)(AbunxIon)*(r01+r02) - (double)(p2)*(r20+r01+r02);
00452
00453
00454
00455 if( fabs(alpha)/(MAX2(AbunxIon*(r01+r02),p2*(r20+r01+r02))) < 1e-9 )
00456 {
00457 alpha = 0.;
00458 nLev3Fail = 2;
00459 }
00460
00461 beta = r10 + r01 + r02;
00462 p1 = alpha/beta;
00463 atoms.PopLevels[1] = p1;
00464
00465 if( p1 < 0. )
00466 {
00467 if( p1 > -1e-37 )
00468 {
00469
00470 p1 = 0.;
00471 atoms.PopLevels[1] = p1;
00472 nLev3Fail = 3;
00473 }
00474
00475 else
00476 {
00477
00478 fprintf( ioQQQ, " atom_level3: insane n1 pop alf, bet, p1=%10.2e%10.2e%10.2e %10.10s%5f\n",
00479 alpha, beta, p1, chLab, t10.WLAng() );
00480 fprintf( ioQQQ, " local electron density and temperature were%10.2e%10.2e\n",
00481 dense.eden, phycon.te );
00482 ShowMe();
00483 cdEXIT(EXIT_FAILURE);
00484 }
00485 }
00486
00487 p0 = AbunxIon - p2 - p1;
00488
00489
00490 atoms.PopLevels[0] = p0;
00491 if( p0 <= 0. )
00492 {
00493 fprintf( ioQQQ, " atom_level3: insane n0 pop p1, 2, abun=%10.2e%10.2e%10.2e \n",
00494 p1, p2, AbunxIon );
00495 ShowMe();
00496 cdEXIT(EXIT_FAILURE);
00497 }
00498
00499
00500 (*t21.Lo()).Pop() = p1;
00501 (*t10.Lo()).Pop() = p0;
00502 (*t20.Lo()).Pop() = p0;
00503 t21.Emis().PopOpc() = (p1 - p2*g1/g2);
00504 t10.Emis().PopOpc() = (p0 - p1*g0/g1);
00505 t20.Emis().PopOpc() = (p0 - p2*g0/g2);
00506 (*t21.Hi()).Pop() = p2;
00507 (*t10.Hi()).Pop() = p1;
00508 (*t20.Hi()).Pop() = p2;
00509
00510
00511 t21.Emis().phots() = t21.Emis().Aul() * (t21.Emis().Pesc() + t21.Emis().Pelec_esc())*p2;
00512 t21.Emis().xIntensity() = t21.Emis().phots() * t21.EnergyErg();
00513
00514 t20.Emis().phots() = t20.Emis().Aul() * (t20.Emis().Pesc() + t20.Emis().Pelec_esc())*p2;
00515 t20.Emis().xIntensity() = t20.Emis().phots() * t20.EnergyErg();
00516
00517 t10.Emis().phots() = t10.Emis().Aul() * (t10.Emis().Pesc() + t10.Emis().Pelec_esc())*p1;
00518 t10.Emis().xIntensity() = t10.Emis().phots() * t10.EnergyErg();
00519
00520 # if 0
00521
00522 t21.ots() = p2 * t21.Emis().Aul() * t21.Emis().Pdest();
00523 t20.ots() = p2 * t20.Emis().Aul() * t20.Emis().Pdest();
00524 t10.ots() = p2 * t10.Emis().Aul() * t10.Emis().Pdest();
00525
00526 RT_OTS_AddLine( t21.ots() , t21.ipCont);
00527 RT_OTS_AddLine( t20.ots() , t20.ipCont);
00528 RT_OTS_AddLine( t10.ots() , t10.ipCont);
00529 # endif
00530
00531
00532
00533
00534
00535 TotInten = t21.Emis().phots() * (double)t21.EnergyErg()
00536 + t20.Emis().phots() * (double)t20.EnergyErg() +
00537 t10.Emis().phots() * (double)t10.EnergyErg();
00538
00539
00540 if( r12 > 0. )
00541 {
00542 t21.Emis().ColOvTot() = c12/r12;
00543 }
00544 else
00545 {
00546 t21.Emis().ColOvTot() = 0.;
00547 }
00548
00549 if( r01 > 0. )
00550 {
00551 t10.Emis().ColOvTot() = c01/r01;
00552 }
00553 else
00554 {
00555 t10.Emis().ColOvTot() = 0.;
00556 }
00557
00558 if( r02 > 0. )
00559 {
00560 t20.Emis().ColOvTot() = c02/r02;
00561 }
00562 else
00563 {
00564 t20.Emis().ColOvTot() = 0.;
00565 }
00566
00567
00568 heat20 = p2*c20*ener20;
00569 cool02 = p0*c02*ener20;
00570 heat21 = p2*c21*ener21;
00571 cool12 = p1*c12*ener21;
00572 heat10 = p1*c10*ener10;
00573 cool01 = p0*c01*ener10;
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583 cnet02 = cool02 - heat20*t20.Emis().ColOvTot();
00584 hnet02 = heat20 *(1. - t20.Emis().ColOvTot());
00585 cnet12 = cool12 - heat21*t21.Emis().ColOvTot();
00586 hnet12 = heat21 *(1. - t21.Emis().ColOvTot());
00587 cnet01 = cool01 - heat10*t10.Emis().ColOvTot();
00588 hnet01 = heat10 *(1. - t10.Emis().ColOvTot());
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 if( fabs(cnet01/MAX2(DBL_MIN,cool01)) < 1e-10 )
00608 {
00609 nLev3Fail = 4;
00610 cnet02 = 0.;
00611 hnet02 = 0.;
00612 cnet12 = 0.;
00613 hnet12 = 0.;
00614 cnet01 = 0.;
00615 hnet01 = 0.;
00616 }
00617
00618 TotCool = cnet02 + cnet12 + cnet01;
00619 TotHeat = hnet02 + hnet12 + hnet01;
00620
00621
00622 if( TotInten > 0. )
00623 {
00624 cool02 = TotCool * t20.Emis().phots() * (double)t20.EnergyErg()/TotInten;
00625 cool12 = TotCool * t21.Emis().phots() * (double)t21.EnergyErg()/TotInten;
00626 cool01 = TotCool * t10.Emis().phots() * (double)t10.EnergyErg()/TotInten;
00627 heat20 = TotHeat * t20.Emis().phots() * (double)t20.EnergyErg()/TotInten;
00628 heat21 = TotHeat * t21.Emis().phots() * (double)t21.EnergyErg()/TotInten;
00629 heat10 = TotHeat * t10.Emis().phots() * (double)t10.EnergyErg()/TotInten;
00630 t20.Coll().cool() = cool02;
00631 t21.Coll().cool() = cool12;
00632 t10.Coll().cool() = cool01;
00633 t20.Coll().heat() = heat20;
00634 t21.Coll().heat() = heat21;
00635 t10.Coll().heat() = heat10;
00636 }
00637 else
00638 {
00639 nLev3Fail = 5;
00640 cool02 = 0.;
00641 cool12 = 0.;
00642 cool01 = 0.;
00643 heat20 = 0.;
00644 heat21 = 0.;
00645 heat10 = 0.;
00646 t20.Coll().cool() = 0.;
00647 t21.Coll().cool() = 0.;
00648 t10.Coll().cool() = 0.;
00649 t20.Coll().heat() = 0.;
00650 t21.Coll().heat() = 0.;
00651 t10.Coll().heat() = 0.;
00652 }
00653
00654
00655
00656
00657 CoolAdd(chLab, t21.WLAng() ,cool12);
00658 CoolAdd(chLab, t10.WLAng() ,cool01);
00659 CoolAdd(chLab, t20.WLAng() ,cool02);
00660
00661
00662
00663
00664
00665
00666 thermal.dCooldT += t10.Coll().cool()*(temp01*thermal.tsq1 - thermal.halfte) +
00667 (t20.Coll().cool() + t21.Coll().cool())*(temp02*thermal.tsq1 - thermal.halfte);
00668
00669
00670 {
00671 enum{DEBUG_LOC=false};
00672 if( DEBUG_LOC )
00673 {
00674 fprintf(ioQQQ,"atom_level3 nLev3Fail %i\n",
00675 nLev3Fail );
00676 }
00677 }
00678 return;
00679 }