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_level3 compute three level atom, 10, 21, and 20 are line */ 00004 #include "cddefines.h" 00005 #include "phycon.h" 00006 #include "physconst.h" 00007 #include "dense.h" 00008 #include "lines_service.h" 00009 #include "thermal.h" 00010 #include "cooling.h" 00011 /* NB - use the nLev3Fail failures mode indicators when debugging this routine */ 00012 #include "atoms.h" 00013 #include "rfield.h" 00014 00015 void atom_level3(transition * t10, 00016 transition * t21, 00017 transition * 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 /* compute three level atom, 10, 21, and 20 are line 00096 * arrays for 10, 21, and 20 transitions. 00097 * one can be a dummy */ 00098 /* >>chng 96 dec 06, to double precision due to round off problems below */ 00099 00100 /* generalized three level atom for any ion 00101 * sum of three levels normalized to total abundance 00102 * 00103 * stat weights of all three lines 00104 * sanity check will confirm ok */ 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 /* these are statistical weights */ 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 /* abundances from the elements grid 00172 * one of these must be a true line */ 00173 if( g010 > 0. ) 00174 { 00175 /* put null terminated line label into chLab using line array*/ 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 /* put null terminated line label into chLab using line array*/ 00183 chIonLbl(chLab, t21); 00184 AbunxIon = dense.xIonDense[t21->Hi->nelem -1][t21->Hi->IonStg-1]; 00185 } 00186 00187 else 00188 /* this cannot possibly happen */ 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 /* if still neg at end, then success!, so possible to 00206 * to check why zero returned */ 00207 nLev3Fail = -1; 00208 00209 /* if two of the lines are below the plasma frequency, hand the third in atom_level2 */ 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 /* all populations are zero */ 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 /* level populations */ 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 /* line heating */ 00254 t20->Coll.heat = 0.; 00255 t21->Coll.heat = 0.; 00256 t10->Coll.heat = 0.; 00257 00258 /* intensity of line */ 00259 t21->Emis->xIntensity = 0.; 00260 t10->Emis->xIntensity = 0.; 00261 t20->Emis->xIntensity = 0.; 00262 00263 /* line cooling */ 00264 t20->Coll.cool = 0.; 00265 t21->Coll.cool = 0.; 00266 t10->Coll.cool = 0.; 00267 00268 /* local ots rates */ 00269 # if 0 00270 /* >>chng 03 oct 04, move to RT_OTS */ 00271 t20->ots = 0.; 00272 t21->ots = 0.; 00273 t10->ots = 0.; 00274 # endif 00275 00276 /* number of photons in line zero */ 00277 t21->Emis->phots = 0.; 00278 t10->Emis->phots = 0.; 00279 t20->Emis->phots = 0.; 00280 00281 /* ratio coll over total excitation */ 00282 t21->Emis->ColOvTot = 0.; 00283 t10->Emis->ColOvTot = 0.; 00284 t20->Emis->ColOvTot = 0.; 00285 00286 /* add zero to cooling */ 00287 CoolAdd(chLab, t21->WLAng ,0.); 00288 CoolAdd(chLab, t10->WLAng ,0.); 00289 CoolAdd(chLab, t20->WLAng ,0.); 00290 return; 00291 } 00292 00293 /* collision strengths */ 00294 o10 = t10->Coll.cs; 00295 o21 = t21->Coll.cs; 00296 o20 = t20->Coll.cs; 00297 00298 /* begin sanity checks, check statistic weights, 00299 * first check is protection against dummy lines */ 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 /* both abundances must be same, equal abundance 00307 * dense.xIonDense(nelem,i) is density of ith ionization stage (cm^-3) */ 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 /*end sanity checks */ 00319 00320 /* net loss of line escape and destruction */ 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 /* find energies of all transitions - one line could be a dummy 00326 * also find Boltzmann factors */ 00327 if( a10 == 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( a21 == 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( a20 == 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 /* all lines are ok */ 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 /* check all energies positive */ 00381 ASSERT( ener10 > 0. && ener20 > 0. && ener21 > 0. ); 00382 00383 /* check if energy order is ok */ 00384 ASSERT( ener10 < ener20 && ener21 < ener20 ); 00385 00386 /* check if energy scale is ok */ 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 /* cdsqte is 8.629E-6 / sqrte * eden */ 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 /* 1( DBLE(r01+r02)/DBLE(r10+r01+r02) - DBLE(r01)/DBLE(r10+r12+r01) ) 00415 * beta is factor with n2 */ 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 /* this catches both negative and round off */ 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 /* guard against roundoff - this should really have been zero 00453 * >>chng 96 nov 14, protection against round-off to zero 00454 * >>chng 96 dec 03, made r01, etc, double, and changed limit to 1e-9 */ 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 /* slightly negative solution, probably just round-off, zero it */ 00470 p1 = 0.; 00471 atoms.PopLevels[1] = p1; 00472 nLev3Fail = 3; 00473 } 00474 00475 else 00476 { 00477 /* very negative solution, better punt */ 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 /* population of lowest level */ 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 /* level populations for line opacities */ 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 /* line emission - net emission in line */ 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 /* >>chng 03 oct 04, move to RT_OTS */ 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 /* now add thess lines to ots field, routine works on f not c scale */ 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 /* total intensity used to divide line up - one may be 0 */ 00532 /* >>chng 99 nov 30, rewrite algebra so double prec throughout, 00533 * for very low density models */ 00534 /*TotInten = t21->Emis->xIntensity + t20->xIntensity + t10->xIntensity;*/ 00535 TotInten = t21->Emis->phots * (double)t21->EnergyErg 00536 + t20->Emis->phots * (double)t20->EnergyErg + 00537 t10->Emis->phots * (double)t10->EnergyErg; 00538 00539 /* fraction that was collisionally excited */ 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 /* heating or cooling due to each line */ 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 /* two cases - collisionally excited (usual case) or 00576 * radiatively excited - in which case line can be a heat source 00577 * following are correct heat exchange, they will mix to get correct deriv 00578 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to 00579 * keep stable solution by effectively dividing up heating and cooling, 00580 * so that negative cooling does not occur */ 00581 00582 /* now get net heating or cooling */ 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 /*TotalCooling = p0*(c01*ener10 + c02*ener20) + p1*c12*ener21 - 00591 (p2*(c21*ener21 + c20*ener20) + p1*c10*ener10);*/ 00592 /* >>chng 96 nov 22, very dense cool models, roundoff error 00593 *could cause [OI] 63 mic to be dominant heating, cooling, or 00594 *just zero 00595 * >>chng 96 dec 06, above change caused o iii 1666 cooling to 00596 * be zeroed when important in a model in kirk's grid. was at 1e-6, 00597 * set to 1e-7 00598 * >>chng 96 dec 17, from 1e-7 to 1e-8, caused temp fail */ 00599 /* >>chng 99 nov 29, had been 1e-30 to prevent div by zero (?), 00600 * change to dble max since 1e-30 was very large compared to 00601 * cooling for 1e-10 cm-3 den models */ 00602 /*if( fabs(cnet01/MAX2(1e-30,cool01)) < 1e-8 )*/ 00603 00604 /* >>chng 02 jan 28, min from 1e-8 to 1e-10, conserve.in had massive 00605 * temp failures when no molecules turned on, due to this tripping */ 00606 /*if( fabs(cnet01/MAX2(DBL_MIN,cool01)) < 1e-8 )*/ 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 /* add cooling due to each line, 00655 * heating broken out above, will be added to thermal.heating[0][22] in CoolEvaluate*/ 00656 /* >>chng 99 nov 30, rewrite algebra to keep precision on very low density models*/ 00657 CoolAdd(chLab, t21->WLAng ,cool12); 00658 CoolAdd(chLab, t10->WLAng ,cool01); 00659 CoolAdd(chLab, t20->WLAng ,cool02); 00660 00661 /* derivative of cooling function 00662 * dC/dT = Cooling * ( T(excitation)/T_e^2 - 1/(2T) ) 00663 * in following I assume that a 1-2 exciation will have the 0-2 00664 * exponential in the dcdt term = NOT A TYPO */ 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 /* two t20->t's above are not a typo! 00669 * */ 00670 { 00671 /*@-redef@*/ 00672 enum{DEBUG_LOC=false}; 00673 /*@+redef@*/ 00674 if( DEBUG_LOC ) 00675 { 00676 fprintf(ioQQQ,"atom_level3 nLev3Fail %i\n", 00677 nLev3Fail ); 00678 } 00679 } 00680 return; 00681 }