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 /*AtomSeqBeryllium compute level populations and emissivity for Be-sequence ions */ 00004 #include "cddefines.h" 00005 #include "phycon.h" 00006 #include "lines_service.h" 00007 #include "thirdparty.h" 00008 #include "dense.h" 00009 #include "cooling.h" 00010 #include "thermal.h" 00011 #include "atoms.h" 00012 00013 void AtomSeqBeryllium(double cs12, 00014 double cs13, 00015 double cs23, 00016 transition * t, 00017 double a30) 00018 { 00019 00020 char chLab[5]; 00021 00022 bool lgNegPop; 00023 00024 int32 ipiv[4], nerror; 00025 long int i, j; 00026 00027 double AbunxIon, 00028 Enr01, 00029 Enr10, 00030 a20, 00031 boltz, 00032 c01, 00033 c02, 00034 c03, 00035 c10, 00036 c12, 00037 c13, 00038 c20, 00039 c21, 00040 c23, 00041 c30, 00042 c31, 00043 c32, 00044 coolng, 00045 cs1u, 00046 excit, 00047 r01, 00048 r02, 00049 r03, 00050 r10, 00051 r12, 00052 r13, 00053 r20, 00054 r21, 00055 r23, 00056 r30, 00057 r31, 00058 r32, 00059 ri02, 00060 ri20, 00061 tot20; 00062 00063 double amat[4][4], 00064 bvec[4], 00065 zz[5][5]; 00066 00067 DEBUG_ENTRY( "AtomSeqBeryllium()" ); 00068 00069 /* function returns total emission in both components of line 00070 * destruction by background opacity computed and added to otslin field 00071 * here, 00072 * total cooling computed and added via CoolAdd 00073 * 00074 * finds population of four level Be-like atom 00075 * 00076 * level 1 = ground, J=0 00077 * levels 2,3,4 are J=0,1,2, OF 3P. 00078 * levels 3-1 is the fast one, j=1 to ground j=0 00079 * 00080 * cs1u is coll str to all J sub-levels of 3P; C.S. goes as 2J+1 00081 * when routine exits the collision strength in the fast line is 1/3 of the entry value, 00082 * so that punch lines data does give the right critical density */ 00083 00084 /* AtomSeqBeryllium(cs12,cs13,cs23,t->t,a30,chLab) 00085 * dense.xIonDense[Hi.nelem,i] is density of ith ionization stage (cm^-3) */ 00086 AbunxIon = dense.xIonDense[t->Hi->nelem -1][t->Hi->IonStg-1]; 00087 00088 /* put null terminated line label into chLab using line array*/ 00089 chIonLbl(chLab, t); 00090 boltz = t->EnergyK/phycon.te; 00091 00092 /* set the cs before if below, since we must reset to line cs in all cases */ 00093 cs1u = t->Coll.cs; 00094 /* >>chng 01 sep 09, reset cs coming in, to propoer cs for the ^3P_1 level only 00095 * so that critical density is printed correctly with punch lines data command */ 00096 t->Coll.cs /= 3.f; 00097 00098 /* low density cutoff to keep matrix happy */ 00099 if( AbunxIon <= 1e-20 || boltz > 30. ) 00100 { 00101 /* put in pop since possible just too cool */ 00102 t->Lo->Pop = AbunxIon; 00103 t->Emis->PopOpc = AbunxIon; 00104 t->Hi->Pop = 0.; 00105 t->Emis->xIntensity = 0.; 00106 t->Coll.cool = 0.; 00107 t->Emis->phots = 0.; 00108 /*>>chng 03 oct 04, move to RT_OTS */ 00109 /*t->ots = 0.;*/ 00110 t->Emis->ColOvTot = 0.; 00111 t->Coll.heat = 0.; 00112 /* level populations */ 00113 atoms.PopLevels[0] = AbunxIon; 00114 atoms.PopLevels[1] = 0.; 00115 atoms.PopLevels[2] = 0.; 00116 atoms.PopLevels[3] = 0.; 00117 atoms.DepLTELevels[0] = 1.; 00118 atoms.DepLTELevels[1] = 0.; 00119 atoms.DepLTELevels[2] = 0.; 00120 atoms.DepLTELevels[3] = 0.; 00121 CoolAdd(chLab, t->WLAng ,0.); 00122 return; 00123 } 00124 excit = exp(-boltz); 00125 00126 /* these must be the statistical weights */ 00127 ASSERT( t->Lo->g == 1. ); 00128 ASSERT( t->Hi->g == 3. ); 00129 00130 /* collision strength must be positive */ 00131 ASSERT( cs1u > 0. ); 00132 00133 /* incuded rates for fastest transition */ 00134 ri02 = t->Emis->pump; 00135 00136 /* back reaction has ratio of stat weights */ 00137 ri20 = t->Emis->pump*1./3.; 00138 00139 /* net rate out of level 3, with destruction */ 00140 a20 = t->Emis->Aul*(t->Emis->Pesc + t->Emis->Pelec_esc + t->Emis->Pdest); 00141 tot20 = a20 + ri20; 00142 00143 /* rates between j=0, lowest 3P level, 00144 * 1/9 is ratio of level stat weight to term stat weight */ 00145 c10 = cs1u*dense.cdsqte/9.; 00146 c01 = c10*excit; 00147 r01 = c01; 00148 r10 = c10; 00149 00150 /* stat weights canceled out here */ 00151 c20 = c10; 00152 c02 = c01*3.; 00153 r02 = c02 + ri02; 00154 r20 = c20 + tot20; 00155 00156 c30 = c10; 00157 c03 = c01*5.; 00158 r30 = c30 + a30; 00159 r03 = c03; 00160 00161 c21 = cs12*dense.cdsqte/3.; 00162 c12 = c21*3.; 00163 r12 = c12; 00164 r21 = c21; 00165 00166 c31 = cs13*dense.cdsqte/5.; 00167 c13 = c31*5.; 00168 r31 = c31; 00169 r13 = c13; 00170 00171 c32 = cs23*dense.cdsqte/5.; 00172 c23 = c32*1.667; 00173 r32 = c32; 00174 r23 = c23; 00175 00176 /* set up matrix */ 00177 for( i=0; i <= 3; i++ ) 00178 { 00179 /* first equation will be sum to abund */ 00180 zz[i][0] = 1.; 00181 zz[4][i] = 0.; 00182 } 00183 00184 /* zz(0,4) = AbunxIon */ 00185 zz[4][0] = 1.; 00186 00187 /* ground level 0 is implicit 00188 * level 1 balance equation */ 00189 zz[0][1] = -r01; 00190 zz[1][1] = r10 + r12 + r13; 00191 zz[2][1] = -r21; 00192 zz[3][1] = -r31; 00193 00194 /* level 2 balance equation */ 00195 zz[0][2] = -r02; 00196 zz[1][2] = -r12; 00197 zz[2][2] = r20 + r21 + r23; 00198 zz[3][2] = -r32; 00199 00200 /* level 3 balance equation */ 00201 zz[0][3] = -r03; 00202 zz[1][3] = -r13; 00203 zz[2][3] = -r23; 00204 zz[3][3] = r30 + r31 + r32; 00205 00206 /* this one may be more robust */ 00207 for( j=0; j <= 3; j++ ) 00208 { 00209 for( i=0; i <= 3; i++ ) 00210 { 00211 amat[i][j] = zz[i][j]; 00212 } 00213 bvec[j] = zz[3+1][j]; 00214 } 00215 00216 nerror = 0; 00217 00218 getrf_wrapper(4, 4, (double*)amat, 4, ipiv, &nerror); 00219 getrs_wrapper('N', 4, 1, (double*)amat, 4, ipiv, bvec, 4, &nerror); 00220 00221 if( nerror != 0 ) 00222 { 00223 fprintf( ioQQQ, " AtomSeqBeryllium: dgetrs finds singular or ill-conditioned matrix\n" ); 00224 cdEXIT(EXIT_FAILURE); 00225 } 00226 /* now put results back into z so rest of code treates only 00227 * one case - as if matin1 had been used */ 00228 for( i=0; i <= 3; i++ ) 00229 { 00230 zz[3+1][i] = bvec[i]; 00231 } 00232 00233 lgNegPop = false; 00234 for( i=0; i <= 3; i++ ) 00235 { 00236 atoms.PopLevels[i] = zz[4][i]*AbunxIon; 00237 if( atoms.PopLevels[i] < 0. ) 00238 lgNegPop = true; 00239 } 00240 00241 /* check for negative level populations, this would be a major error */ 00242 if( lgNegPop ) 00243 { 00244 fprintf( ioQQQ, " AtomSeqBeryllium finds non-positive pop,=" ); 00245 for( i=0; i <= 3; i++ ) 00246 { 00247 fprintf( ioQQQ, "%g ", atoms.PopLevels[i] ); 00248 } 00249 fprintf( ioQQQ, "%s \n", chLab ); 00250 fprintf( ioQQQ, " te=%g total abund=%g boltz=%g \n", 00251 phycon.te, AbunxIon, boltz ); 00252 cdEXIT(EXIT_FAILURE); 00253 } 00254 00255 /* convert level populations over to departure coeficients relative to ground */ 00256 atoms.DepLTELevels[0] = 1.; 00257 atoms.DepLTELevels[1] = (atoms.PopLevels[1]/atoms.PopLevels[0])/ 00258 excit; 00259 atoms.DepLTELevels[2] = (atoms.PopLevels[2]/atoms.PopLevels[0])/ 00260 (excit*3.); 00261 atoms.DepLTELevels[3] = (atoms.PopLevels[3]/atoms.PopLevels[0])/ 00262 (excit*5.); 00263 00264 /* compute ratio Aul/(Aul+Cul) */ 00265 t->Emis->ColOvTot = c02/r02; 00266 00267 t->Lo->Pop = AbunxIon; 00268 00269 /* >>chng 96 sep 12, ipLnPopu had not been set before */ 00270 t->Hi->Pop = atoms.PopLevels[2]; 00271 00272 t->Emis->PopOpc = AbunxIon - atoms.PopLevels[2]*1./3.; 00273 00274 /* this will be escaping part of line 00275 * number of escaping line photons, used elsewhere for outward beam */ 00276 t->Emis->phots = atoms.PopLevels[2] * t->Emis->Aul * ( t->Emis->Pesc + t->Emis->Pelec_esc ); 00277 00278 t->Emis->xIntensity = t->Emis->phots * t->EnergyErg; 00279 00280 /* two cases - collisionally excited (usual case) or 00281 * radiatively excited - in which case line can be a heat source 00282 * following are correct heat exchange, they will mix to get correct deriv 00283 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to 00284 * keep stable solution by effectively dividing up heating and cooling, 00285 * so that negative cooling does not occur */ 00286 Enr01 = atoms.PopLevels[0]*(c01 + c02 + c03)*t->EnergyErg; 00287 Enr10 = (atoms.PopLevels[1]*c10 + atoms.PopLevels[2]*c20 + 00288 atoms.PopLevels[3]*c30)*t->EnergyErg; 00289 00290 /* net cooling due to excit minus part of de-excit */ 00291 t->Coll.cool = Enr01 - Enr10*t->Emis->ColOvTot; 00292 00293 /* net heating is remainder */ 00294 t->Coll.heat = Enr10*(1. - t->Emis->ColOvTot); 00295 00296 /* put line cooling into cooling stack */ 00297 coolng = t->Coll.cool; 00298 CoolAdd(chLab, t->WLAng ,coolng); 00299 00300 /* derivative of cooling function */ 00301 thermal.dCooldT += coolng*(t->EnergyK*thermal.tsq1 - thermal.halfte); 00302 00303 /* >>chhg 03 oct 04, move to RT_OTS */ 00304 /* destroyed part of line 00305 dest = atoms.PopLevels[2]*t->Emis->Aul*t->Pdest; 00306 t->ots = dest; */ 00307 00308 /* now add to ots field 00309 * >>chng 03 oct 03, moved to RT_OTS 00310 RT_OTS_AddLine(dest , t->ipCont );*/ 00311 return; 00312 }