cloudy trunk
Loading...
Searching...
No Matches
hypho.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*hypho - create hydrogenic photoionization cross sections */
4#include "cddefines.h"
5#include "hypho.h"
6/* some numbers used below */
7#define NCM 3000
8#define NFREQ NCM
9/*exp1 routine from ucl group to compute 1-exp(-x)
10 * mod implicit none, and f77 control, gfj '96 jul 11 */
11
12STATIC double exp1(double x)
13{
14 double dx,
15 exp1_v;
16
17 DEBUG_ENTRY( "exp1()" );
18
19 dx = fabs(x);
20 if( dx < 1.0e-9 )
21 {
22 exp1_v = -x;
23 }
24 else if( dx < 1.0e-5 )
25 {
26 exp1_v = ((-x*0.5) - 1.0)*x;
27 }
28 else if( dx < 1.0e-3 )
29 {
30 exp1_v = (((-x*0.1666666667) - 0.5)*x - 1.0)*x;
31 }
32 else
33 {
34 exp1_v = 1.0 - exp(x);
35 }
36
37 return( exp1_v );
38}
39
40/*hypho create hydrogenic photoionization cross sections */
41void hypho(
42 /* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */
43 double zed,
44 /* principal quantum number */
45 long int n,
46 /* lowest angular momentum */
47 long int lmin,
48 /* highest angular momentum */
49 long int lmax,
50 /* scaled lowest energy, in units of zed^2/n^2, at which the cs will be done */
51 double en,
52 /* number of points to do */
53 long int ncell,
54 /* array of energies (in units given above) */
55 realnum anu[],
56 /* calculated photoionization cross section in cm^-2 */
57 realnum bfnu[])
58{
59 long int i,
60 ifp,
61 il,
62 ilmax,
63 j,
64 k,
65 l,
66 lg,
67 ll,
68 llk,
69 lll,
70 llm,
71 lm,
72 lmax1,
73 m,
74 mulp,
75 mulr,
76 muls;
77 /* >>chng 01 sep 11, replace array allocation on stack with
78 * MALLOC to avoid bug in gcc 3.0 on Sparc platforms, PvH */
79/* mm[NFREQ], */
80 long *mm;
81 double a,
82 ai,
83 al,
84 alfac,
85 con2,
86 e,
87 ee,
88 fll,
89 flm,
90 fn,
91 fth,
92 g11,
93 g12,
94 g21,
95 g22,
96 g31,
97 g32,
98 gl,
99 gn0,
100 gn1e,
101 gne,
102 gnt,
103 p,
104 p1,
105 p2,
106 se,
107 sl,
108 sl4,
109 sm,
110 sm4,
111 sn,
112 sn4,
113 sum,
114 zed2;
115
116 static double ab[NFREQ],
117 alo[7000],
118 fal[7000],
119 freq[NFREQ],
120 g2[2][NFREQ],
121 g3[2][NFREQ];
122
123 double anl,
124 con3,
125 fac,
126 x;
127 static int first = true;
128 static double zero = 0.;
129 static double one = 1.;
130 static double two = 2.;
131 static double four = 4.;
132 static double ten = 10.;
133 static double con1 = 0.8559;
134
135 DEBUG_ENTRY( "hypho()" );
136
137 /*generate hydrogenic photoionization cross sections */
138 /* *************************************************************
139 * GENERATE HYDROGENIC PHOTOIONIZATION CROSS SECTIONS
140 *
141 * Hypho calculates Hydrogenic bound-free (BF) xsectns for each n,l .
142 * For Opacity work we calculate l-weighted average for each n .
143 * BF xsectns for Hydorgenic ion are tabulated at E/z**2. E is
144 * the energy in Ry corresponding to the frequencies on the Opacity
145 * mesh. It can changed accoding to need.
146 *
147 * zed=Z-Nelc,
148 * n=principal quantum number
149 * lmin=lowest angular momentum considered of level n
150 * lmx=highest angular momentum considered of level n
151 * en=zed*zed/(n*n) lowest energy at which the hydrogenic cross sections are
152 * to be calculated
153 * anu = photon energy
154 * bfnu=photoionization cross section in cm^-2
155 *
156 * local variables */
157 /* CON1=conversion factor from a0**2 to Mb, x-sectns in megabarns
158 * CON1=8.5594e-19 * 1.e+18 */
159
160 mm = (long*)MALLOC((size_t)(NFREQ*sizeof(long)));
161
162 if( ncell > NCM )
163 {
164 fprintf( stderr, " increase ncm in hypho to at least%5ld\n",
165 ncell );
167 }
168
169 if( first )
170 {
171 fal[0] = zero;
172 for( i=1; i < 7000; i++ )
173 {
174 ai = (double)(i);
175 alo[i-1] = log10(ai);
176 fal[i] = alo[i-1] + fal[i-1];
177 }
178 first = false;
179 }
180
181 /* these may not be redefined, and so will crash */
182 ll = -INT_MAX;
183 lm = -INT_MAX;
184 anl = -DBL_MAX;
185
186 /* CALCULATION OF HYDROGENIC PHOTOIONIZATION CROSS SECTION
187 *
188 * INITIALIZATION FOR N
189 * * con2=(a0**2*1.e+18)*(n/zed)**2, zed=z-nelc
190 * * for hydrogen, the threshold, fth=1/n**2, and for hydrogenic atom,it is
191 * zed**2/n**2. */
192
193 zed2 = zed*zed;
194 fn = (double)(n);
195 sn = fn*fn;
196 sn4 = four*sn;
197 con2 = con1*POW2(fn/ zed);
198 fth = one/sn;
199
200 gn0 = 2.3052328943 - 2.302585093*fal[n+n-1] - fn*0.61370563888 +
201 alo[n-1]*(fn + one)*2.30258093;
202 lmax1 = MIN2(lmax,n-1);
203 ilmax = n - lmin;
204
205 /* * calculate top-up st.wt for given n and lmin=lmax+1 (R-matrix). subtract
206 * the sum from 2n**2. */
207 gl = 0.;
208 if( lmin != 0 )
209 {
210 for( i=0; i < lmin; i++ )
211 {
212 lg = i;
213 gl += 2.*lg + 1;
214 }
215 }
216
217 gnt = 2.*(POW2( (double)n ) - gl);
218
219 /* define photon energies epnu corresponding to hydrogenic atom with charge
220 * z and freq to hydrogen atom; INITIALIZE G'S */
221 for( i=0; i < ncell; i++ )
222 {
223 mm[i] = 0;
224 bfnu[i] = (realnum)zero;
225 freq[i] = anu[i]*zed2;
226 for( j=0; j < 2; j++ )
227 {
228 g2[j][i] = zero;
229 g3[j][i] = zero;
230 }
231 }
232
233 /* gjf Aug 95 change
234 * original code returned total &(*(*$# if freq(1)<en */
235 freq[0] = MAX2(freq[0],en);
236
237 /* calculate cross section: L LOOP */
238
239 for( il=1; il <= ilmax; il++ )
240 {
241 l = n - il;
242 m = 0;
243 al = (double)(l);
244 k = n - l - 2;
245 con3 = con2/(two*al + one);
246
247 /* FREQUENCY LOOP (FREQ UNITS ARE RYD*ZED**2) */
248 for( ifp=0; ifp < ncell; ifp++ )
249 {
250 if( freq[ifp] < fth )
251 {
252 if( l <= lmax1 )
253 anl = 0.;
254 }
255 else
256 {
257 g11 = zero;
258 g12 = zero;
259 /* >>chng 99 jun 24, move m from below to end, change m-1 to m */
260 /*m += 1;*/
261 se = freq[ifp] - fth;
262 if( se < 0. )
263 {
264 fprintf( stderr, " %4ld%12.4e%12.4e\n", ifp,
265 freq[ifp], fth );
266 }
267
268 /* if ( se .lt. -1.e-30) se = 1.e-06 */
269 e = sqrt(se);
270 x = one + sn*se;
271 if( k < 0 )
272 {
273 if( e >= 0.314 )
274 {
275 ee = 6.2831853/e;
276 p = 0.5*log10(exp1(-ee));
277 }
278 else
279 {
280 p = zero;
281 }
282
283 if( e > 1.0e-6 )
284 {
285 a = two*(fn - atan(fn*e)/e);
286 }
287 else
288 {
289 a = zero;
290 }
291
292 ab[m] = (gn0 + a)/2.302585 - p - (fn + two)*
293 log10(x);
294 gne = 0.1;
295 gn1e = x*gne/(fn + fn);
296 g3[1][m] = gne;
297 g3[0][m] = gn1e;
298 g2[1][m] = gne*fn*x*(fn + fn - one);
299 g2[0][m] = gn1e*(fn + fn - one)*(four +
300 (fn - one)*x);
301 }
302
303 g22 = g2[1][m];
304 g32 = g3[1][m];
305 g21 = g2[0][m];
306 g31 = g3[0][m];
307 muls = mm[m];
308
309 if( k < 0 )
310 {
311
312 /* l.eq.n-1
313 * */
314 g11 = g31;
315 if( l == 0 )
316 g11 = zero;
317 g12 = g32;
318 }
319
320 else if( k == 0 )
321 {
322
323 /* l.eq.n-2
324 * */
325 g11 = g21;
326 if( l == 0 )
327 g11 = zero;
328 g12 = g22;
329 }
330
331 else
332 {
333
334 /* l.lt.n-2
335 * */
336 if( k <= 1 )
337 {
338 ll = n - 1;
339 lm = n - 2;
340 }
341 sl = POW2( (double)ll );
342 sl4 = four*sl;
343 fll = (double)(ll);
344 g12 = (sn4 - sl4 + (two*sl - fll)*x)*g22 -
345 sn4*(sn - sl)*(one + POW2(fll + one)*se)*
346 g32;
347
348 if( l != 0 )
349 {
350 sm = POW2( (double)lm );
351 sm4 = four*sm;
352 flm = (double)(lm);
353 g11 = (sn4 - sm4 + (two*sm + flm)*x)*g21 -
354 sn4*(sn - POW2(flm + one))*(one +
355 sm*se)*g31;
356 g31 = g21;
357 g21 = g11;
358 }
359
360 g32 = g22;
361 g22 = g12;
362
363 if( (m+1) == ncell )
364 {
365 ll -= 1;
366 if( l != 0 )
367 lm -= 1;
368 }
369
370 if( g12 >= 1.e20 )
371 {
372 muls += 35;
373 g22 *= 1.e-35;
374 g32 *= 1.e-35;
375 g12 *= 1.e-35;
376
377 if( l != 0 )
378 {
379 g11 *= 1.e-35;
380 g21 *= 1.e-35;
381 g31 *= 1.e-35;
382 }
383
384 }
385 }
386
387 mm[m] = muls;
388 g2[1][m] = g22;
389 g3[1][m] = g32;
390 g2[0][m] = g21;
391 g3[0][m] = g31;
392
393 alfac = fal[n+l] - fal[n-l-1] + two*(al - fn)*
394 alo[n*2-1];
395
396 p1 = one;
397 lll = l + 1;
398 llm = l - 1;
399 mulr = 0;
400
401 if( llm >= 1 )
402 {
403 for( i=1; i <= llm; i++ )
404 {
405 ai = (double)(i);
406 p1 *= one + ai*ai*se;
407 if( p1 >= 1.e20 )
408 {
409 p1 *= 1.e-10;
410 mulr += 10;
411 }
412 }
413 }
414
415 p2 = p1;
416 llk = llm + 1;
417 if( llk < 1 )
418 llk = 1;
419
420 for( i=llk; i <= lll; i++ )
421 {
422 ai = (double)(i);
423 p2 *= one + ai*ai*se;
424 }
425
426 mulp = 0;
427 while( g12 >= one )
428 {
429 mulp -= 10;
430 g12 *= 1.e-10;
431 if( l != 0 )
432 g11 *= 1.e-10;
433 }
434
435 sum = alfac + (realnum)(mulr) + two*(ab[m] + (realnum)(muls-
436 mulp+1));
437
438 fac = zero;
439
440 /* >>chng 96 apr 19, 35 had been 50 */
441 if( fabs(sum) < 35. )
442 fac = pow(ten,sum);
443 if( l != 0 )
444 g11 *= g11*p1;
445 g12 *= g12*p2;
446
447 /* anl=con3*(g11 *al + g12 *(al+one))
448 * *x*fac */
449
450 if( l <= lmax1 )
451 anl = fac*x*con3*(g11*al + g12*(al + 1.));
452 anl *= 2.*(2.*al + 1.);
453
454 bfnu[ifp] += (realnum)(anl*1e-18);
455 /* >>chng 99 jun 24, move m inc to here */
456 ++m;
457 }
458
459 }
460
461 }
462
463 /* bfmin = 1.e-03 * bfnu(1) / gnt */
464 for( i=0; i < ncell; i++ )
465 {
466 bfnu[i] /= (realnum)gnt;
467 }
468
469 free( mm );
470 return;
471}
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
#define NFREQ
Definition hypho.cpp:8
#define NCM
Definition hypho.cpp:7
STATIC double exp1(double x)
Definition hypho.cpp:12
void hypho(double zed, long int n, long int lmin, long int lmax, double en, long int ncell, realnum anu[], realnum bfnu[])
Definition hypho.cpp:41
void zero(void)
Definition zero.cpp:73