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 /*CoolDima compute cooling due to level 2 lines */ 00004 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "taulines.h" 00008 #include "dense.h" 00009 #include "phycon.h" 00010 #include "conv.h" 00011 #include "thermal.h" 00012 #include "opacity.h" 00013 #include "lines_service.h" 00014 #include "rfield.h" 00015 #include "mewecoef.h" 00016 #include "atoms.h" 00017 #include "cooling.h" 00018 00019 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */ 00020 STATIC double ColStrGBar(transition * t , realnum cs1 ); 00021 00022 void CoolDima(void) 00023 { 00024 long int i, 00025 ion, 00026 nelem; 00027 double cs; 00028 double OTSLevel2, 00029 sumots, 00030 CoolLevel2; 00031 static double TeEvalCS = -1.; 00032 static long int nzoneEval=0; 00033 00034 DEBUG_ENTRY( "CoolDima()" ); 00035 00036 /* >>chng 03 nov 29, add this option to skip rest of routine */ 00037 /* no level2 command sets nWindLine to -1 */ 00038 if( nWindLine<0 ) 00039 return; 00040 00041 /* only force evaluation of cooling when significant, or first call 00042 * in this zone, or not into first zone */ 00043 if( nzone != nzoneEval || conv.lgLevel2_Cool_Imp || !nzone ) 00044 { 00045 nzoneEval = nzone; 00046 00047 /* check whether we need to reevaluate the collision strenghts. 00048 * Do so if large change in temperature or any stage of ionization has been lowered. */ 00049 if( (conv.lgSearch || conv.lgIonStageTrimed || 00050 fabs(phycon.te-TeEvalCS)/phycon.te > 0.05) ) 00051 { 00052 for( i=0; i < nWindLine; i++ ) 00053 { 00054 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00055 { 00056 /* only evaluate cs if positive abundance */ 00057 ion = TauLine2[i].Hi->IonStg; 00058 nelem = TauLine2[i].Hi->nelem; 00059 00060 if( dense.xIonDense[nelem-1][ion-1] > 0. ) 00061 { 00062 /* now generate the collision strength */ 00063 cs = ColStrGBar(&TauLine2[i] , cs1_flag_lev2[i] ); 00064 } 00065 else 00066 { 00067 cs = 1.; 00068 } 00069 /* now put the cs into the line array */ 00070 PutCS(cs,&TauLine2[i] ); 00071 } 00072 } 00073 TeEvalCS = phycon.te; 00074 } 00075 00076 /* now always update the cooling since this adds ots flux */ 00077 for( i=0; i < nWindLine; i++ ) 00078 { 00079 /* we need to call this even if nothing present since then sets zero 00080 * ignore all lines with negative atomic number 00081 * >>chng 97 aug 30, get rid of test for non-pos ipLnNelem 00082 * pointer tells atom_level2 to use existing WineEtc labels */ 00083 /* only call non-hydrogenic, non-he-like lines */ 00084 /* >>chng 02 sug 11, do not include he-like in sum */ 00085 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00086 { 00087 atom_level2(&TauLine2[i] ); 00088 } 00089 } 00090 } 00091 else 00092 { 00093 /* now use old cooling and ots flux */ 00094 for( i=0; i < nWindLine; i++ ) 00095 { 00096 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00097 { 00098 /* recall that these should be trivial cooling, ots rates */ 00099 CoolAdd( " ", TauLine2[i].WLAng , TauLine2[i].Coll.cool); 00100 /*>>chng 03 oct 04, move to RT_OTS */ 00101 /*RT_OTS_AddLine( TauLine2[i].ots , TauLine2[i].ipCont);*/ 00102 } 00103 } 00104 } 00105 00106 /* find ratio of level 1 to level 2 ots to see how important level 2 is. 00107 * if level 2 ots not important then will update dest probs only first time 00108 * in each zone. If significant then continuously update */ 00109 OTSLevel2 = 0.; 00110 CoolLevel2 = 0.; 00111 if( nzone > 1 ) 00112 { 00113 for( i=0; i < nWindLine; i++ ) 00114 { 00115 long int ip = TauLine2[i].ipCont-1; 00116 CoolLevel2 += TauLine2[i].Coll.cool; 00117 if( opac.opacity_abs[ip] > SMALLFLOAT ) 00118 { 00119 OTSLevel2 += TauLine2[i].Emis->ots / ( (realnum)opac.opacity_abs[ip]); 00120 ASSERT( fabs(OTSLevel2) < 1e37 ); 00121 } 00122 } 00123 00124 /* now sum all ots rates */ 00125 sumots = 0.; 00126 for( i=0; i < rfield.nflux; i++ ) 00127 { 00128 sumots += rfield.otslin[i]; 00129 } 00130 /* is the ots contribution significant ? */ 00131 if( OTSLevel2/SDIV(sumots) > 1e-4 ) 00132 { 00133 /* true if level 2 lines were contributors to the ots rates, set in dimacool */ 00134 conv.lgLevel2_OTS_Imp = true; 00135 } 00136 else 00137 { 00138 /* true if level 2 lines were contributors to the ots rates, set in dimacool */ 00139 conv.lgLevel2_OTS_Imp = false; 00140 } 00141 00142 /* set flags saying if important */ 00143 if( CoolLevel2/SDIV(thermal.ctot) > 1e-4 ) 00144 { 00145 /* true if level 2 lines were contributors to the cooling, set in dimacool */ 00146 conv.lgLevel2_Cool_Imp = true; 00147 } 00148 else 00149 { 00150 /* true if level 2 lines were contributors to the cooling, set in dimacool */ 00151 conv.lgLevel2_Cool_Imp = false; 00152 } 00153 00154 { 00155 /* debugging code for Lya problems */ 00156 /*@-redef@*/ 00157 enum {DEBUG_LOC=false}; 00158 /*@+redef@*/ 00159 if( DEBUG_LOC ) 00160 { 00161 fprintf(ioQQQ,"%li\t%.2e\t%.2e\n", 00162 nzone, OTSLevel2/SDIV(sumots) , CoolLevel2/SDIV(thermal.ctot) ); 00163 } 00164 } 00165 } 00166 else 00167 { 00168 /* true if level 2 lines were contributors to the ots rates, set in dimacool */ 00169 conv.lgLevel2_OTS_Imp = true; 00170 /* true if level 2 lines were contributors to the cooling, set in dimacool */ 00171 conv.lgLevel2_Cool_Imp = true; 00172 } 00173 return; 00174 } 00175 00176 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */ 00177 STATIC double ColStrGBar(transition * t , realnum cs1 ) 00178 { 00179 long int i, 00180 j; 00181 double ColStrGBar_v, 00182 a, 00183 b, 00184 c, 00185 d, 00186 e1, 00187 gb, 00188 x, 00189 y; 00190 double xx, 00191 yy; 00192 00193 DEBUG_ENTRY( "ColStrGBar()" ); 00194 00195 /* Calculation of the collision strengths of multiplets. 00196 * Neutrals are recalculated from 00197 * >>refer cs gbar Fisher et al. (1996) 00198 * >>refer cs gbar Gaetz & Salpeter (1983, ApJS 52, 155) and 00199 * >>refer cs gbar Mewe (1972, A&A 20, 215) 00200 * fits for ions. */ 00201 00202 /* routine to implement g-bar data taken from 00203 *>>refer cs gbar Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985, 00204 *>>refercon A&AS, 62, 197 */ 00205 00206 /* zero hydrogenic lines since they are done by iso-sequence */ 00207 if( t->Hi->nelem == t->Hi->IonStg ) 00208 { 00209 ColStrGBar_v = 0.0; 00210 return( ColStrGBar_v ); 00211 } 00212 00213 /*was the block data linked in? */ 00214 ASSERT( MeweCoef.g[1][0] != 0.); 00215 00216 /* which type of transition is this? cs1 is the flag */ 00217 00218 /* >>chng 01 may 30 - cs1 < 0 means a forced collision strength */ 00219 if( cs1 < 0. ) 00220 { 00221 ColStrGBar_v = -cs1; 00222 return( ColStrGBar_v ); 00223 } 00224 00225 /* >>chng 99 feb 27, change to assert */ 00226 ASSERT( cs1 >= 0.05 ); 00227 00228 /* excitation energy over kT */ 00229 y = t->EnergyK/phycon.te; 00230 if( cs1 < 1.5 ) 00231 { 00232 xx = -log10(y); 00233 00234 if( cs1 < 0.5 ) 00235 { 00236 yy = (1.398813573838321 + xx*(0.02943050869177121 + xx* 00237 (-0.4439783893114510 + xx*(0.2316073358577902 + xx*(0.001870493481643103 + 00238 xx*(-0.008227246351067403))))))/(1.0 + xx*(-0.6064792600526370 + 00239 xx*(0.1958559534507252 + xx*(-0.02110452007196644 + 00240 xx*(0.01348743933722316 + xx*(-0.0001944731334371711)))))); 00241 } 00242 00243 else 00244 { 00245 yy = (1.359675968512206 + xx*(0.04636500015069853 + xx* 00246 (-0.4491620298246676 + xx*(0.2498199231048967 + xx*(0.005053803073345794 + 00247 xx*(-0.01015647880244268))))))/(1.0 + xx*(-0.5904799485819767 + 00248 xx*(0.1877833737815317 + xx*(-0.01536634911179847 + 00249 xx*(0.01530712091180953 + xx*(-0.0001909176790831023)))))); 00250 } 00251 00252 ColStrGBar_v = pow((double)10,yy)*t->Emis->gf/(t->EnergyWN * WAVNRYD* 13.6); 00253 } 00254 else 00255 { 00256 i = (long int)cs1; 00257 00258 if( i < 26 ) 00259 { 00260 e1 = log(1.0+1.0/y) - 0.4/POW2(y + 1.0); 00261 a = MeweCoef.g[i-1][0]; 00262 b = MeweCoef.g[i-1][1]; 00263 c = MeweCoef.g[i-1][2]; 00264 d = MeweCoef.g[i-1][3]; 00265 x = (double)t->Hi->nelem - 3.0; 00266 00267 if( i == 14 ) 00268 { 00269 a *= 1.0 - 0.5/x; 00270 b = 1.0 - 0.8/x; 00271 c *= 1.0 - 1.0/x; 00272 } 00273 00274 else if( i == 16 ) 00275 { 00276 a *= 1.0 - 0.9/x; 00277 b *= 1.0 - 1.7/x; 00278 c *= 1.0 - 2.1/x; 00279 } 00280 00281 else if( i == 18 ) 00282 { 00283 a *= 1.0 + 2.0/x; 00284 b *= 1.0 - 0.7/x; 00285 } 00286 00287 gb = a + (b*y - c*y*y + d)*e1 + c*y; 00288 00289 /* ipLnRyd is exciation energy in Rydbergs */ 00290 ColStrGBar_v = 14.510395*t->Emis->gf*gb/(t->EnergyWN * WAVNRYD); 00291 /* following i>=26 */ 00292 } 00293 00294 else 00295 { 00296 /* 210 is the dimem of g, so [209] is largest val */ 00297 if( i < 210 ) 00298 { 00299 j = (long)(MeweCoef.g[i-1][3]); 00300 if( j == 1 ) 00301 { 00302 ColStrGBar_v = t->Lo->g*MeweCoef.g[i-1][0]* 00303 pow(phycon.te/pow(10.,(double)MeweCoef.g[i-1][2]),(double)MeweCoef.g[i-1][1]); 00304 } 00305 else 00306 { 00307 ColStrGBar_v = t->Lo->g*MeweCoef.g[i-1][0]* 00308 sexp(MeweCoef.g[i-1][1]*(pow(10.,(double)MeweCoef.g[i-1][2])/ 00309 phycon.te)); 00310 } 00311 } 00312 00313 else 00314 { 00315 /* This is for AlII 1670 line only! 00316 * ColStrGBar=0.0125*te**0.603 */ 00317 /* 98 dec 27, this is still in use */ 00318 ColStrGBar_v = 0.0125*phycon.sqrte*phycon.te10* 00319 phycon.te003; 00320 } 00321 } 00322 } 00323 00324 /* following to make sure that negative values not returned */ 00325 ColStrGBar_v = MAX2(ColStrGBar_v,1e-10); 00326 return( ColStrGBar_v ); 00327 }