cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_adfa.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 /*phfit derive photoionization cross sections for first 30 elements */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "atmdat.h"
7 #include "iso.h"
8 
11 {
12  DEBUG_ENTRY( "t_ADfA()" );
13 
14  /* this option, use the new atmdat_rad_rec recombination rates */
16 
17  double help[9];
18  const long VERSION_MAGIC = 20061204L;
19 
20  const static char chFile[] = "phfit.dat";
21 
22  FILE *io = open_data( chFile, "r" );
23 
24  bool lgErr = false;
25  long i=-1, j=-1, k=-1, n;
26 
27  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
28  if( lgErr || i != VERSION_MAGIC )
29  {
30  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
31  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
32  cdEXIT(EXIT_FAILURE);
33  }
34 
35  for( n=0; n < 7; n++ )
36  lgErr = lgErr || ( fscanf( io, "%ld", &L[n] ) != 1 );
37  for( n=0; n < 30; n++ )
38  lgErr = lgErr || ( fscanf( io, "%ld", &NINN[n] ) != 1 );
39  for( n=0; n < 30; n++ )
40  lgErr = lgErr || ( fscanf( io, "%ld", &NTOT[n] ) != 1 );
41  while( true )
42  {
43  lgErr = lgErr || ( fscanf( io, "%ld %ld %ld", &i, &j, &k ) != 3 );
44  if( i == -1 && j == -1 && k == -1 )
45  break;
46  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le", &help[0], &help[1],
47  &help[2], &help[3], &help[4], &help[5] ) != 6 );
48  for( int l=0; l < 6; ++l )
49  PH1[i][j][k][l] = (realnum)help[l];
50  }
51  while( true )
52  {
53  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
54  if( i == -1 && j == -1 )
55  break;
56  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le", &help[0], &help[1],
57  &help[2], &help[3], &help[4], &help[5], &help[6] ) != 7 );
58  for( int l=0; l < 7; ++l )
59  PH2[i][j][l] = (realnum)help[l];
60  }
61  fclose( io );
62 
63  ASSERT( !lgErr );
64 
65  const static char chFile2[] = "hpfit.dat";
66 
67  io = open_data( chFile2, "r" );
68 
69  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
70  if( lgErr || i != VERSION_MAGIC )
71  {
72  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile2, i );
73  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
74  cdEXIT(EXIT_FAILURE);
75  }
76 
77  for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
78  {
79  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &help[0], &help[1],
80  &help[2], &help[3], &help[4] ) != 5 );
81  for( int l=0; l < 5; ++l )
82  PHH[i][l] = (realnum)help[l];
83  }
84 
85  fclose( io );
86 
87  ASSERT( !lgErr );
88 
89  const static char chFile3[] = "rec_lines.dat";
90 
91  io = open_data( chFile3, "r" );
92 
93  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
94  if( lgErr || i != VERSION_MAGIC )
95  {
96  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile3, i );
97  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
98  cdEXIT(EXIT_FAILURE);
99  }
100 
101  for( i=0; i < 110; i++ )
102  {
103  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
104  &help[3], &help[4], &help[5], &help[6], &help[7] ) != 8 );
105  for( int l=0; l < 8; ++l )
106  P[l][i] = (realnum)help[l];
107  }
108 
109 
110  for( i=0; i < 405; i++ )
111  {
112  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
113  &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
114  for( int l=0; l < 9; ++l )
115  ST[l][i] = (realnum)help[l];
116  }
117 
118  fclose( io );
119 
120  ASSERT( !lgErr );
121 
122  const static char chFile4[] = "rad_rec.dat";
123 
124  io = open_data( chFile4, "r" );
125 
126  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
127  if( lgErr || i != VERSION_MAGIC )
128  {
129  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile4, i );
130  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
131  cdEXIT(EXIT_FAILURE);
132  }
133 
134  while( true )
135  {
136  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
137  if( i == -1 && j == -1 )
138  break;
139  lgErr = lgErr || ( fscanf( io, "%le %le", &help[0], &help[1] ) != 2 );
140  for( int l=0; l < 2; ++l )
141  rrec[i][j][l] = (realnum)help[l];
142  }
143  while( true )
144  {
145  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
146  if( i == -1 && j == -1 )
147  break;
148  lgErr = lgErr || ( fscanf( io, "%le %le %le %le", &help[0], &help[1],
149  &help[2], &help[3] ) != 4 );
150  for( int l=0; l < 4; ++l )
151  rnew[i][j][l] = (realnum)help[l];
152  }
153  while( true )
154  {
155  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
156  if( i == -1 )
157  break;
158  lgErr = lgErr || ( fscanf( io, "%le %le %le", &help[0], &help[1], &help[2] ) != 3 );
159  for( int l=0; l < 3; ++l )
160  fe[i][l] = (realnum)help[l];
161  }
162 
163  fclose( io );
164 
165  ASSERT( !lgErr );
166 
167  const static char chFile5[] = "h_rad_rec.dat";
168 
169  io = open_data( chFile5, "r" );
170 
171  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
172  if( lgErr || i != VERSION_MAGIC )
173  {
174  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile5, i );
175  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
176  cdEXIT(EXIT_FAILURE);
177  }
178 
179  for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
180  {
181  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
182  &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
183  for( int l=0; l < 9; ++l )
184  HRF[i][l] = (realnum)help[l];
185  }
186 
187  fclose( io );
188 
189  ASSERT( !lgErr );
190 
191  const static char chFile6[] = "h_phot_cs.dat";
192 
193  io = open_data( chFile6, "r" );
194 
195  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
196  if( lgErr || i != VERSION_MAGIC )
197  {
198  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile6, i );
199  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
200  cdEXIT(EXIT_FAILURE);
201  }
202 
203  for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
204  {
205  lgErr = lgErr || ( fscanf( io, "%le", &help[0] ) != 1 );
206  STH[i] = (realnum)help[0];
207  }
208 
209  fclose( io );
210 
211  ASSERT( !lgErr );
212 
213  const static char chFile7[] = "coll_ion.dat";
214 
215  io = open_data( chFile7, "r" );
216 
217  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
218  if( lgErr || i != VERSION_MAGIC )
219  {
220  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile7, i );
221  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
222  cdEXIT(EXIT_FAILURE);
223  }
224 
225  while( true )
226  {
227  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
228  if( i == -1 && j == -1 )
229  break;
230  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &CF[i][j][0], &CF[i][j][1],
231  &CF[i][j][2], &CF[i][j][3], &CF[i][j][4] ) != 5 );
232  }
233 
234  fclose( io );
235 
236  ASSERT( !lgErr );
237 
238  const static char chFile8[] = "h_coll_str.dat";
239 
240  io = open_data( chFile8, "r" );
241 
242  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
243  if( lgErr || i != VERSION_MAGIC )
244  {
245  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile8, i );
246  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
247  cdEXIT(EXIT_FAILURE);
248  }
249 
250  while( true )
251  {
252  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
253  if( i == -1 && j == -1 )
254  break;
255  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &HCS[i-2][j-1][0], &HCS[i-2][j-1][1],
256  &HCS[i-2][j-1][2], &HCS[i-2][j-1][3], &HCS[i-2][j-1][4], &HCS[i-2][j-1][5],
257  &HCS[i-2][j-1][6], &HCS[i-2][j-1][7] ) != 8 );
258  }
259 
260  fclose( io );
261 
262  ASSERT( !lgErr );
263 }
264 
265 double t_ADfA::phfit(long int nz,
266  long int ne,
267  long int is,
268  double e)
269 {
270  long int nint,
271  nout;
272  double a,
273  b,
274  crs,
275  einn,
276  p1,
277  q,
278  x,
279  y,
280  z;
281 
282  DEBUG_ENTRY( "phfit()" );
283 
284  /*** Version 3. October 8, 1996.
285  *** Written by D. A. Verner, verner@pa.uky.edu
286  *** Inner-shell ionization energies of some low-ionized species are slightly
287  *** improved to fit smoothly the experimental inner-shell ionization energies
288  *** of neutral atoms.
289  ******************************************************************************
290  *** This subroutine calculates partial photoionization cross sections
291  *** for all ionization stages of all atoms from H to Zn (Z=30) by use of
292  *** the following fit parameters:
293  *** Outer shells of the Opacity Project (OP) elements:
294  *** >>refer all photocs Verner, Ferland, Korista, Yakovlev, 1996, ApJ, in press.
295  *** Inner shells of all elements, and outer shells of the non-OP elements:
296  *** Verner and Yakovlev, 1995, A&AS, 109, 125
297  *** Input parameters: nz - atomic number from 1 to 30 (integer)
298  *** ne - number of electrons from 1 to iz (integer)
299  *** is - shell number (integer)
300  *** e - photon energy, eV
301  *** version - enum, PHFIT96 (default): calculates
302  *** new cross sections, PHFIT95: calculates
303  *** only old Hartree-Slater cross sections
304  *** Output parameter: photoionization cross section, Mb
305  *** Shell numbers:
306  *** 1 - 1s, 2 - 2s, 3 - 2p, 4 - 3s, 5 - 3p, 6 - 3d, 7 - 4s.
307  *** If a species in the ground state has no electrons on the given shell,
308  *** the subroutine returns 0.
309  ****************************************************************************** */
310 
311  crs = 0.0;
312  if( nz < 1 || nz > 30 )
313  {
314  return(crs);
315  }
316 
317  if( ne < 1 || ne > nz )
318  {
319  return(crs);
320  }
321 
322  nout = NTOT[ne-1];
323  if( nz == ne && nz > 18 )
324  nout = 7;
325  if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) ||
326  nz == 25) || nz == 26) )
327  nout = 7;
328  if( is > nout )
329  {
330  return(crs);
331  }
332 
333  if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 )
334  {
335  return(crs);
336  }
337 
338  ASSERT( is >= 1 && is <= 7 );
339 
340  if( e < PH1[is-1][ne-1][nz-1][0] )
341  {
342  return(crs);
343  }
344 
345  nint = NINN[ne-1];
346  if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) )
347  {
348  einn = 0.0;
349  }
350  else
351  {
352  if( ne < 3 )
353  {
354  einn = 1.0e30;
355  }
356  else
357  {
358  einn = PH1[nint-1][ne-1][nz-1][0];
359  }
360  }
361 
362  if( (is <= nint || e >= einn) || version == PHFIT95 )
363  {
364  p1 = -PH1[is-1][ne-1][nz-1][4];
365  y = e/PH1[is-1][ne-1][nz-1][1];
366  q = -0.5*p1 - L[is-1] - 5.5;
367  a = PH1[is-1][ne-1][nz-1][2]*(POW2(y - 1.0) +
368  POW2(PH1[is-1][ne-1][nz-1][5]));
369  b = sqrt(y/PH1[is-1][ne-1][nz-1][3]) + 1.0;
370  crs = a*pow(y,q)*pow(b,p1);
371  }
372  else
373  {
374  if( (is < nout && is > nint) && e < einn )
375  {
376  return(crs);
377  }
378  p1 = -PH2[ne-1][nz-1][3];
379  q = -0.5*p1 - 5.5;
380  x = e/PH2[ne-1][nz-1][0] - PH2[ne-1][nz-1][5];
381  z = sqrt(x*x+POW2(PH2[ne-1][nz-1][6]));
382  a = PH2[ne-1][nz-1][1]*(POW2(x - 1.0) +
383  POW2(PH2[ne-1][nz-1][4]));
384  b = 1.0 + sqrt(z/PH2[ne-1][nz-1][2]);
385  crs = a*pow(z,q)*pow(b,p1);
386  }
387  return(crs);
388 }
389 
390 /* same as hpfit, but energy is relative to threshold */
391 double t_ADfA::hpfit_rel(long int iz,
392  long int n,
393  double e)
394 {
395  long m;
396  double crs , ex , eth;
397 
398  if( n == 0 )
399  {
400  m = 1;
401  }
402  else
403  {
404  if( n == 1 )
405  {
406  m = 2;
407  }
408  else
409  {
410  m = n;
411  }
412  }
413 
414  eth = ph1(0,0,iz-1,0)/POW2((double)m);
415  ex = MAX2(1. , e/eth );
416 
417  crs = hpfit( iz , n , ex );
418  ASSERT( crs > 0. );
419 
420  return crs;
421 }
422 
423 double t_ADfA::hpfit(long int iz,
424  long int n,
425  double e)
426 {
427  long int l,
428  m;
429  double cs,
430  eth,
431  ex,
432  q,
433  x;
434 
435  DEBUG_ENTRY( "hpfit()" );
436 
437  /*state specific photoionization cross sections for model hydrogen atom
438  * Version 1, September 23, 1997
439  ******************************************************************************
440  *** This subroutine calculates state-specific photoionization cross sections
441  *** for hydrogen and hydrogen-like ions.
442  *** Input parameters: iz - atomic number
443  *** n - shell number, from 0 to 400:
444  *** 0 - 1s
445  *** 1 - 2s
446  *** 2 - 2p
447  *** 3 - 3
448  *** ......
449  *** e - photon energy, eV
450  *** return value - cross section, cm^(-2)
451  *******************************************************************************/
452 
453  cs = 0.0;
454 
455  ASSERT( iz > 0 && e>0. );
456 
457  if( n >= NHYDRO_MAX_LEVEL )
458  {
459  fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n );
460  cdEXIT(EXIT_FAILURE);
461  }
462 
463  l = 0;
464  if( n == 2 )
465  {
466  l = 1;
467  }
468  q = 3.5 + l - 0.5*PHH[n][1];
469 
470  if( n == 0 )
471  {
472  m = 1;
473  }
474  else
475  {
476  if( n == 1 )
477  {
478  m = 2;
479  }
480  else
481  {
482  m = n;
483  }
484  }
485 
486  eth = ph1(0,0,iz-1,0)/POW2((double)m);
487  ex = MAX2(1. , e/eth );
488 
489  /* Don't just force to be at least one...make sure e/eth is close to one or greater. */
490  ASSERT( e/eth > 0.95 );
491 
492  if( ex < 1.0 )
493  {
494  return(0.);
495  }
496 
497  x = ex/PHH[n][0];
498  cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/
499  pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/
500  POW2((double)iz));
501  return(cs);
502 }
503 
504 void t_ADfA::rec_lines(double t,
505  realnum r[][471])
506 {
507  long int i,
508  j,
509  ipj1,
510  ipj2;
511 
512  double a,
513  a1,
514  dr[4][405],
515  p1,
516  p2,
517  p3,
518  p4,
519  p5,
520  p6,
521  rr[4][110],
522  te,
523  x,
524  z;
525 
526  static long jd[6]={143,145,157,360,376,379};
527 
528  static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50,
529  52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100,
530  101,103,104};
531 
532  static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112,
533  120,114,143,145,157,152,169,183,200,163,225,223,237,232,235,
534  249,247,300,276,278,376,360,379,384};
535 
536  DEBUG_ENTRY( "rec_lines()" );
537 
538  /*effective recombination coefficients for lines of C, N, O, by D. Verner
539  * Version 2, April 30, 1997
540  ******************************************************************************
541  *** This subroutine calculates effective recombination coefficients
542  *** for 110 permitted recombination lines of C, N, O (Pequignot, Petitjean,
543  *** & Boisson, 1991, A&A, 251, 680) and 405 permitted dielectronic
544  *** recombination lines (Nussbaumer & Storey, 1984, A&AS, 56, 293)
545  *** Input parameter: t - temperature, K
546  *** Output parameters: r(i,j), i=1,471
547  *** r(i,1) - atomic number
548  *** r(i,2) - number of electrons
549  *** r(i,3) - wavelength, angstrom
550  *** r(i,4) - rate coefficient, cm^3 s^(-1)
551  ****************************************************************************** */
552 
553  for( i=0; i < 110; i++ )
554  {
555  rr[0][i] = P[0][i];
556  rr[1][i] = P[1][i];
557  rr[2][i] = P[2][i];
558  z = P[0][i] - P[1][i] + 1.0;
559  te = 1.0e-04*t/z/z;
560  p1 = P[3][i];
561  p2 = P[4][i];
562  p3 = P[5][i];
563  p4 = P[6][i];
564  if( te < 0.004 )
565  {
566  a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
567  a = a1/sqrt(te/0.004);
568  }
569  else
570  {
571  if( te > 2.0 )
572  {
573  a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
574  a = a1/pow(te/2.0,1.5);
575  }
576  else
577  {
578  a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
579  }
580  }
581  rr[3][i] = 1.0e-13*z*a*P[7][i];
582  }
583 
584  for( i=0; i < 405; i++ )
585  {
586  dr[0][i] = ST[0][i];
587  dr[1][i] = ST[1][i];
588  dr[2][i] = ST[2][i];
589  te = 1.0e-04*t;
590  p1 = ST[3][i];
591  p2 = ST[4][i];
592  p3 = ST[5][i];
593  p4 = ST[6][i];
594  p5 = ST[7][i];
595  p6 = ST[8][i];
596  if( te < p6 )
597  {
598  x = p5*(1.0/te - 1.0/p6);
599  if( x > 80.0 )
600  {
601  a = 0.0;
602  }
603  else
604  {
605  a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/
606  p6);
607  a = a1/exp(x);
608  }
609  }
610  else
611  {
612  if( te > 6.0 )
613  {
614  a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/
615  exp(p5/6.0);
616  a = a1/pow(te/6.0,1.5);
617  }
618  else
619  {
620  a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/
621  te);
622  }
623  }
624  dr[3][i] = 1.0e-12*a;
625  }
626 
627  for( i=0; i < 6; i++ )
628  {
629  ipj1 = jd[i];
630  ipj2 = ipj1 + 1;
631  dr[3][ipj1-1] += dr[3][ipj2-1];
632  dr[0][ipj2-1] = 0.0;
633  }
634 
635  for( i=0; i < 38; i++ )
636  {
637  ipj1 = ip[i];
638  ipj2 = id[i];
639  rr[3][ipj1-1] += dr[3][ipj2-1];
640  dr[0][ipj2-1] = 0.0;
641  }
642 
643  for( i=0; i < 110; i++ )
644  {
645  r[0][i] = (realnum)rr[0][i];
646  r[1][i] = (realnum)rr[1][i];
647  r[2][i] = (realnum)rr[2][i];
648  r[3][i] = (realnum)rr[3][i];
649  }
650 
651  j = 110;
652  for( i=0; i < 405; i++ )
653  {
654  if( dr[0][i] > 1.0 )
655  {
656  j += 1;
657  r[0][j-1] = (realnum)dr[0][i];
658  r[1][j-1] = (realnum)dr[1][i];
659  r[2][j-1] = (realnum)dr[2][i];
660  r[3][j-1] = (realnum)dr[3][i];
661  }
662  }
663  return;
664 }
665 
666 double t_ADfA::rad_rec(long int iz,
667  long int in,
668  double t)
669 {
670  /*
671  *** Version 4. June 29, 1999.
672  *** Written by D. A. Verner, verner@pa.uky.edu
673  ******************************************************************************
674  *** This subroutine calculates rates of radiative recombination for all ions
675  *** of all elements from H through Zn by use of the following fits:
676  *** H-like, He-like, Li-like, Na-like -
677  *** >>refer all reccoef Verner & Ferland, 1996, ApJS, 103, 467
678  *** Other ions of C, N, O, Ne - Pequignot et al. 1991, A&A, 251, 680,
679  *** refitted by Verner & Ferland formula to ensure correct asymptotes
680  *** Fe XVII-XXIII -
681  *** >>refer Fe17-23 recom Arnaud & Raymond, 1992, ApJ, 398, 394
682  *** Fe I-XV - refitted by Verner & Ferland formula to ensure correct asymptotes
683  *** Other ions of Mg, Si, S, Ar, Ca, Fe, Ni -
684  *** -
685  *** >>refer all recom Shull & Van Steenberg, 1982, ApJS, 48, 95
686  *** Other ions of Na, Al -
687  *** >>refer Na, Al recom Landini & Monsignori Fossi, 1990, A&AS, 82, 229
688  *** Other ions of F, P, Cl, K, Ti, Cr, Mn, Co (excluding Ti I-II, Cr I-IV,
689  *** Mn I-V, Co I) -
690  *** >>refer many recom Landini & Monsignori Fossi, 1991, A&AS, 91, 183
691  *** All other species - interpolations of the power-law fits
692  *** Input parameters: iz - atomic number
693  *** in - number of electrons from 1 to iz
694  *** t - temperature, K
695  *** return result: - rate coefficient, cm^3 s^(-1)
696  ******************************************************************************
697  */
698  double tt;
699  double rate;
700 
701  DEBUG_ENTRY( "rad_rec()" );
702 
703  rate = 0.0;
704  if( iz < 1 || iz > 30 )
705  {
706  fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n",
707  iz );
708  cdEXIT(EXIT_FAILURE);
709  }
710  if( in < 1 || in > iz )
711  {
712  fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n",
713  in );
714  cdEXIT(EXIT_FAILURE);
715  }
716  if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) ||
717  (iz == 26 && in > 11) )
718  {
719  tt = sqrt(t/rnew[in-1][iz-1][2]);
720  rate =
721  rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])*
722  pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1]));
723  }
724  else
725  {
726  tt = t*1.0e-04;
727  if( iz == 26 && in <= 13 )
728  {
729  rate = fe[in-1][0]/pow(tt,fe[in-1][1] +
730  fe[in-1][2]*log10(tt));
731  }
732  else
733  {
734  rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]);
735  }
736  }
737 
738  return rate;
739 }
740 
741 double t_ADfA::H_rad_rec(long int iz,
742  long int n,
743  double t)
744 {
745  /*
746  * Version 4, October 9, 1997
747  ******************************************************************************
748  *** This subroutine calculates state-specific recombination rates
749  *** for hydrogen and hydrogen-like ions.
750  *** Input parameters: iz - atomic number
751  *** n - shell number, from 0 to 400:
752  *** 0 - 1s
753  *** 1 - 2s
754  *** 2 - 2p
755  *** 3 - 3
756  *** ......
757  *** t - temperature, K
758  *** Output parameter: r - rate coefficient, cm^3 s^(-1)
759  *** If n is negative, the subroutine returns the total recombination
760  *** rate coefficient
761  ******************************************************************************
762  */
763  double rate,
764  TeScaled,
765  x,
766  x1,
767  x2;
768 
769  DEBUG_ENTRY( "H_rad_rec()" );
770 
771  rate = 0.0;
772 
773  /* iz is charge, must be 1 or greater */
774  ASSERT( iz > 0 );
775 
776  /* n is level number, must be less than dim or hydro vectors */
777  ASSERT( n < NHYDRO_MAX_LEVEL );
778 
779  TeScaled = t/POW2((double)iz);
780 
781  if( n < 0 )
782  {
783  x1 = sqrt(TeScaled/3.148);
784  x2 = sqrt(TeScaled/7.036e05);
785  rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748);
786  }
787  else
788  {
789  x = log10(TeScaled);
790  rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]*
791  x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/
792  (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]*
793  powi(x,3) + HRF[n][7]*powi(x,4));
794  rate = pow(10.0,rate)/TeScaled;
795  }
796  rate *= iz;
797 
798  return rate;
799 }
800 
801 /*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
802  * returns collisional ionization rate coefficient cm^3 s^-1*/
804  /* atomic number, 1 for hydrogen */
805  long int iz,
806  /* stage of ionization, 1 for atom */
807  long int in,
808  /* temperature */
809  double t)
810 {
811  double rate, te, u;
812 
813  DEBUG_ENTRY( "coll_ion()" );
814  /*D Verner's routine to compute collisional ionization rate coefficients
815  * Version 3, April 21, 1997
816  * Cu (Z=29) and Zn (Z=30) are added (fits from Ni, correct thresholds).
817  ******************************************************************************
818  *** This subroutine calculates rates of direct collisional ionization
819  *** for all ionization stages of all elements from H to Ni (Z=28)
820  *** by use of the fits from
821  *>>refer all collion Voronov, G. S., 1997, At. Data Nucl. Data Tables, 65, 1
822  *** Input parameters: iz - atomic number on pphysical scale, H is 1
823  *** in - number of electrons from 1 to iz
824  *** t - temperature, K
825  *** Output parameter: c - rate coefficient, cm^3 s^(-1)
826  ****************************************************************************** */
827  rate = 0.0;
828 
829  if( iz < 1 || iz > 30 )
830  {
831  /* return zero rate is atomic number outside range of code */
832  return( 0. );
833  }
834 
835  if( in < 1 || in > iz )
836  {
837  /* return zero rate is ion stage is impossible */
838  return( 0. );
839  }
840 
841  te = t*EVRYD/TE1RYD;
842  u = CF[in-1][iz-1][0]/te;
843  if( u > 80.0 )
844  {
845  return( 0. );
846  }
847 
848  rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]*
849  sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])*
850  exp(-u));
851 
852  return(rate);
853 }
854 
855 realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe )
856 {
857  realnum rate;
858 
859  DEBUG_ENTRY( "h_coll_str()" );
860 
861  ASSERT( ipLo < ipHi );
862 
863 #ifndef NDEBUG
864  long ipISO = ipH_LIKE;
865  long nelem = ipHYDROGEN;
866 #endif
867  ASSERT( N_(ipLo) < N_(ipHi) );
868  ASSERT( N_(ipHi) <= 5 );
869 
870  rate = (realnum)HCS[ipHi-1][ipLo][ipTe];
871 
872  return rate;
873 }

Generated for cloudy by doxygen 1.8.1.2