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 /*hypho - create hydrogenic photoionization cross sections */ 00004 #include "cddefines.h" 00005 #include "hypho.h" 00006 /* some numbers used below */ 00007 #define NCM 3000 00008 #define NFREQ NCM 00009 /*exp1 routine from ucl group to compute 1-exp(-x) 00010 * mod implicit none, and f77 control, gfj '96 jul 11 */ 00011 00012 STATIC double exp1(double x) 00013 { 00014 double dx, 00015 exp1_v; 00016 00017 DEBUG_ENTRY( "exp1()" ); 00018 00019 dx = fabs(x); 00020 if( dx < 1.0e-9 ) 00021 { 00022 exp1_v = -x; 00023 } 00024 else if( dx < 1.0e-5 ) 00025 { 00026 exp1_v = ((-x*0.5) - 1.0)*x; 00027 } 00028 else if( dx < 1.0e-3 ) 00029 { 00030 exp1_v = (((-x*0.1666666667) - 0.5)*x - 1.0)*x; 00031 } 00032 else 00033 { 00034 exp1_v = 1.0 - exp(x); 00035 } 00036 00037 return( exp1_v ); 00038 } 00039 00040 /*hypho create hydrogenic photoionization cross sections */ 00041 void hypho( 00042 /* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */ 00043 double zed, 00044 /* principal quantum number */ 00045 long int n, 00046 /* lowest angular momentum */ 00047 long int lmin, 00048 /* highest angular momentum */ 00049 long int lmax, 00050 /* scaled lowest energy, in units of zed^2/n^2, at which the cs will be done */ 00051 double en, 00052 /* number of points to do */ 00053 long int ncell, 00054 /* array of energies (in units given above) */ 00055 realnum anu[], 00056 /* calculated photoionization cross section in cm^-2 */ 00057 realnum bfnu[]) 00058 { 00059 long int i, 00060 ifp, 00061 il, 00062 ilmax, 00063 j, 00064 k, 00065 l, 00066 lg, 00067 ll, 00068 llk, 00069 lll, 00070 llm, 00071 lm, 00072 lmax1, 00073 m, 00074 mulp, 00075 mulr, 00076 muls; 00077 /* >>chng 01 sep 11, replace array allocation on stack with 00078 * MALLOC to avoid bug in gcc 3.0 on Sparc platforms, PvH */ 00079 /* mm[NFREQ], */ 00080 long *mm; 00081 double a, 00082 ai, 00083 al, 00084 alfac, 00085 con2, 00086 e, 00087 ee, 00088 fll, 00089 flm, 00090 fn, 00091 fth, 00092 g11, 00093 g12, 00094 g21, 00095 g22, 00096 g31, 00097 g32, 00098 gl, 00099 gn0, 00100 gn1e, 00101 gne, 00102 gnt, 00103 p, 00104 p1, 00105 p2, 00106 se, 00107 sl, 00108 sl4, 00109 sm, 00110 sm4, 00111 sn, 00112 sn4, 00113 sum, 00114 zed2; 00115 00116 static double ab[NFREQ], 00117 alo[7000], 00118 fal[7000], 00119 freq[NFREQ], 00120 g2[2][NFREQ], 00121 g3[2][NFREQ]; 00122 00123 double anl, 00124 con3, 00125 fac, 00126 x; 00127 static int first = true; 00128 static double zero = 0.; 00129 static double one = 1.; 00130 static double two = 2.; 00131 static double four = 4.; 00132 static double ten = 10.; 00133 static double con1 = 0.8559; 00134 00135 DEBUG_ENTRY( "hypho()" ); 00136 00137 /*generate hydrogenic photoionization cross sections */ 00138 /* ************************************************************* 00139 * GENERATE HYDROGENIC PHOTOIONIZATION CROSS SECTIONS 00140 * 00141 * Hypho calculates Hydrogenic bound-free (BF) xsectns for each n,l . 00142 * For Opacity work we calculate l-weighted average for each n . 00143 * BF xsectns for Hydorgenic ion are tabulated at E/z**2. E is 00144 * the energy in Ry corresponding to the frequencies on the Opacity 00145 * mesh. It can changed accoding to need. 00146 * 00147 * zed=Z-Nelc, 00148 * n=principal quantum number 00149 * lmin=lowest angular momentum considered of level n 00150 * lmx=highest angular momentum considered of level n 00151 * en=zed*zed/(n*n) lowest energy at which the hydrogenic cross sections are 00152 * to be calculated 00153 * anu = photon energy 00154 * bfnu=photoionization cross section in cm^-2 00155 * 00156 * local variables */ 00157 /* CON1=conversion factor from a0**2 to Mb, x-sectns in megabarns 00158 * CON1=8.5594e-19 * 1.e+18 */ 00159 00160 mm = (long*)MALLOC((size_t)(NFREQ*sizeof(long))); 00161 00162 if( ncell > NCM ) 00163 { 00164 fprintf( stderr, " increase ncm in hypho to at least%5ld\n", 00165 ncell ); 00166 cdEXIT(EXIT_FAILURE); 00167 } 00168 00169 if( first ) 00170 { 00171 fal[0] = zero; 00172 for( i=1; i < 7000; i++ ) 00173 { 00174 ai = (double)(i); 00175 alo[i-1] = log10(ai); 00176 fal[i] = alo[i-1] + fal[i-1]; 00177 } 00178 first = false; 00179 } 00180 00181 /* these may not be redefined, and so will crash */ 00182 ll = -INT_MAX; 00183 lm = -INT_MAX; 00184 anl = -DBL_MAX; 00185 00186 /* CALCULATION OF HYDROGENIC PHOTOIONIZATION CROSS SECTION 00187 * 00188 * INITIALIZATION FOR N 00189 * * con2=(a0**2*1.e+18)*(n/zed)**2, zed=z-nelc 00190 * * for hydrogen, the threshold, fth=1/n**2, and for hydrogenic atom,it is 00191 * zed**2/n**2. */ 00192 00193 zed2 = zed*zed; 00194 fn = (double)(n); 00195 sn = fn*fn; 00196 sn4 = four*sn; 00197 con2 = con1*POW2(fn/ zed); 00198 fth = one/sn; 00199 00200 gn0 = 2.3052328943 - 2.302585093*fal[n+n-1] - fn*0.61370563888 + 00201 alo[n-1]*(fn + one)*2.30258093; 00202 lmax1 = MIN2(lmax,n-1); 00203 ilmax = n - lmin; 00204 00205 /* * calculate top-up st.wt for given n and lmin=lmax+1 (R-matrix). subtract 00206 * the sum from 2n**2. */ 00207 gl = 0.; 00208 if( lmin != 0 ) 00209 { 00210 for( i=0; i < lmin; i++ ) 00211 { 00212 lg = i; 00213 gl += 2.*lg + 1; 00214 } 00215 } 00216 00217 gnt = 2.*(POW2( (double)n ) - gl); 00218 00219 /* define photon energies epnu corresponding to hydrogenic atom with charge 00220 * z and freq to hydrogen atom; INITIALIZE G'S */ 00221 for( i=0; i < ncell; i++ ) 00222 { 00223 mm[i] = 0; 00224 bfnu[i] = (realnum)zero; 00225 freq[i] = anu[i]*zed2; 00226 for( j=0; j < 2; j++ ) 00227 { 00228 g2[j][i] = zero; 00229 g3[j][i] = zero; 00230 } 00231 } 00232 00233 /* gjf Aug 95 change 00234 * original code returned total &(*(*$# if freq(1)<en */ 00235 freq[0] = MAX2(freq[0],en); 00236 00237 /* calculate cross section: L LOOP */ 00238 00239 for( il=1; il <= ilmax; il++ ) 00240 { 00241 l = n - il; 00242 m = 0; 00243 al = (double)(l); 00244 k = n - l - 2; 00245 con3 = con2/(two*al + one); 00246 00247 /* FREQUENCY LOOP (FREQ UNITS ARE RYD*ZED**2) */ 00248 for( ifp=0; ifp < ncell; ifp++ ) 00249 { 00250 if( freq[ifp] < fth ) 00251 { 00252 if( l <= lmax1 ) 00253 anl = 0.; 00254 } 00255 else 00256 { 00257 g11 = zero; 00258 g12 = zero; 00259 /* >>chng 99 jun 24, move m from below to end, change m-1 to m */ 00260 /*m += 1;*/ 00261 se = freq[ifp] - fth; 00262 if( se < 0. ) 00263 { 00264 fprintf( stderr, " %4ld%12.4e%12.4e\n", ifp, 00265 freq[ifp], fth ); 00266 } 00267 00268 /* if ( se .lt. -1.e-30) se = 1.e-06 */ 00269 e = sqrt(se); 00270 x = one + sn*se; 00271 if( k < 0 ) 00272 { 00273 if( e >= 0.314 ) 00274 { 00275 ee = 6.2831853/e; 00276 p = 0.5*log10(exp1(-ee)); 00277 } 00278 else 00279 { 00280 p = zero; 00281 } 00282 00283 if( e > 1.0e-6 ) 00284 { 00285 a = two*(fn - atan(fn*e)/e); 00286 } 00287 else 00288 { 00289 a = zero; 00290 } 00291 00292 ab[m] = (gn0 + a)/2.302585 - p - (fn + two)* 00293 log10(x); 00294 gne = 0.1; 00295 gn1e = x*gne/(fn + fn); 00296 g3[1][m] = gne; 00297 g3[0][m] = gn1e; 00298 g2[1][m] = gne*fn*x*(fn + fn - one); 00299 g2[0][m] = gn1e*(fn + fn - one)*(four + 00300 (fn - one)*x); 00301 } 00302 00303 g22 = g2[1][m]; 00304 g32 = g3[1][m]; 00305 g21 = g2[0][m]; 00306 g31 = g3[0][m]; 00307 muls = mm[m]; 00308 00309 if( k < 0 ) 00310 { 00311 00312 /* l.eq.n-1 00313 * */ 00314 g11 = g31; 00315 if( l == 0 ) 00316 g11 = zero; 00317 g12 = g32; 00318 } 00319 00320 else if( k == 0 ) 00321 { 00322 00323 /* l.eq.n-2 00324 * */ 00325 g11 = g21; 00326 if( l == 0 ) 00327 g11 = zero; 00328 g12 = g22; 00329 } 00330 00331 else 00332 { 00333 00334 /* l.lt.n-2 00335 * */ 00336 if( k <= 1 ) 00337 { 00338 ll = n - 1; 00339 lm = n - 2; 00340 } 00341 sl = POW2( (double)ll ); 00342 sl4 = four*sl; 00343 fll = (double)(ll); 00344 g12 = (sn4 - sl4 + (two*sl - fll)*x)*g22 - 00345 sn4*(sn - sl)*(one + POW2(fll + one)*se)* 00346 g32; 00347 00348 if( l != 0 ) 00349 { 00350 sm = POW2( (double)lm ); 00351 sm4 = four*sm; 00352 flm = (double)(lm); 00353 g11 = (sn4 - sm4 + (two*sm + flm)*x)*g21 - 00354 sn4*(sn - POW2(flm + one))*(one + 00355 sm*se)*g31; 00356 g31 = g21; 00357 g21 = g11; 00358 } 00359 00360 g32 = g22; 00361 g22 = g12; 00362 00363 if( (m+1) == ncell ) 00364 { 00365 ll -= 1; 00366 if( l != 0 ) 00367 lm -= 1; 00368 } 00369 00370 if( g12 >= 1.e20 ) 00371 { 00372 muls += 35; 00373 g22 *= 1.e-35; 00374 g32 *= 1.e-35; 00375 g12 *= 1.e-35; 00376 00377 if( l != 0 ) 00378 { 00379 g11 *= 1.e-35; 00380 g21 *= 1.e-35; 00381 g31 *= 1.e-35; 00382 } 00383 00384 } 00385 } 00386 00387 mm[m] = muls; 00388 g2[1][m] = g22; 00389 g3[1][m] = g32; 00390 g2[0][m] = g21; 00391 g3[0][m] = g31; 00392 00393 alfac = fal[n+l] - fal[n-l-1] + two*(al - fn)* 00394 alo[n*2-1]; 00395 00396 p1 = one; 00397 lll = l + 1; 00398 llm = l - 1; 00399 mulr = 0; 00400 00401 if( llm >= 1 ) 00402 { 00403 for( i=1; i <= llm; i++ ) 00404 { 00405 ai = (double)(i); 00406 p1 *= one + ai*ai*se; 00407 if( p1 >= 1.e20 ) 00408 { 00409 p1 *= 1.e-10; 00410 mulr += 10; 00411 } 00412 } 00413 } 00414 00415 p2 = p1; 00416 llk = llm + 1; 00417 if( llk < 1 ) 00418 llk = 1; 00419 00420 for( i=llk; i <= lll; i++ ) 00421 { 00422 ai = (double)(i); 00423 p2 *= one + ai*ai*se; 00424 } 00425 00426 mulp = 0; 00427 while( g12 >= one ) 00428 { 00429 mulp -= 10; 00430 g12 *= 1.e-10; 00431 if( l != 0 ) 00432 g11 *= 1.e-10; 00433 } 00434 00435 sum = alfac + (realnum)(mulr) + two*(ab[m] + (realnum)(muls- 00436 mulp+1)); 00437 00438 fac = zero; 00439 00440 /* >>chng 96 apr 19, 35 had been 50 */ 00441 if( fabs(sum) < 35. ) 00442 fac = pow(ten,sum); 00443 if( l != 0 ) 00444 g11 *= g11*p1; 00445 g12 *= g12*p2; 00446 00447 /* anl=con3*(g11 *al + g12 *(al+one)) 00448 * *x*fac */ 00449 00450 if( l <= lmax1 ) 00451 anl = fac*x*con3*(g11*al + g12*(al + 1.)); 00452 anl *= 2.*(2.*al + 1.); 00453 00454 bfnu[ifp] += (realnum)(anl*1e-18); 00455 /* >>chng 99 jun 24, move m inc to here */ 00456 ++m; 00457 } 00458 00459 } 00460 00461 } 00462 00463 /* bfmin = 1.e-03 * bfnu(1) / gnt */ 00464 for( i=0; i < ncell; i++ ) 00465 { 00466 bfnu[i] /= (realnum)gnt; 00467 } 00468 00469 free( mm ); 00470 return; 00471 }