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 /* HyperfineCreat establish space for hf arrays, reads atomic data from hyperfine.dat */ 00004 /* HyperfineCS - returns collision strengths for hyperfine struc transitions */ 00005 /*H21cm computes rate for H 21 cm from upper to lower excitation by atomic hydrogen */ 00006 /*h21_t_ge_20 compute rate for H 21 cm from upper to lower excitation by atomic hydrogen */ 00007 /*h21_t_lt_20 compute rate for H 21 cm from upper to lower excitation by atomic hydrogen */ 00008 /*H21cm_electron compute H 21 cm rate from upper to lower excitation by electrons - call by CoolEvaluate */ 00009 /*H21cm_H_atom - evaluate H atom spin changing collision rate, called by CoolEvaluate */ 00010 /*H21cm_proton - evaluate proton spin changing H atom collision rate, */ 00011 #include "cddefines.h" 00012 #include "lines_service.h" 00013 #include "phycon.h" 00014 #include "dense.h" 00015 #include "rfield.h" 00016 #include "taulines.h" 00017 #include "iso.h" 00018 #include "trace.h" 00019 #include "hyperfine.h" 00020 #include "physconst.h" 00021 00022 /* H21_cm_pops - fine level populations for 21 cm with Lya pumping included 00023 * called in CoolEvaluate */ 00024 void H21_cm_pops( void ) 00025 { 00026 /*atom_level2( HFLines[0] );*/ 00027 /*return;*/ 00028 /* 00029 things we know on entry to this routine: 00030 total population of 2p: StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop 00031 total population of 1s: StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop 00032 continuum pumping rate (lo-up) inside 21 cm line: HFLines[0].pump 00033 upper to lower collision rate inside 21 cm line: HFLines[0].cs*dense.cdsqte 00034 occupation number inside Lya: OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) 00035 00036 level populations (cm-3) must be computed: 00037 population of upper level of 21cm: HFLines[0].Hi->Pop 00038 population of lower level of 21cm: HFLines[0].Lo->Pop 00039 stimulated emission corrected population of lower level: HFLines[0].Emis->PopOpc 00040 */ 00041 00042 double x, 00043 PopTot; 00044 double a32,a31,a41,a42,a21, occnu_lya , 00045 rate12 , rate21 , pump12 , pump21 , coll12 , coll21, 00046 texc , occnu_lya_23 , occnu_lya_13,occnu_lya_24,occnu_lya_14, texc1, texc2; 00047 a31 = 2.08e8; /* Einstein co-efficient for transition 1p1/2 to 0s1/2 */ 00048 a32 = 4.16e8; /* Einstein co-efficient for transition 1p1/2 to 1s1/2 */ 00049 a41 = 4.16e8; /* Einstein co-efficient for transition 1p3/2 to 0s1/2 */ 00050 a42 = 2.08e8; /* Einstein co-efficient for transition 1p3/2 to 1s1/2 */ 00051 /* These A values are determined from eqn. 17.64 of "The theory of Atomic structure 00052 * and Spectra" by R. D. Cowan 00053 * A hyperfine level has degeneracy Gf=(2F + 1) 00054 * a2p1s = 6.24e8; Einstein co-efficient for transition 2p to 1s */ 00055 a21 = 2.85e-15; /* Einstein co-efficient for transition 1s1/2 to 0s1/2 */ 00056 00057 /* above is spontaneous rate - the net rate is this times escape and destruction 00058 * probabilities */ 00059 a21 *= (HFLines[0].Emis->Pdest + HFLines[0].Emis->Pesc + HFLines[0].Emis->Pelec_esc); 00060 00061 /* >>chng 04 dec 01, add hyperfine.lgLya_pump_21cm, option to turn off Lya pump 00062 * of 21 cm, with no 21cm lya pump */ 00063 occnu_lya = OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) * 00064 hyperfine.lgLya_pump_21cm; 00065 00066 /* >>chng 05 apr 20, GS, the lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/ 00067 texc = TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ); 00068 /* >>chng 05 apr 21, GS, Energy difference between 2p1/2 and 2p3/2 is taken from NSRDS */ 00069 if( texc > 0. ) 00070 { 00071 /* convert to boltz factor, which will applied to occupation number of hiher energy transition */ 00072 texc1 = sexp(0.068/texc); 00073 texc2 = sexp(((82259.272-82258.907)*T1CM)/texc); 00074 } 00075 else 00076 { 00077 texc1 = 0.; 00078 texc2 = 0.; 00079 } 00080 00081 /* the continuum within Lya seen by the two levels is not exactly the same brightness. They 00082 * differ by the exp when Lya is on Wein tail of black body, which must be true if 21 cm is important */ 00083 00084 occnu_lya_23 = occnu_lya; 00085 occnu_lya_13 = occnu_lya*texc1; 00086 occnu_lya_14 = occnu_lya_13*texc2; 00087 occnu_lya_24 = occnu_lya*texc2; 00088 00089 /* this is the 21 cm upward continuum pumping rate [s-1] for the attenuated incident and 00090 * local continuum and including line optical depths */ 00091 pump12 = HFLines[0].Emis->pump; 00092 pump21 = pump12 * HFLines[0].Lo->g / HFLines[0].Hi->g; 00093 00094 /* collision rates s-1 within 1s, 00095 * were multiplied by collider density when evaluated in CoolEvaluate */ 00096 /* ContBoltz is Boltzmann factor for wavelength of line */ 00097 coll12 = HFLines[0].Coll.cs*dense.cdsqte/HFLines[0].Lo->g*rfield.ContBoltz[HFLines[0].ipCont-1]; 00098 coll21 = HFLines[0].Coll.cs*dense.cdsqte/HFLines[0].Hi->g; 00099 00100 /* set up rate (s-1) equations 00101 * all process out of 1 that eventually go to 2 */ 00102 rate12 = 00103 /* collision rate (s-1) from 1 to 2 */ 00104 coll12 + 00105 /* direct external continuum pumping (s-1) in 21 cm line - usually dominated by CMB */ 00106 pump12 + 00107 /* pump rate (s-1) up to 3, times fraction that decay to 2, hence net 1-2 */ 00108 3.*a31*occnu_lya_13 *a32/(a31+a32)+ 00109 /* pump rate (s-1) up to 4, times fraction that decay to 2, hence net 1-2 */ 00110 /* >>chng 05 apr 04, GS, degeneracy corrected from 6 to 3 */ 00111 3.*a41*occnu_lya_14 *a42/(a41+a42); 00112 00113 /* set up rate (s-1) equations 00114 * all process out of 2 that eventually go to 1 */ 00115 /* spontaneous + induced 2 -> 1 by external continuum inside 21 cm line */ 00116 /* >>chng 04 dec 03, do not include spontaneous decay, for numerical stability */ 00117 rate21 = 00118 /* collisional deexcitation */ 00119 coll21 + 00120 /* net spontaneous decay plus external continuum pumping in 21 cm line */ 00121 pump21 + 00122 /* rate from 2 to 3 time fraction that go back to 1, hence net 2 - 1 */ 00123 /* >>chng 05 apr 04,GS, degeneracy corrected from 2 to unity */ 00124 occnu_lya_23*a32 * a31/(a31+a32)+ 00125 occnu_lya_24*a42*a41/(a41+a42); 00126 00127 /* x = HFLines[0].Hi->Pop/HFLines[0].Lo->Pop */ 00128 x = rate12 / SDIV(a21 + rate21); 00129 00130 PopTot = StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]; 00131 /* the Transitions term is the total population of 1s */ 00132 HFLines[0].Hi->Pop = (x/(1.+x))* PopTot; 00133 HFLines[0].Lo->Pop = (1./(1.+x))* PopTot; 00134 00135 /* the population with correction for stimulated emission */ 00136 HFLines[0].Emis->PopOpc = HFLines[0].Lo->Pop*((3*rate21- rate12) + 3*a21)/SDIV(3*(a21+ rate21)); 00137 00138 /* number of escaping line photons, used elsewhere for outward beam */ 00139 HFLines[0].Emis->phots = MAX2(0.,a21*HFLines[0].Hi->Pop); 00140 00141 /* intensity of line */ 00142 HFLines[0].Emis->xIntensity = HFLines[0].Emis->phots*HFLines[0].EnergyErg; 00143 00144 /* ratio of collisional to total (collisional + pumped) excitation */ 00145 HFLines[0].Emis->ColOvTot = coll12 / rate12; 00146 00147 /* finally save the spin temperature */ 00148 if( HFLines[0].Hi->Pop > SMALLFLOAT ) 00149 { 00150 hyperfine.Tspin21cm = TexcLine( &HFLines[0] ); 00151 /* this line must be non-zero - it does strongly mase in limit_compton_hi_t sim - 00152 * in that sim pop ratio goes to unity for a float and TexcLine ret zero */ 00153 if( hyperfine.Tspin21cm == 0. ) 00154 hyperfine.Tspin21cm = phycon.te; 00155 } 00156 else 00157 { 00158 hyperfine.Tspin21cm = phycon.te; 00159 } 00160 00161 return; 00162 } 00163 00164 /*H21cm_electron computes rate for H 21 cm from upper to lower excitation by electrons - call by CoolEvaluate 00165 * >>refer H1 cs Smith, F.J., 1966, Planet. Space Sci 14, 929 */ 00166 double H21cm_electron( double temp ) 00167 { 00168 double hold; 00169 temp = MIN2(1e4 , temp ); 00170 /* following fit is from */ 00171 /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */ 00172 00173 hold = -9.607 + log10( sqrt(temp)) * sexp( pow(log10(temp) , 4.5 ) / 1800. ); 00174 hold = pow(10.,hold ); 00175 return( hold ); 00176 } 00177 00178 /* computes rate for H 21 cm from upper to lower excitation by atomic hydrogen 00179 * from 00180 * >>refer H1 cs Allison, A.C., & Dalgarno A., 1969, ApJ 158, 423 */ 00181 /* the following is the best current survey of 21 cm excitation */ 00182 /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */ 00183 #if 0 00184 STATIC double h21_t_ge_20( double temp ) 00185 { 00186 double y; 00187 double x1, 00188 teorginal = temp; 00189 /* data go up to 1,000K must not go above this */ 00190 temp = MIN2( 1000.,temp ); 00191 x1 =1.0/sqrt(temp); 00192 y =-21.70880995483007-13.76259674006133*x1; 00193 y = exp(y); 00194 00195 /* >>chng 02 feb 14, extrapolate above 1e3 K as per Liszt 2001 recommendation 00196 * page 699 of */ 00197 /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */ 00198 if( teorginal > 1e3 ) 00199 { 00200 y *= pow(teorginal/1e3 , 0.33 ); 00201 } 00202 00203 return( y ); 00204 } 00205 00206 /* this branch for T < 20K, data go down to 1 K */ 00207 STATIC double h21_t_lt_20( double temp ) 00208 { 00209 double y; 00210 double x1; 00211 00212 /* must not go below 1K */ 00213 temp = MAX2( 1., temp ); 00214 x1 =temp*log(temp); 00215 y =9.720710314268267E-08+6.325515312006680E-08*x1; 00216 return(y*y); 00217 } 00218 #endif 00219 00220 /* >> chng 04 dec 15, GS. The fitted rate co-efficients (cm3s-1) in the temperature range 1K to 300K is from 00221 * >>refer H1 cs Zygelman, B. 2005, ApJ preprint doi:10.1086/'427682". 00222 * The rate is 4/3 times the Dalgarno (1969) rate for the 00223 temperature range 300K to 1000K. Above 1000K, the rate is extrapolated according to Liszt 2001.*/ 00224 STATIC double h21_t_ge_10( double temp ) 00225 { 00226 double y; 00227 double x1,x2,x3, 00228 teorginal = temp; 00229 /* data go up to 300K */ 00230 temp = MIN2( 300., temp ); 00231 x1 =temp; 00232 y =1.4341127e-9+9.4161077e-15*x1-9.2998995e-9/(log(x1))+6.9539411e-9/sqrt(x1)+1.7742293e-8*(log(x1))/pow(x1,2); 00233 if( teorginal > 300. ) 00234 { 00235 /* data go up to 1000*/ 00236 x3 = MIN2( 1000., teorginal ); 00237 x2 =1.0/sqrt(x3); 00238 y =-21.70880995483007-13.76259674006133*x2; 00239 y = 1.236686*exp(y); 00240 00241 } 00242 if( teorginal > 1e3 ) 00243 { 00244 /*data go above 1000*/ 00245 y *= pow(teorginal/1e3 , 0.33 ); 00246 } 00247 return( y ); 00248 } 00249 /* this branch for T < 10K, data go down to 1 K */ 00250 STATIC double h21_t_lt_10( double temp ) 00251 { 00252 double y; 00253 double x1; 00254 00255 /* must not go below 1K */ 00256 temp = MAX2(1., temp ); 00257 x1 =temp; 00258 y =8.5622857e-10+2.331358e-11*x1+9.5640586e-11*pow((log(x1)),2)-4.6220869e-10*sqrt(x1)-4.1719545e-10/sqrt(x1); 00259 return(y); 00260 } 00261 00262 /*H21cm_H_atom - evaluate H atom spin changing H atom collision rate, 00263 * called by CoolEvaluate 00264 * >>refer H1 cs Allison, A.C. & Dalgarno, A., 1969, ApJ 158, 423 00265 */ 00266 double H21cm_H_atom( double temp ) 00267 { 00268 double hold; 00269 if( temp >= 10. ) 00270 { 00271 hold = h21_t_ge_10( temp ); 00272 } 00273 else 00274 { 00275 hold = h21_t_lt_10( temp ); 00276 } 00277 00278 return hold; 00279 } 00280 00281 /*H21cm_proton - evaluate proton spin changing H atom collision rate, 00282 * called by CoolEvaluate */ 00283 double H21cm_proton( double temp ) 00284 { 00285 /*>>refer 21cm p coll Furlanetto, S.R. & Furlanetto, M.R. 2007, 00286 *>>refcon MNRAS, doi:10.1111/j.1365-2966.2007.11921.x 00287 * previously had used proton rate, which is 3.2 times H0 rate according to 00288 *>>refer 21cm cs Liszt, H. A&A, 371, 698 */ 00289 /* fit to table 1 of first paper */ 00290 /*--------------------------------------------------------------* 00291 TableCurve Function: c:\storage\litera~1\21cm\atomic~1\p21cm.c Jun 20, 2007 3:37:50 PM 00292 proton coll deex 00293 X= temperature (K) 00294 Y= rate coefficient (1e-9 cm3 s-1) 00295 Eqn# 4419 y=a+bx+cx^2+dx^(0.5)+elnx/x 00296 r2=0.9999445384690351 00297 r2adj=0.9999168077035526 00298 StdErr=5.559328579039901E-12 00299 Fstat=49581.16793656295 00300 a= 9.588389834316704E-11 00301 b= -5.158891920816405E-14 00302 c= 5.895348443553458E-19 00303 d= 2.05304960232429E-11 00304 e= 9.122617940315725E-10 00305 *--------------------------------------------------------------*/ 00306 00307 double hold; 00308 /* only fit this range, did not include T = 1K point which 00309 * causes an inflection */ 00310 temp = MAX2( 2. , temp ); 00311 temp = MIN2( 2e4 , temp ); 00312 00313 /* within range of fitted rate coefficients */ 00314 double x1,x2,x3,x4; 00315 x1 = temp; 00316 x2 = temp*temp; 00317 x3 = sqrt(temp); 00318 x4 = log(temp)/temp; 00319 hold =9.588389834316704E-11 - 5.158891920816405E-14*x1 00320 +5.895348443553458E-19*x2 + 2.053049602324290E-11*x3 00321 +9.122617940315725E-10*x4; 00322 00323 return hold; 00324 } 00325 00326 /* 00327 * HyperfineCreate, HyperfineCS written July 2001 00328 * William Goddard for Gary Ferland 00329 * This code calculates line intensities for known 00330 * hyperfine transitions. 00331 */ 00332 00333 /* two products, the transition structure HFLines, which contains all information for the lines, 00334 * and nHFLines, the number of these lines. 00335 * 00336 * these are in taulines.h 00337 * 00338 * info to create them contained in hyperfine.dat 00339 * 00340 * abundances of nuclei are also in hyperfine.dat, stored in 00341 */ 00342 00343 /* Ion contains twelve varying temperatures, specified above, used for */ 00344 /* calculating collision strengths. */ 00345 typedef struct 00346 { 00347 double strengths[12]; 00348 } Ion; 00349 00350 static Ion *Strengths; 00351 00352 /* HyperfineCreat establish space for hf arrays, reads atomic data from hyperfine.dat */ 00353 void HyperfineCreate(void) 00354 { 00355 FILE *ioDATA; 00356 char chLine[INPUT_LINE_LENGTH]; 00357 bool lgEOL; 00358 /*double c, h, k, N, Ne, q12, q21, upsilon, x;*/ 00359 realnum spin, wavelength; 00360 long int i, j, mass, nelec, ion, nelem; 00361 00362 DEBUG_ENTRY( "HyperfineCreate()" ); 00363 00364 /* list of ion collision strengths for the temperatures listed in table */ 00365 /* HFLines containing all the data in Hyperfine.dat, and transition is */ 00366 /* defined in cddefines.h */ 00367 00368 /*transition *HFLines;*/ 00369 00370 /* get the line data for the hyperfine lines */ 00371 if( trace.lgTrace ) 00372 fprintf( ioQQQ," Hyperfine opening hyperfine.dat:"); 00373 00374 ioDATA = open_data( "hyperfine.dat", "r" ); 00375 00376 /* first line is a version number and does not count */ 00377 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00378 { 00379 fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n"); 00380 cdEXIT(EXIT_FAILURE); 00381 } 00382 /* count how many lines are in the file, ignoring all lines 00383 * starting with '#' */ 00384 nHFLines = 0; 00385 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00386 { 00387 /* we want to count the lines that do not start with # 00388 * since these contain data */ 00389 if( chLine[0] != '#') 00390 ++nHFLines; 00391 } 00392 00393 /* MALLOC the transition HFLines array */ 00394 HFLines = (transition *)MALLOC( (size_t)(nHFLines)*sizeof(transition) ); 00395 /* initialize array to impossible values to make sure eventually done right */ 00396 for( i=0; i< nHFLines; ++i ) 00397 { 00398 TransitionJunk( &HFLines[i] ); 00399 HFLines[i].Hi = AddState2Stack(); 00400 HFLines[i].Lo = AddState2Stack(); 00401 HFLines[i].Emis = AddLine2Stack( true ); 00402 } 00403 00404 Strengths = (Ion *)MALLOC( (size_t)(nHFLines)*sizeof(Ion) ); 00405 hyperfine.HFLabundance = (realnum *)MALLOC( (size_t)(nHFLines)*sizeof(realnum) ); 00406 00407 /* now rewind the file so we can read it a second time*/ 00408 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 ) 00409 { 00410 fprintf( ioQQQ, " Hyperfine could not rewind hyperfine.dat.\n"); 00411 cdEXIT(EXIT_FAILURE); 00412 } 00413 00414 /* check that magic number is ok, read the line */ 00415 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00416 { 00417 fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n"); 00418 cdEXIT(EXIT_FAILURE); 00419 } 00420 00421 /* check that magic number is ok, scan it in */ 00422 i = 1; 00423 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00424 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00425 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00426 00427 { 00428 /* the following is the set of numbers that appear at the start of hyperfine.dat 01 07 12 */ 00429 static int iYR=2, iMN=10, iDY=28; 00430 if( ( nelem != iYR ) || ( nelec != iMN ) || ( ion != iDY ) ) 00431 { 00432 fprintf( ioQQQ, 00433 " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" ); 00434 fprintf( ioQQQ, 00435 " I expected to find the number %i %i %i and got %li %li %li instead.\n" , 00436 iYR, iMN , iDY , 00437 nelem , nelec , ion ); 00438 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00439 cdEXIT(EXIT_FAILURE); 00440 } 00441 } 00442 00443 /* 00444 * scan the string taken from Hyperfine.dat, parsing into 00445 * needed variables. 00446 * nelem is the atomic number. 00447 * IonStg is the ionization stage. Atom = 1, Z+ = 2, Z++ = 3, etc. 00448 * Aul is used to find the einstein A in the function GetGF. 00449 * most of the variables are floats. 00450 */ 00451 00452 /* this will count the number of lines */ 00453 j = 0; 00454 00455 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00456 { 00457 /* only look at lines without '#' in first col */ 00458 /* make sure line in data file is < 140 char */ 00459 if( chLine[0] != '#') 00460 { 00461 double Aul; 00462 ASSERT( j < nHFLines ); 00463 00464 double help[3]; 00465 sscanf(chLine, "%li%i%le%i%le%le%le%le%le%le%le%le%le%le%le%le%le%le%le", 00466 &mass, 00467 &HFLines[j].Hi->nelem, 00468 &help[0], 00469 &HFLines[j].Hi->IonStg, 00470 &help[1], 00471 &help[2], 00472 &Aul, 00473 &Strengths[j].strengths[0], 00474 &Strengths[j].strengths[1], 00475 &Strengths[j].strengths[2], 00476 &Strengths[j].strengths[3], 00477 &Strengths[j].strengths[4], 00478 &Strengths[j].strengths[5], 00479 &Strengths[j].strengths[6], 00480 &Strengths[j].strengths[7], 00481 &Strengths[j].strengths[8], 00482 &Strengths[j].strengths[9], 00483 &Strengths[j].strengths[10], 00484 &Strengths[j].strengths[11]); 00485 spin = (realnum)help[0]; 00486 hyperfine.HFLabundance[j] = (realnum)help[1]; 00487 wavelength = (realnum)help[2]; 00488 HFLines[j].Emis->Aul = (realnum)Aul; 00489 HFLines[j].Hi->g = (realnum) (2*(spin + .5) + 1); 00490 HFLines[j].Lo->g = (realnum) (2*(spin - .5) + 1); 00491 HFLines[j].WLAng = wavelength * 1e8f; 00492 HFLines[j].EnergyWN = (realnum) (1. / wavelength); 00493 HFLines[j].Emis->gf = (realnum)(GetGF(HFLines[j].Emis->Aul, HFLines[j].EnergyWN, HFLines[j].Hi->g)); 00494 /*fprintf(ioQQQ,"HFLinesss1\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00495 j,HFLines[j].opacity , HFLines[j].Emis->gf , HFLines[j].Emis->Aul , HFLines[j].EnergyWN,HFLines[j].Lo->g);*/ 00496 00497 ASSERT(HFLines[j].Emis->gf > 0.); 00498 /* increment the counter */ 00499 j++; 00500 } 00501 00502 } 00503 00504 /* closing the file */ 00505 fclose(ioDATA); 00506 00507 # if 0 00508 /* for debugging and developing only */ 00509 /* calculating the luminosity for each isotope */ 00510 for(i = 0; i < nHFLines; i++) 00511 { 00512 N = dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1]; 00513 Ne = dense.eden; 00514 00515 h = 6.626076e-27; /* erg * sec */ 00516 c = 3e10; /* cm / sec */ 00517 k = 1.380658e-16; /* erg / K */ 00518 00519 upsilon = HyperfineCS(i); 00520 /*statistical weights must still be identified */ 00521 q21 = COLL_CONST * upsilon / (phycon.sqrte * HFLines[i].Hi->g); 00522 00523 q12 = HFLines[i].Hi->g/ HFLines[i].Lo->g * q21 * exp(-1 * h * c * HFLines[i].EnergyWN / (k * phycon.te)); 00524 00525 x = Ne * q12 / (HFLines[i].Emis->Aul * (1 + Ne * q21 / HFLines[i].Aul)); 00526 HFLines[i].xIntensity = N * HFLines[i].Emis->Aul * x / (1.0 + x) * h * c / (HFLines[i].WLAng / 1e8); 00527 00528 } 00529 # endif 00530 00531 return; 00532 } 00533 00534 00535 /*HyperfineCS returns interpolated collision strength for element nelem and ion ion */ 00536 double HyperfineCS( long i ) 00537 { 00538 00539 /*Search HFLines to find i of ion */ 00540 int /*i = 0,*/ j = 0; 00541 # define N_TE_TABLE 12 00542 double slope, upsilon, TemperatureTable[N_TE_TABLE] = {.1e6, .15e6, .25e6, .4e6, .6e6, 00543 1.0e6, 1.5e6, 2.5e6, 4e6, 6e6, 10e6, 15e6}; 00544 00545 DEBUG_ENTRY( "HyperfineCS()" ); 00546 00547 /*while(i < nHFLines+1 && HFLines[i].Hi->nelem != nelem && HFLines[i].Hi->IonStg != ion) 00548 ++i;*/ 00549 ASSERT( i >= 0. && i <= nHFLines ); 00550 00551 /*j = temperature */ 00552 /* calculate actual cooling rate for isotope. */ 00553 /* The temperature-dependent collision strength must first be calculated. */ 00554 /* phycon.te is compared to the first, last and intermediate values of strengths[]*/ 00555 /* and interpolated by taking the log of the collision strength */ 00556 /* and temperature. */ 00557 00558 if( phycon.te <= TemperatureTable[0]) 00559 { 00560 /* temperature below bounds of table */ 00561 upsilon = Strengths[i].strengths[0]; 00562 } 00563 else if( phycon.te >= TemperatureTable[N_TE_TABLE-1]) 00564 { 00565 /* temperature above bounds of table */ 00566 upsilon = Strengths[i].strengths[N_TE_TABLE-1]; 00567 } 00568 else 00569 { 00570 j = 1; 00571 /* want Table[j-1] < te < Table[j] */ 00572 /*while ( phycon.te >= TemperatureTable[j] && j < N_TE_TABLE )*/ 00573 /*while ( TemperatureTable[j] < phycon.te && j < N_TE_TABLE )*/ 00574 while ( j < N_TE_TABLE && TemperatureTable[j] < phycon.te ) 00575 j++; 00576 00577 ASSERT( j >= 0 && j < N_TE_TABLE); 00578 ASSERT( (TemperatureTable[j-1] <= phycon.te ) && (TemperatureTable[j] >= phycon.te) ); 00579 00580 if( fp_equal( phycon.te, TemperatureTable[j] ) ) 00581 { 00582 upsilon = Strengths[i].strengths[j]; 00583 } 00584 /*phycon.te must be less than TemperatureTable[j], greater than TemperatureTable[j-1][0]*/ 00585 else if( phycon.te < TemperatureTable[j]) 00586 { 00587 slope = (log10(Strengths[i].strengths[j-1]) - log10(Strengths[i].strengths[j])) / 00588 (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j])); 00589 00590 upsilon = (log10((phycon.te)) - (log10(TemperatureTable[j-1])-6.))*slope + log10(Strengths[i].strengths[j-1]); 00591 upsilon = pow(10., upsilon); 00592 } 00593 else 00594 { 00595 upsilon = Strengths[i].strengths[j-1]; 00596 } 00597 } 00598 00599 return upsilon; 00600 }