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_oi drive the solution of OI level populations, Ly-beta pumping */ 00004 /*oi_level_pops get OI level population with Ly-beta pumping */ 00005 #include "cddefines.h" 00006 #include "taulines.h" 00007 #include "doppvel.h" 00008 #include "iso.h" 00009 #include "trace.h" 00010 #include "dense.h" 00011 #include "rt.h" 00012 #include "rfield.h" 00013 #include "phycon.h" 00014 #include "lines_service.h" 00015 #include "thirdparty.h" 00016 #include "atoms.h" 00017 00018 /*oi_level_pops get OI level population with Ly-beta pumping */ 00019 STATIC void oi_level_pops(double abundoi, 00020 double *coloi); 00021 00022 /*atom_oi drive the solution of OI level populations, Ly-beta pumping */ 00023 void atom_oi_calc(double *coloi) 00024 { 00025 long int i; 00026 double esin, 00027 eslb, 00028 esoi, 00029 flb, 00030 foi, 00031 opaclb, 00032 opacoi, 00033 tin, 00034 tout, 00035 xlb, 00036 xoi; 00037 static double esab; 00038 static double aoi = TauLines[ipTO1025].Emis->Aul; 00039 static double alb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul; 00040 00041 DEBUG_ENTRY( "atom_oi_calc()" ); 00042 00043 /* A's from Pradhan; OI pump line; Ly beta, 8446 */ 00044 00045 /* Notation; 00046 * PMPH31 net rate hydrogen level 3 depopulated to 1 00047 * PMPO15 and PMPO51 are similar, but for oxygen 00048 * ESCH31 is emission in 1025 transition */ 00049 00050 /* called by LINES before calc really start, protect here 00051 * also used for cases where OI not present */ 00052 if( dense.xIonDense[ipOXYGEN][0] <= 0. ) 00053 { 00054 for( i=0; i < 6; i++ ) 00055 { 00056 atoms.popoi[i] = 0.; 00057 } 00058 /* return zero */ 00059 *coloi = 0.; 00060 atoms.pmpo15 = 0.; 00061 atoms.pmpo51 = 0.; 00062 atoms.pmph31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul* 00063 (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc + 00064 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc); 00065 00066 atoms.esch31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul* 00067 (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc + 00068 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc); 00069 00070 /* all trace output turned on with "trace ly beta command' */ 00071 if( trace.lgTr8446 && trace.lgTrace ) 00072 { 00073 fprintf( ioQQQ, 00074 " P8446 called for first time, finds A*escape prob 3-1 =%10.2e\n", 00075 atoms.esch31 ); 00076 } 00077 return; 00078 } 00079 00080 if( iteration > 1 ) 00081 { 00082 /* two sided escape prob */ 00083 tin = TauLines[ipTO1025].Emis->TauIn; 00084 esin = esc_CRDwing_1side(tin,1e-4); 00085 tout = (TauLines[ipTO1025].Emis->TauTot)*0.9 - 00086 TauLines[ipTO1025].Emis->TauIn; 00087 00088 if( trace.lgTr8446 && trace.lgTrace ) 00089 { 00090 fprintf( ioQQQ, " P8446 tin, tout=%10.2e%10.2e\n", tin, tout ); 00091 } 00092 00093 if( tout > 0. ) 00094 { 00095 tout = TauLines[ipTO1025].Emis->TauTot - TauLines[ipTO1025].Emis->TauIn; 00096 00097 /* do not update esab if we overran optical depth scale */ 00098 esab = 0.5*(esin + esc_CRDwing_1side(tout,1e-4)); 00099 } 00100 } 00101 else 00102 { 00103 /* one-sided escape probability */ 00104 esab = esc_CRDwing_1side(TauLines[ipTO1025].Emis->TauIn,1e-4); 00105 } 00106 00107 esoi =TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc; 00108 eslb =Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + 00109 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc; 00110 00111 /* all trace output turned on with "trace ly beta command' */ 00112 if( trace.lgTr8446 && trace.lgTrace ) 00113 { 00114 fprintf( ioQQQ, 00115 " P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n", 00116 DoppVel.doppler[0], DoppVel.doppler[7], eslb, esoi, esab ); 00117 } 00118 00119 /* find relative opacities for two lines */ 00120 opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/DoppVel.doppler[7]; 00121 opaclb = 1.22e-8*dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/DoppVel.doppler[0]; 00122 00123 /* these are x sub a (OI) and x sub b (ly beta) defined in Elitz+Netz */ 00124 xoi = opacoi/(opacoi + opaclb); 00125 xlb = opaclb/(opacoi + opaclb); 00126 00127 /* find relative line-widths, assume same rest freq */ 00128 foi = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipOXYGEN]; 00129 flb = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipHYDROGEN]* 00130 MAX2(0.,1.- TauLines[ipTO1025].Emis->Pelec_esc - TauLines[ipTO1025].Emis->Pesc); 00131 00132 if( trace.lgTr8446 && trace.lgTrace ) 00133 { 00134 fprintf( ioQQQ, 00135 " P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n", 00136 opaclb, opacoi, xlb, xoi, flb, foi ); 00137 } 00138 00139 /* pumping of OI by L-beta - this goes into OI matrix as 1-5 rate 00140 * lgInducProcess set false with no induced command, usually true */ 00141 if( rfield.lgInducProcess ) 00142 { 00143 atoms.pmpo15 = (realnum)(flb*alb*(StatesElem[ipH_LIKE][ipHYDROGEN][ipH3p].Pop*dense.xIonDense[ipHYDROGEN][1])*xoi*(1. - esab)/ 00144 dense.xIonDense[ipOXYGEN][0]); 00145 /* net decay rate from upper level */ 00146 atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - 00147 esab)*foi)); 00148 } 00149 else 00150 { 00151 atoms.pmpo15 = 0.; 00152 atoms.pmpo51 = 0.; 00153 } 00154 00155 /* find level populations for OI */ 00156 oi_level_pops(dense.xIonDense[ipOXYGEN][0],coloi); 00157 00158 /* escape term from n=3 of H; followed by pump term to n=3 */ 00160 //atoms.pmph31 = alb*(1. - (1. - flb)*(1. - eslb) - xlb*(1. - esab)*flb); 00161 00162 /*pmph13 = xlb*(1. - esab)*foi*aoi*atoms.popoi[4];*/ 00163 atoms.pmph31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul* 00164 (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc); 00165 00166 /* following is actual emission rate, used to predict Ly-beta in lines.for */ 00167 atoms.esch31 = alb*(eslb*(1. - flb) + esab*flb); 00168 00169 /* all trace output turned on with "trace ly beta command' */ 00170 if( trace.lgTr8446 && trace.lgTrace ) 00171 { 00172 fprintf( ioQQQ, " P8446 PMPH31=%10.2e EscP3-1%10.2e ESCH31=%10.2e\n", 00173 atoms.pmph31, atoms.pmph31/alb, atoms.esch31/alb ); 00174 } 00175 00176 /* continuum pumping due to J=1, 0 sub states. 00177 * neglect J=2 since L-Beta very optically thick */ 00178 00179 /* lower level populations */ 00180 TauLines[ipT1304].Lo->Pop = atoms.popoi[0]; 00181 TauLines[ipTO1025].Lo->Pop = atoms.popoi[0]; 00182 TauLines[ipT1039].Lo->Pop = atoms.popoi[0]; 00183 TauLines[ipT8446].Lo->Pop = atoms.popoi[1]; 00184 TauLines[ipT4368].Lo->Pop = atoms.popoi[1]; 00185 TauLines[ipTOI13].Lo->Pop = atoms.popoi[2]; 00186 TauLines[ipTOI11].Lo->Pop = atoms.popoi[2]; 00187 TauLines[ipTOI29].Lo->Pop = atoms.popoi[3]; 00188 TauLines[ipTOI46].Lo->Pop = atoms.popoi[4]; 00189 00190 /* upper level populations */ 00191 TauLines[ipT1304].Hi->Pop = atoms.popoi[1]; 00192 TauLines[ipTO1025].Hi->Pop = atoms.popoi[4]; 00193 TauLines[ipT1039].Hi->Pop = atoms.popoi[3]; 00194 TauLines[ipT8446].Hi->Pop = atoms.popoi[2]; 00195 TauLines[ipT4368].Hi->Pop = atoms.popoi[5]; 00196 TauLines[ipTOI13].Hi->Pop = atoms.popoi[3]; 00197 TauLines[ipTOI11].Hi->Pop = atoms.popoi[4]; 00198 TauLines[ipTOI29].Hi->Pop = atoms.popoi[5]; 00199 TauLines[ipTOI46].Hi->Pop = atoms.popoi[5]; 00200 00201 /* opacity factor 00202 * t1304(ipLnEmis.PopOpc) = (popoi(1)-popoi(2)*3.0) */ 00204 TauLines[ipT1304].Emis->PopOpc = atoms.popoi[0]; 00205 TauLines[ipTO1025].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[4]*0.6); 00206 TauLines[ipT1039].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[3]*3.0); 00207 00208 /* t8446(ipLnEmis.PopOpc) = (popoi(2)-popoi(3)*.33) */ 00210 TauLines[ipT8446].Emis->PopOpc = 00211 (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33)); 00212 TauLines[ipT4368].Emis->PopOpc = (atoms.popoi[1] - atoms.popoi[5]*.33); 00213 TauLines[ipTOI13].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[3]*3.0); 00214 TauLines[ipTOI11].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[4]*0.6); 00215 TauLines[ipTOI29].Emis->PopOpc = (atoms.popoi[3] - atoms.popoi[5]*.33); 00216 TauLines[ipTOI46].Emis->PopOpc = (atoms.popoi[4] - atoms.popoi[5]*1.67); 00217 return; 00218 } 00219 00220 /*oi_level_pops get OI level population with Ly-beta pumping */ 00221 STATIC void oi_level_pops(double abundoi, 00222 double *coloi) 00223 { 00224 bool lgNegPop; 00225 00226 long int i, j; 00227 00228 int32 ipiv[6], ner; 00229 00230 double a21, 00231 a32, 00232 a41, 00233 a43, 00234 a51, 00235 a53, 00236 a62, 00237 a64, 00238 a65, 00239 c12, 00240 c13, 00241 c14, 00242 c15, 00243 c16, 00244 c21, 00245 c23, 00246 c24, 00247 c25, 00248 c26, 00249 c31, 00250 c32, 00251 c34, 00252 c35, 00253 c36, 00254 c41, 00255 c42, 00256 c43, 00257 c45, 00258 c46, 00259 c51, 00260 c52, 00261 c53, 00262 c54, 00263 c56, 00264 c61, 00265 c62, 00266 c63, 00267 c64, 00268 c65, 00269 cs, 00270 deptoi[6], 00271 e12, 00272 e23, 00273 e34, 00274 e45, 00275 e56, 00276 simple; 00277 00278 double amat[6][6], 00279 bvec[6], 00280 zz[7][7]; 00281 00282 static realnum g[6]={9.,3.,9.,3.,15.,9}; 00283 00284 DEBUG_ENTRY( "oilevl()" ); 00285 00286 /* following used for linpac matrix inversion */ 00287 00288 /* compute emission from six level OI atom*/ 00289 00290 /* Boltzmann factors for ContBoltz since collisions not dominant in UV tran 00291 * ipoiex is array lof dEnergy for each level, set in DoPoint */ 00292 e12 = rfield.ContBoltz[atoms.ipoiex[0]-1]; 00293 e23 = rfield.ContBoltz[atoms.ipoiex[1]-1]; 00294 e34 = rfield.ContBoltz[atoms.ipoiex[2]-1]; 00295 e45 = rfield.ContBoltz[atoms.ipoiex[3]-1]; 00296 e56 = rfield.ContBoltz[atoms.ipoiex[4]-1]; 00297 00298 /* total rad rates here have dest by background continuum */ 00299 a21 = TauLines[ipT1304].Emis->Aul*(TauLines[ipT1304].Emis->Pdest+ TauLines[ipT1304].Emis->Pesc + TauLines[ipT1304].Emis->Pelec_esc); 00300 a41 = TauLines[ipT1039].Emis->Aul*(TauLines[ipT1039].Emis->Pdest+ TauLines[ipT1039].Emis->Pesc + TauLines[ipT1039].Emis->Pelec_esc); 00301 a51 = TauLines[ipTO1025].Emis->Aul*(TauLines[ipTO1025].Emis->Pdest+ TauLines[ipTO1025].Emis->Pesc + TauLines[ipTO1025].Emis->Pelec_esc); 00302 a51 = atoms.pmpo51; 00303 a32 = TauLines[ipT8446].Emis->Aul*(TauLines[ipT8446].Emis->Pdest+ TauLines[ipT8446].Emis->Pesc + TauLines[ipT8446].Emis->Pelec_esc); 00304 a62 = TauLines[ipT4368].Emis->Aul*(TauLines[ipT4368].Emis->Pdest+ TauLines[ipT4368].Emis->Pesc + TauLines[ipT4368].Emis->Pelec_esc); 00305 a43 = TauLines[ipTOI13].Emis->Aul*(TauLines[ipTOI13].Emis->Pdest+ TauLines[ipTOI13].Emis->Pesc + TauLines[ipTOI13].Emis->Pelec_esc); 00306 a53 = TauLines[ipTOI11].Emis->Aul*(TauLines[ipTOI11].Emis->Pdest+ TauLines[ipTOI11].Emis->Pesc + TauLines[ipTOI11].Emis->Pelec_esc); 00307 a64 = TauLines[ipTOI29].Emis->Aul*(TauLines[ipTOI29].Emis->Pdest+ TauLines[ipTOI29].Emis->Pesc + TauLines[ipTOI29].Emis->Pelec_esc); 00308 a65 = TauLines[ipTOI46].Emis->Aul*(TauLines[ipTOI46].Emis->Pdest+ TauLines[ipTOI46].Emis->Pesc + TauLines[ipTOI46].Emis->Pelec_esc); 00309 00310 /* even at density of 10^17 excited states not in lte due 00311 * to fast transitions down - just need to kill a21 to get to unity at 10^17*/ 00312 00313 /* the 2-1 transition is 1302, cs Wang and McConkey '92 Jphys B 25, 5461 */ 00314 cs = 2.151e-5*phycon.te/phycon.te03; 00315 PutCS(cs,&TauLines[ipT1304]); 00316 00317 /* the 5-1 transition is 1027, cs Wang and McConkey '92 Jphys B 25, 5461 */ 00318 cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01; 00319 PutCS(cs,&TauLines[ipTO1025]); 00320 c21 = dense.cdsqte*TauLines[ipT1304].Coll.cs/g[1]; 00321 c51 = dense.cdsqte*TauLines[ipTO1025].Coll.cs/g[4]; 00322 00323 /* all following are g-bar approx, g-bar = 0.2 */ 00324 c31 = dense.cdsqte*1.0/g[2]; 00325 PutCS(0.27,&TauLines[ipT1039]); 00326 c41 = dense.cdsqte*TauLines[ipT1039].Coll.cs/g[3]; 00327 c61 = dense.cdsqte*1./g[5]; 00328 00329 c12 = c21*g[1]/g[0]*e12; 00330 c13 = c31*g[2]/g[0]*e12*e23; 00331 c14 = c41*g[3]/g[0]*e12*e23*e34; 00332 c15 = c51*g[4]/g[0]*e12*e23*e34*e45; 00333 c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56; 00334 00335 c32 = dense.cdsqte*85./g[2]; 00336 c42 = dense.cdsqte*85./g[3]; 00337 c52 = dense.cdsqte*85./g[4]; 00338 c62 = dense.cdsqte*85./g[5]; 00339 00340 c23 = c32*g[2]/g[1]*e23; 00341 c24 = c42*g[3]/g[1]*e23*e34; 00342 c25 = c52*g[4]/g[1]*e23*e34*e45; 00343 c26 = c62*g[5]/g[1]*e23*e34*e45*e56; 00344 00345 c43 = dense.cdsqte*70./g[3]; 00346 c53 = dense.cdsqte*312./g[4]; 00347 c63 = dense.cdsqte*1./g[5]; 00348 00349 c34 = c43*g[3]/g[2]*e34; 00350 c35 = c53*g[4]/g[2]*e34*e45; 00351 c36 = c63*g[5]/g[2]*e34*e45*e56; 00352 00353 c54 = dense.cdsqte*50./g[4]; 00354 c64 = dense.cdsqte*415./g[5]; 00355 00356 c45 = c54*g[4]/g[3]*e45; 00357 c46 = c64*g[5]/g[3]*e45*e56; 00358 00359 c65 = dense.cdsqte*400./g[5]; 00360 c56 = c65*g[5]/g[4]*e56; 00361 00364 /* this is check for whether matrix inversion likely to fail */ 00365 simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62); 00366 if( simple < 1e-19 ) 00367 { 00368 atoms.popoi[0] = (realnum)abundoi; 00369 for( i=1; i < 6; i++ ) 00370 { 00371 atoms.popoi[i] = 0.; 00372 } 00373 *coloi = 0.; 00374 return; 00375 } 00376 00377 /*--------------------------------------------------------- */ 00378 00379 for( i=0; i < 6; i++ ) 00380 { 00381 zz[i][0] = 1.0; 00382 zz[6][i] = 0.; 00383 } 00384 00385 /* first equation is sum of populations */ 00386 zz[6][0] = abundoi; 00387 00388 /* level two, 3s 3So */ 00389 zz[0][1] = -c12; 00390 zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21; 00391 zz[2][1] = -c32 - a32; 00392 zz[3][1] = -c42; 00393 zz[4][1] = -c52; 00394 zz[5][1] = -c62 - a62; 00395 00396 /* level three */ 00397 zz[0][2] = -c13; 00398 zz[1][2] = -c23; 00399 zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32; 00400 zz[3][2] = -c43 - a43; 00401 zz[4][2] = -c53 - a53; 00402 zz[5][2] = -c63; 00403 00404 /* level four */ 00405 zz[0][3] = -c14; 00406 zz[1][3] = -c24; 00407 zz[2][3] = -c34; 00408 zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43; 00409 zz[4][3] = -c54; 00410 zz[5][3] = -c64 - a64; 00411 00412 /* level five */ 00413 zz[0][4] = -c15 - atoms.pmpo15; 00414 zz[1][4] = -c25; 00415 zz[2][4] = -c35; 00416 zz[3][4] = -c45; 00417 zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53; 00418 zz[5][4] = -c65 - a65; 00419 00420 /* level six */ 00421 zz[0][5] = -c16; 00422 zz[1][5] = -c26; 00423 zz[2][5] = -c36; 00424 zz[3][5] = -c46; 00425 zz[4][5] = -c56; 00426 zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62; 00427 00428 /* this one may be more robust */ 00429 for( j=0; j < 6; j++ ) 00430 { 00431 for( i=0; i < 6; i++ ) 00432 { 00433 amat[i][j] = zz[i][j]; 00434 } 00435 bvec[j] = zz[6][j]; 00436 } 00437 00438 ner = 0; 00439 00440 getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner); 00441 getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner); 00442 00443 /*DGETRF(6,6,(double*)amat,6,ipiv,&ner);*/ 00444 /*DGETRS('N',6,1,(double*)amat,6,ipiv,bvec,6,&ner);*/ 00445 if( ner != 0 ) 00446 { 00447 fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" ); 00448 cdEXIT(EXIT_FAILURE); 00449 } 00450 00451 /* now put results back into z so rest of code treates only 00452 * one case - as if matin1 had been used */ 00453 for( i=0; i < 6; i++ ) 00454 { 00455 zz[6][i] = bvec[i]; 00456 } 00457 00458 lgNegPop = false; 00459 for( i=0; i < 6; i++ ) 00460 { 00461 atoms.popoi[i] = (realnum)zz[6][i]; 00462 if( atoms.popoi[i] < 0. ) 00463 lgNegPop = true; 00464 } 00465 00466 /* following used to confirm that all dep coef are unity at 00467 * density of 1e17, t=10,000, and all A's set to zero */ 00468 if( trace.lgTrace && trace.lgTr8446 ) 00469 { 00470 deptoi[0] = 1.; 00471 deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]* 00472 e12); 00473 deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]* 00474 e12*e23); 00475 deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]* 00476 e12*e23*e34); 00477 deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]* 00478 e12*e23*e34*e45); 00479 deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]* 00480 e12*e23*e34*e45*e56); 00481 00482 fprintf( ioQQQ, " oilevl finds levl pop" ); 00483 for(i=0; i < 6; i++) 00484 fprintf( ioQQQ, "%11.3e", atoms.popoi[i] ); 00485 fprintf( ioQQQ, "\n" ); 00486 00487 fprintf( ioQQQ, " oilevl finds dep coef" ); 00488 for(i=0; i < 6; i++) 00489 fprintf( ioQQQ, "%11.3e", deptoi[i] ); 00490 fprintf( ioQQQ, "\n" ); 00491 } 00492 00493 /* this happens due to numerical instability in matrix inversion routine */ 00494 if( lgNegPop ) 00495 { 00496 atoms.nNegOI += 1; 00497 00498 fprintf( ioQQQ, " OILEVL finds negative population" ); 00499 for(i=0;i < 6; i++) 00500 fprintf( ioQQQ, "%10.2e", atoms.popoi[i] ); 00501 fprintf( ioQQQ, "\n" ); 00502 00503 fprintf( ioQQQ, " simple 5 =%10.2e\n", simple ); 00504 00505 atoms.popoi[5] = 0.; 00506 atoms.popoi[4] = (realnum)(abundoi*(c15 + atoms.pmpo15)/(a51 + a53 + 00507 c51 + c53)); 00508 atoms.popoi[3] = 0.; 00509 atoms.popoi[2] = (realnum)(atoms.popoi[4]*(a53 + c53)/(a32 + c32)); 00510 atoms.popoi[1] = (realnum)((atoms.popoi[2]*(a32 + c32) + abundoi* 00511 c12)/(a21 + c21)); 00512 00513 atoms.popoi[0] = (realnum)abundoi; 00514 /* write(QQ,'('' OILEVL resets this to simple pop'',1P,6E10.2)') 00515 * 1 popoi */ 00516 } 00517 00518 /* this is total cooling due to model atom, can be neg (heating) */ 00519 *coloi = 00520 (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 + 00521 (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 + 00522 (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 + 00523 (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 + 00524 (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 + 00525 (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 + 00526 (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 + 00527 (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 + 00528 (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13; 00529 00530 return; 00531 }