cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hypho.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 
12 STATIC 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 */
41 void 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 );
166  cdEXIT(EXIT_FAILURE);
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 }

Generated for cloudy by doxygen 1.8.1.1