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