cloudy  trunk
 All Data Structures Namespaces 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-2022 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 
6 #include "atmdat_adfa.h"
7 
8 #include "atmdat.h"
9 #include "iso.h"
10 
12 t_ADfA::t_ADfA()
13 {
14  DEBUG_ENTRY( "t_ADfA()" );
15 
16  /* this option, use the new atmdat_rad_rec recombination rates */
17  version = PHFIT_UNDEF;
18 
19  double help[9];
20  const long VERSION_MAGIC = 20061204L;
21 
22  static const char chFile[] = "phfit.dat";
23 
24  FILE *io = open_data( chFile, "r" );
25 
26  bool lgErr = false;
27  long i=-1, j=-1, k=-1, n;
28 
29  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
30  if( lgErr || i != VERSION_MAGIC )
31  {
32  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
33  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
34  cdEXIT(EXIT_FAILURE);
35  }
36 
37  for( n=0; n < 7; n++ )
38  lgErr = lgErr || ( fscanf( io, "%ld", &L[n] ) != 1 );
39  for( n=0; n < 30; n++ )
40  lgErr = lgErr || ( fscanf( io, "%ld", &NINN[n] ) != 1 );
41  for( n=0; n < 30; n++ )
42  lgErr = lgErr || ( fscanf( io, "%ld", &NTOT[n] ) != 1 );
43  while( true )
44  {
45  lgErr = lgErr || ( fscanf( io, "%ld %ld %ld", &i, &j, &k ) != 3 );
46  if( i == -1 && j == -1 && k == -1 )
47  break;
48  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le", &help[0], &help[1],
49  &help[2], &help[3], &help[4], &help[5] ) != 6 );
50  for( int l=0; l < 6; ++l )
51  PH1[i][j][k][l] = (realnum)help[l];
52  }
53  while( true )
54  {
55  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
56  if( i == -1 && j == -1 )
57  break;
58  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le", &help[0], &help[1],
59  &help[2], &help[3], &help[4], &help[5], &help[6] ) != 7 );
60  for( int l=0; l < 7; ++l )
61  PH2[i][j][l] = (realnum)help[l];
62  }
63  fclose( io );
64 
65  ASSERT( !lgErr );
66 
67  static const char chFile2[] = "hpfit.dat";
68 
69  io = open_data( chFile2, "r" );
70 
71  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
72  if( lgErr || i != VERSION_MAGIC )
73  {
74  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile2, i );
75  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
76  cdEXIT(EXIT_FAILURE);
77  }
78 
79  for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
80  {
81  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &help[0], &help[1],
82  &help[2], &help[3], &help[4] ) != 5 );
83  for( int l=0; l < 5; ++l )
84  PHH[i][l] = (realnum)help[l];
85  }
86 
87  fclose( io );
88 
89  ASSERT( !lgErr );
90 
91  static const char chFile3[] = "rec_lines.dat";
92 
93  io = open_data( chFile3, "r" );
94 
95  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
96  if( lgErr || i != VERSION_MAGIC )
97  {
98  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile3, i );
99  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
100  cdEXIT(EXIT_FAILURE);
101  }
102 
103  for( i=0; i < 110; i++ )
104  {
105  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
106  &help[3], &help[4], &help[5], &help[6], &help[7] ) != 8 );
107  for( int l=0; l < 8; ++l )
108  P[l][i] = (realnum)help[l];
109  }
110 
111 
112  for( i=0; i < 405; i++ )
113  {
114  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
115  &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
116  for( int l=0; l < 9; ++l )
117  ST[l][i] = (realnum)help[l];
118  }
119 
120  fclose( io );
121 
122  ASSERT( !lgErr );
123 
124  static const char chFile4[] = "rad_rec.dat";
125 
126  io = open_data( chFile4, "r" );
127 
128  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
129  if( lgErr || i != VERSION_MAGIC )
130  {
131  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile4, i );
132  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
133  cdEXIT(EXIT_FAILURE);
134  }
135 
136  while( true )
137  {
138  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
139  if( i == -1 && j == -1 )
140  break;
141  lgErr = lgErr || ( fscanf( io, "%le %le", &help[0], &help[1] ) != 2 );
142  for( int l=0; l < 2; ++l )
143  rrec[i][j][l] = (realnum)help[l];
144  }
145  while( true )
146  {
147  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
148  if( i == -1 && j == -1 )
149  break;
150  lgErr = lgErr || ( fscanf( io, "%le %le %le %le", &help[0], &help[1],
151  &help[2], &help[3] ) != 4 );
152  for( int l=0; l < 4; ++l )
153  rnew[i][j][l] = (realnum)help[l];
154  }
155  while( true )
156  {
157  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
158  if( i == -1 )
159  break;
160  lgErr = lgErr || ( fscanf( io, "%le %le %le", &help[0], &help[1], &help[2] ) != 3 );
161  for( int l=0; l < 3; ++l )
162  fe[i][l] = (realnum)help[l];
163  }
164 
165  fclose( io );
166 
167  ASSERT( !lgErr );
168 
169  static const char chFile5[] = "h_rad_rec.dat";
170 
171  io = open_data( chFile5, "r" );
172 
173  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
174  if( lgErr || i != VERSION_MAGIC )
175  {
176  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile5, i );
177  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
178  cdEXIT(EXIT_FAILURE);
179  }
180 
181  for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
182  {
183  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le %le", &help[0], &help[1], &help[2],
184  &help[3], &help[4], &help[5], &help[6], &help[7], &help[8] ) != 9 );
185  for( int l=0; l < 9; ++l )
186  HRF[i][l] = (realnum)help[l];
187  }
188 
189  fclose( io );
190 
191  ASSERT( !lgErr );
192 
193  static const char chFile6[] = "h_phot_cs.dat";
194 
195  io = open_data( chFile6, "r" );
196 
197  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
198  if( lgErr || i != VERSION_MAGIC )
199  {
200  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile6, i );
201  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
202  cdEXIT(EXIT_FAILURE);
203  }
204 
205  for( i=0; i < NHYDRO_MAX_LEVEL; i++ )
206  {
207  lgErr = lgErr || ( fscanf( io, "%le", &help[0] ) != 1 );
208  STH[i] = (realnum)help[0];
209  }
210 
211  fclose( io );
212 
213  ASSERT( !lgErr );
214 
215  static const char chFile7[] = "coll_ion.dat";
216 
217  io = open_data( chFile7, "r" );
218 
219  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
220  if( lgErr || i != VERSION_MAGIC )
221  {
222  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile7, i );
223  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
224  cdEXIT(EXIT_FAILURE);
225  }
226 
227  while( true )
228  {
229  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
230  if( i == -1 && j == -1 )
231  break;
232  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le", &CF[i][j][0], &CF[i][j][1],
233  &CF[i][j][2], &CF[i][j][3], &CF[i][j][4] ) != 5 );
234  }
235 
236  fclose( io );
237 
238  ASSERT( !lgErr );
239 
240  /*refer HI cs Anderson, H., Ballance, C.P., Badnell, N.R.,
241  *refercon & Summers, H.P 2000, J Phys B, 33, 1255 */
242  static const char chFile8[] = "h_coll_str.dat";
243 
244  io = open_data( chFile8, "r" );
245 
246  lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
247  if( lgErr || i != VERSION_MAGIC )
248  {
249  fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile8, i );
250  fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
251  cdEXIT(EXIT_FAILURE);
252  }
253 
254  while( true )
255  {
256  lgErr = lgErr || ( fscanf( io, "%ld %ld", &i, &j ) != 2 );
257  if( i == -1 && j == -1 )
258  break;
259  lgErr = lgErr || ( fscanf( io, "%le %le %le %le %le %le %le %le", &HCS[i-2][j-1][0], &HCS[i-2][j-1][1],
260  &HCS[i-2][j-1][2], &HCS[i-2][j-1][3], &HCS[i-2][j-1][4], &HCS[i-2][j-1][5],
261  &HCS[i-2][j-1][6], &HCS[i-2][j-1][7] ) != 8 );
262  }
263 
264  fclose( io );
265 
266  ASSERT( !lgErr );
267 }
268 
269 double t_ADfA::phfit(long int nz,
270  long int ne,
271  long int is,
272  double e)
273 {
274  long int nint,
275  nout;
276  double a,
277  b,
278  crs,
279  einn,
280  p1,
281  q,
282  x,
283  y,
284  z;
285 
286  DEBUG_ENTRY( "phfit()" );
287 
288  /*** Version 3. October 8, 1996.
289  *** Written by D. A. Verner, verner@pa.uky.edu
290  *** Inner-shell ionization energies of some low-ionized species are slightly
291  *** improved to fit smoothly the experimental inner-shell ionization energies
292  *** of neutral atoms.
293  ******************************************************************************
294  *** This subroutine calculates partial photoionization cross sections
295  *** for all ionization stages of all atoms from H to Zn (Z=30) by use of
296  *** the following fit parameters:
297  *** Outer shells of the Opacity Project (OP) elements:
298  *** >>refer all photo_cs Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487.
299  *** Inner shells of all elements, and outer shells of the non-OP elements:
300  *** Verner and Yakovlev, 1995, A&AS, 109, 125
301  *** Input parameters: nz - atomic number from 1 to 30 (integer)
302  *** ne - number of electrons from 1 to iz (integer)
303  *** is - shell number (integer)
304  *** e - photon energy, eV
305  *** version - enum, PHFIT96 (default): calculates
306  *** new cross sections, PHFIT95: calculates
307  *** only old Hartree-Slater cross sections
308  *** Output parameter: photoionization cross section, Mb
309  *** Shell numbers:
310  *** 1 - 1s, 2 - 2s, 3 - 2p, 4 - 3s, 5 - 3p, 6 - 3d, 7 - 4s.
311  *** If a species in the ground state has no electrons on the given shell,
312  *** the subroutine returns 0.
313  ****************************************************************************** */
314 
315  crs = 0.0;
316  if( nz < 1 || nz > 30 )
317  {
318  return crs;
319  }
320 
321  if( ne < 1 || ne > nz )
322  {
323  return crs;
324  }
325 
326  nout = NTOT[ne-1];
327  if( nz == ne && nz > 18 )
328  nout = 7;
329  if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) ||
330  nz == 25) || nz == 26) )
331  nout = 7;
332  if( is > nout )
333  {
334  return crs;
335  }
336 
337  if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 )
338  {
339  return crs;
340  }
341 
342  ASSERT( is >= 1 && is <= 7 );
343 
344  if( e < PH1[is-1][ne-1][nz-1][0] )
345  {
346  return crs;
347  }
348 
349  nint = NINN[ne-1];
350  if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) )
351  {
352  einn = 0.0;
353  }
354  else
355  {
356  if( ne < 3 )
357  {
358  einn = 1.0e30;
359  }
360  else
361  {
362  einn = PH1[nint-1][ne-1][nz-1][0];
363  }
364  }
365 
366  if( (is <= nint || e >= einn) || version == PHFIT95 )
367  {
368  p1 = -PH1[is-1][ne-1][nz-1][4];
369  y = e/PH1[is-1][ne-1][nz-1][1];
370  q = -0.5*p1 - L[is-1] - 5.5;
371  a = PH1[is-1][ne-1][nz-1][2]*(POW2(y - 1.0) +
372  POW2(PH1[is-1][ne-1][nz-1][5]));
373  b = sqrt(y/PH1[is-1][ne-1][nz-1][3]) + 1.0;
374  crs = a*pow(y,q)*pow(b,p1);
375  }
376  else
377  {
378  if( (is < nout && is > nint) && e < einn )
379  {
380  return crs;
381  }
382  p1 = -PH2[ne-1][nz-1][3];
383  q = -0.5*p1 - 5.5;
384  x = e/PH2[ne-1][nz-1][0] - PH2[ne-1][nz-1][5];
385  z = sqrt(x*x+POW2(PH2[ne-1][nz-1][6]));
386  a = PH2[ne-1][nz-1][1]*(POW2(x - 1.0) +
387  POW2(PH2[ne-1][nz-1][4]));
388  b = 1.0 + sqrt(z/PH2[ne-1][nz-1][2]);
389  crs = a*pow(z,q)*pow(b,p1);
390  }
391  return crs;
392 }
393 
394 double t_ADfA::hpfit(long int iz,
395  long int n,
396  double e)
397 {
398  long int l,
399  m;
400  double cs,
401  eth,
402  ex,
403  q,
404  x;
405 
406  DEBUG_ENTRY( "hpfit()" );
407 
408  /*state specific photoionization cross sections for model hydrogen atom
409  * Version 1, September 23, 1997
410  ******************************************************************************
411  *** This subroutine calculates state-specific photoionization cross sections
412  *** for hydrogen and hydrogen-like ions.
413  *** Input parameters: iz - atomic number
414  *** n - shell number, from 0 to 400:
415  *** 0 - 1s
416  *** 1 - 2s
417  *** 2 - 2p
418  *** 3 - 3
419  *** ......
420  *** e - photon energy, eV
421  *** return value - cross section, cm^(-2)
422  *******************************************************************************/
423 
424  ASSERT( iz > 0 && e>0. );
425 
426  if( n >= NHYDRO_MAX_LEVEL )
427  {
428  fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n );
429  cdEXIT(EXIT_FAILURE);
430  }
431 
432  l = 0;
433  if( n == 2 )
434  {
435  l = 1;
436  }
437  q = 3.5 + l - 0.5*PHH[n][1];
438 
439  if( n == 0 )
440  {
441  m = 1;
442  }
443  else
444  {
445  if( n == 1 )
446  {
447  m = 2;
448  }
449  else
450  {
451  m = n;
452  }
453  }
454 
455  eth = ph1(0,0,iz-1,0)/POW2((double)m);
456  ex = MAX2(1. , e/eth );
457 
458  /* Don't just force to be at least one...make sure e/eth is close to one or greater. */
459  ASSERT( e/eth > 0.95 );
460 
461  if( ex < 1.0 )
462  {
463  return 0.;
464  }
465 
466  x = ex/PHH[n][0];
467  cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/
468  pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/
469  POW2((double)iz));
470  return cs;
471 }
472 
473 void t_ADfA::rec_lines(double t,
474  realnum r[][NRECCOEFCNO])
475 {
476  long int i,
477  j,
478  ipj1,
479  ipj2;
480 
481  double a,
482  a1,
483  dr[4][405],
484  p1,
485  p2,
486  p3,
487  p4,
488  p5,
489  p6,
490  rr[4][110],
491  te,
492  x,
493  z;
494 
495  static long jd[6]={143,145,157,360,376,379};
496 
497  static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50,
498  52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100,
499  101,103,104};
500 
501  static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112,
502  120,114,143,145,157,152,169,183,200,163,225,223,237,232,235,
503  249,247,300,276,278,376,360,379,384};
504 
505  DEBUG_ENTRY( "rec_lines()" );
506 
507  /*effective recombination coefficients for lines of C, N, O, by D. Verner
508  * Version 2, April 30, 1997
509  ******************************************************************************
510  *** This subroutine calculates effective recombination coefficients
511  *** for 110 permitted recombination lines of C, N, O (Pequignot, Petitjean,
512  *** & Boisson, 1991, A&A, 251, 680) and 405 permitted dielectronic
513  *** recombination lines (Nussbaumer & Storey, 1984, A&AS, 56, 293)
514  *** Input parameter: t - temperature, K
515  *** Output parameters: r(i,j), i=1,471
516  *** r(i,1) - atomic number
517  *** r(i,2) - number of electrons
518  *** r(i,3) - wavelength, angstrom
519  *** r(i,4) - rate coefficient, cm^3 s^(-1)
520  ****************************************************************************** */
521 
522  for( i=0; i < 110; i++ )
523  {
524  rr[0][i] = P[0][i];
525  rr[1][i] = P[1][i];
526  rr[2][i] = P[2][i];
527  z = P[0][i] - P[1][i] + 1.0;
528  te = 1.0e-04*t/z/z;
529  p1 = P[3][i];
530  p2 = P[4][i];
531  p3 = P[5][i];
532  p4 = P[6][i];
533  if( te < 0.004 )
534  {
535  a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
536  a = a1/sqrt(te/0.004);
537  }
538  else
539  {
540  if( te > 2.0 )
541  {
542  a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
543  a = a1/pow(te/2.0,1.5);
544  }
545  else
546  {
547  a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
548  }
549  }
550  rr[3][i] = 1.0e-13*z*a*P[7][i];
551  }
552 
553  for( i=0; i < 405; i++ )
554  {
555  dr[0][i] = ST[0][i];
556  dr[1][i] = ST[1][i];
557  dr[2][i] = ST[2][i];
558  te = 1.0e-04*t;
559  p1 = ST[3][i];
560  p2 = ST[4][i];
561  p3 = ST[5][i];
562  p4 = ST[6][i];
563  p5 = ST[7][i];
564  p6 = ST[8][i];
565  if( te < p6 )
566  {
567  x = p5*(1.0/te - 1.0/p6);
568  if( x > 80.0 )
569  {
570  a = 0.0;
571  }
572  else
573  {
574  a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/powpq(p6,3,2)/exp(p5/p6);
575  a = a1/exp(x);
576  }
577  }
578  else
579  {
580  if( te > 6.0 )
581  {
582  a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/powpq(6.0,3,2)/exp(p5/6.0);
583  a = a1/powpq(te/6.0,3,2);
584  }
585  else
586  {
587  a = (p1/te + p2 + p3*te + p4*te*te)/powpq(te,3,2)/exp(p5/te);
588  }
589  }
590  dr[3][i] = 1.0e-12*a;
591  }
592 
593  for( i=0; i < 6; i++ )
594  {
595  ipj1 = jd[i];
596  ipj2 = ipj1 + 1;
597  dr[3][ipj1-1] += dr[3][ipj2-1];
598  dr[0][ipj2-1] = 0.0;
599  }
600 
601  for( i=0; i < 38; i++ )
602  {
603  ipj1 = ip[i];
604  ipj2 = id[i];
605  rr[3][ipj1-1] += dr[3][ipj2-1];
606  dr[0][ipj2-1] = 0.0;
607  }
608 
609  for( i=0; i < 110; i++ )
610  {
611  r[0][i] = (realnum)rr[0][i];
612  r[1][i] = (realnum)rr[1][i];
613  r[2][i] = (realnum)rr[2][i];
614  r[3][i] = (realnum)rr[3][i];
615  }
616 
617  j = 110;
618  for( i=0; i < 405; i++ )
619  {
620  if( dr[0][i] > 1.0 )
621  {
622  j += 1;
623  r[0][j-1] = (realnum)dr[0][i];
624  r[1][j-1] = (realnum)dr[1][i];
625  r[2][j-1] = (realnum)dr[2][i];
626  r[3][j-1] = (realnum)dr[3][i];
627  }
628  }
629  return;
630 }
631 
632 double t_ADfA::rad_rec(long int iz,
633  long int in,
634  double t)
635 {
636  /*
637  *** Version 4. June 29, 1999.
638  *** Written by D. A. Verner, verner@pa.uky.edu
639  ******************************************************************************
640  *** This subroutine calculates rates of radiative recombination for all ions
641  *** of all elements from H through Zn by use of the following fits:
642  *** H-like, He-like, Li-like, Na-like -
643  *** >>refer all rec_coef Verner, D. A., & Ferland, G. J. 1996, ApJS, 103, 467
644  *** Other ions of C, N, O, Ne - Pequignot et al. 1991, A&A, 251, 680,
645  *** refitted by Verner & Ferland formula to ensure correct asymptotes
646  *** Fe XVII-XXIII -
647  *** >>refer Fe17-23 recom Arnaud, M., & Raymond, J. 1992, ApJ, 398, 394
648  *** Fe I-XV - refitted by Verner & Ferland formula to ensure correct asymptotes
649  *** Other ions of Mg, Si, S, Ar, Ca, Fe, Ni -
650  *** -
651  *** >>refer all recom Shull, M. J., & Van Steenberg, M. 1982, ApJS, 48, 95
652  *** Other ions of Na, Al -
653  *** >>refer Na, Al recom Landini, M., & Monsignori-Fossi, B.C. 1990, A&AS, 82, 229
654  *** Other ions of F, P, Cl, K, Ti, Cr, Mn, Co (excluding Ti I-II, Cr I-IV,
655  *** Mn I-V, Co I) -
656  *** >>refer many recom Landini, M., & Monsignori-Fossi, B.C. 1991, A&AS, 91, 183
657  *** All other species - interpolations of the power-law fits
658  *** Input parameters: iz - atomic number
659  *** in - number of electrons from 1 to iz
660  *** t - temperature, K
661  *** return result: - rate coefficient, cm^3 s^(-1)
662  ******************************************************************************
663  */
664  double tt;
665  double rate;
666 
667  DEBUG_ENTRY( "rad_rec()" );
668 
669  if( iz < 1 || iz > 30 )
670  {
671  fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n",
672  iz );
673  cdEXIT(EXIT_FAILURE);
674  }
675  if( in < 1 || in > iz )
676  {
677  fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n",
678  in );
679  cdEXIT(EXIT_FAILURE);
680  }
681  if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) ||
682  (iz == 26 && in > 11) )
683  {
684  tt = sqrt(t/rnew[in-1][iz-1][2]);
685  rate =
686  rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])*
687  pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1]));
688  }
689  else
690  {
691  tt = t*1.0e-04;
692  if( iz == 26 && in <= 13 )
693  {
694  rate = fe[in-1][0]/pow(tt,fe[in-1][1] +
695  fe[in-1][2]*log10(tt));
696  }
697  else
698  {
699  rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]);
700  }
701  }
702 
703  return rate;
704 }
705 
706 double t_ADfA::H_rad_rec(long int iz,
707  long int n,
708  double t)
709 {
710  /*
711  * Version 4, October 9, 1997
712  ******************************************************************************
713  *** This subroutine calculates state-specific recombination rates
714  *** for hydrogen and hydrogen-like ions.
715  *** Input parameters: iz - atomic number
716  *** n - shell number, from 0 to 400:
717  *** 0 - 1s
718  *** 1 - 2s
719  *** 2 - 2p
720  *** 3 - 3
721  *** ......
722  *** t - temperature, K
723  *** Output parameter: r - rate coefficient, cm^3 s^(-1)
724  *** If n is negative, the subroutine returns the total recombination
725  *** rate coefficient
726  ******************************************************************************
727  */
728  double rate,
729  TeScaled,
730  x,
731  x1,
732  x2;
733 
734  DEBUG_ENTRY( "H_rad_rec()" );
735 
736  /* iz is charge, must be 1 or greater */
737  ASSERT( iz > 0 );
738 
739  /* n is level number, must be less than dim or hydro vectors */
740  ASSERT( n < NHYDRO_MAX_LEVEL );
741 
742  TeScaled = t/POW2((double)iz);
743 
744  if( n < 0 )
745  {
746  x1 = sqrt(TeScaled/3.148);
747  x2 = sqrt(TeScaled/7.036e05);
748  rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748);
749  }
750  else
751  {
752  x = log10(TeScaled);
753  rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]*
754  x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/
755  (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]*
756  powi(x,3) + HRF[n][7]*powi(x,4));
757  rate = exp10(rate)/TeScaled;
758  }
759  rate *= iz;
760 
761  return rate;
762 }
763 
764 /*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
765  * returns collisional ionization rate coefficient cm^3 s^-1*/
766 double t_ADfA::coll_ion(
767  /* atomic number, 1 for hydrogen */
768  long int iz,
769  /* number of bound electrons before ionization*/
770  long int in,
771  /* temperature */
772  double t)
773 {
774  double rate, te, u;
775 
776  DEBUG_ENTRY( "coll_ion()" );
777  /*D Verner's routine to compute collisional ionization rate coefficients
778  * Version 3, April 21, 1997
779  * Cu (Z=29) and Zn (Z=30) are added (fits from Ni, correct thresholds).
780  ******************************************************************************
781  *** This subroutine calculates rates of direct collisional ionization
782  *** for all ionization stages of all elements from H to Ni (Z=28)
783  *** by use of the fits from
784  *>>refer all coll_ion Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1
785  *** Input parameters: iz - atomic number on pphysical scale, H is 1
786  *** in - number of electrons from 1 to iz
787  *** t - temperature, K
788  *** Output parameter: c - rate coefficient, cm^3 s^(-1)
789  ****************************************************************************** */
790 
791  te = t*EVRYD/TE1RYD;
792  u = CF[in-1][iz-1][0]/te;
793  if( u > 80.0 )
794  {
795  return 0.;
796  }
797 
798  rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]*
799  sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])*
800  exp(-u)) * atmdat.lgCollIonOn;
801 
802  return rate;
803 }
804 
805 double t_ADfA::coll_ion_wrapper(
806  /* (atomic number - 1), 0 for hydrogen */
807  long int z,
808  /* stage of ionization, 0 for atom */
809  long int n,
810  /* temperature K */
811  double t)
812 {
813  double rate = 0.0;
814 
815  static const bool DEBUG_COLL_ION = false;
816 
817  DEBUG_ENTRY( "coll_ion_wrapper()" );
818 
819  if( z < 0 || z > LIMELM-1 )
820  {
821  /* return zero rate is atomic number outside range of code */
822  return 0.;
823  }
824 
825  if( n < 0 || n > z )
826  {
827  /* return zero rate is ion stage is impossible */
828  return 0.;
829  }
830 
831  /*Get the rate coefficients from either Dima or Hybrid.
832  * The enum variable CIRCData is set to a default in init_defaults_preparse
833  * It can be changed parse_set via the input command set collisional ionization XXX */
834  if( atmdat.CIRCData == t_atmdat::DIMA)
835  {
836  /* collisional ionization by thermal electrons
837  * >>chng 97 mar 19, to Dima's new routine using
838  * >>refer all coll_ion Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1
839  *coll_ion(atomic number, number of electrons, temperature)*/
840  rate = coll_ion( z+1, z+1-n, t );
841 
842  if( DEBUG_COLL_ION )
843  {
844  fprintf(ioQQQ,"DIMA:%li\t%li\t%e\t%e\n",z,n,t,rate);
845  }
846  }
847  else if( atmdat.CIRCData == t_atmdat::HYBRID)
848  {
849  // Get the rate coefficient by scaling Dima to Dere.
850  rate = coll_ion_hybrid( z, n, t );
851 
852  if( DEBUG_COLL_ION )
853  {
854  fprintf(ioQQQ,"HYBRID:%li\t%li\t%e\t%e\n",z,n,t,rate);
855  }
856  }
857  else
858  {
859  //CIRCData is an enum so it must be equal to one of the enum values.
860  TotalInsanity();
861  }
862 
863  ASSERT( rate >= 0.0 );
864  return rate;
865 }
866 
867 /*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
868  * returns collisional ionization rate coefficient cm^3 s^-1*/
869 double t_ADfA::coll_ion_hybrid(
870  /* (atomic number - 1) 0 for hydrogen */
871  long int nelem,
872  /* stage of ionization, 0 for atom */
873  long int ion,
874  /* temperature */
875  double t)
876 {
877 
878  DEBUG_ENTRY( "coll_ion_hybrid()" );
879 
883  static const double DereRatio[30][30]=
884  {{0.9063,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
885  {1.0389,1.0686,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
886  {0.6466,1.0772,0.963,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
887  {0.9715,0.7079,0.973,0.9899,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
888  {1.0781,1.3336,1.0032,1.1102,0.9874,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
889  {1.0499,0.913,1.0377,1.299,1.2874,1.0656,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
890  {1.0421,1.1966,1.0842,0.9619,1.0583,1.155,1.0648,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
891  {1.041,1.1181,1.0164,1.0271,0.9789,0.6829,1.1805,1.0672,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
892  {0.2636,0.658,1.0431,1.0749,0.894,1.059,0.9888,1.1426,1.0633,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
893  {0.8089,1.1395,1.1041,1.0449,1.1365,1.089,1.0147,0.9135,1.336,1.0429,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
894  {0.9545,1.1302,1.1105,1.2067,1.2569,1.079,0.8866,0.9924,0.9933,1.0488,1.0602,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
895  {0.3793,0.9857,1.1152,1.2826,1.3006,1.2416,1.0893,0.93,0.9953,0.9992,1.3479,1.0622,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
896  {0.7458,0.9975,1.0251,0.9553,1.5023,1.2136,1.2566,1.2417,1.2381,1.3755,1.1113,1.1629,1.0589,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
897  {0.7328,0.8798,0.4492,0.8221,1.7874,1.5418,1.5777,1.3036,1.124,1.0667,1.0009,1.002,1.1284,1.0673,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
898  {0.7075,0.9805,0.9451,0.5937,1.2061,1.1772,1.3818,1.4073,1.2643,1.1095,0.9903,1.3716,1.0058,1.0966,1.0517,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
899  {1.3572,0.8925,0.8119,0.8991,0.6633,1.4519,1.2073,1.4058,1.4519,1.2726,1.1428,0.9818,1.3643,1.0074,1.1836,1.1005,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
900  {0.5412,0.8428,0.9237,0.819,0.9151,0.7909,1.7424,1.2653,1.4513,1.4871,1.2633,1.1373,0.9969,0.9919,1.0091,1.2077,1.105,0,0,0,0,0,0,0,0,0,0,0,0,0},
901  {0.9242,0.8644,0.9752,0.8562,0.6998,0.6273,0.8686,2.0326,1.0364,1.4687,1.4812,1.2701,1.1376,0.9932,1.0001,1.0084,1.2036,1.1074,0,0,0,0,0,0,0,0,0,0,0,0},
902  {0.6381,0.9604,1.3556,0.9364,0.9678,1.0604,1.3067,1.0434,2.0722,0.979,1.5318,1.4968,1.2713,1.1427,1.0036,0.9953,1.03,1.2071,1.1083,0,0,0,0,0,0,0,0,0,0,0},
903  {0.7652,1.1668,1.0422,0.8705,0.9193,0.9982,1.077,1.3875,1.1544,2.4653,1.3188,1.5375,1.5307,1.2776,1.1427,1.0117,0.9872,1.0118,1.1899,1.1119,0,0,0,0,0,0,0,0,0,0},
904  {1.0951,0.5891,0.9222,1.2903,1.287,1.0563,1.0444,1.1405,1.4841,1.2583,2.7347,0.9844,1.034,1.0412,1.1062,1.2211,1.6728,1.6347,1.6554,1.8095,1.9561,0,0,0,0,0,0,0,0,0},
905  {0.9789,0.7389,1.3107,0.9274,0.9921,0.9812,0.8971,0.9936,1.0079,1.0377,1.2073,2.7198,0.9995,1.037,1.0433,1.1093,1.2323,1.6673,1.6231,1.6666,1.8089,1.7302,0,0,0,0,0,0,0,0},
906  {0.6169,0.4618,1.3566,3.3812,1.1951,1.1746,1.0133,0.9195,1.0548,1.0608,1.1016,1.285,3.1081,0.9986,1.0094,1.0379,1.005,0.9932,1.0651,0.8787,0.952,0.9862,1.0083,0,0,0,0,0,0,0},
907  {1.0603,0.262,0.9237,2.4323,1.8345,1.2197,1.0924,1.0205,0.9379,1.0946,1.1001,1.1741,1.3608,3.152,0.9692,0.8866,1.0556,1.1127,1.2289,1.6828,1.6298,1.6469,1.8082,1.0605,0,0,0,0,0,0},
908  {0.9061,0.7287,0.9676,1.9425,1.5785,1.9028,1.3827,1.1136,1.0368,0.9516,1.1389,1.1594,1.1642,1.4161,2.8162,0.9744,0.9246,1.0644,1.1118,1.2332,1.6892,1.6311,1.6347,1.8073,1.0722,0,0,0,0,0},
909  {0.9904,1.0568,1.824,1.6578,1.5033,0.9704,0.8933,0.8579,1.084,1.0308,0.9637,1.1885,1.2179,1.2979,1.4846,3.0344,0.9879,0.8927,1.0534,1.1233,1.2242,1.6794,1.6255,1.6487,1.8141,1.7629,0,0,0,0},
910  {0.9667,1.2622,0.959,2.8444,1.304,1.5632,1.709,2.298,1.427,1.0885,1.0285,0.9767,1.2237,1.2632,1.3585,1.5495,3.1385,0.9959,0.9327,1.0904,1.0357,1.0125,1.0027,0.9788,0.8975,1.0539,1.048,0,0,0},
911  {1.0004,0.8721,1.3073,0.9564,2.7856,1.6197,1.5516,2.229,2.1494,1.4916,1.0871,1.0359,0.9768,1.2269,1.3265,1.4169,1.597,3.4077,0.9979,0.8869,1.0635,1.1063,1.2267,1.6914,1.6401,1.6216,1.0598,1.0363,0,0},
912  {0.728,1.3939,0.4822,0.626,0.5463,2.2491,1.064,1.1998,1.3083,1.6545,1.1149,0.8826,0.8566,0.8196,1.1173,1.17,1.2445,1.394,2.9832,0.8715,0.7703,0.929,0.968,1.0828,1.4973,1.4513,1.4476,0.9621,0.9353,0},
913  {1.0766,0.6459,0.7688,0.2358,0.3709,0.3476,1.6754,0.8288,0.9184,1.0686,1.4546,0.9515,0.7272,0.7048,0.6911,0.9722,1.0308,1.0988,1.1959,2.6404,0.7644,0.6768,0.8181,0.8566,0.9689,1.3344,1.3031,1.3345,0.8676,0.8478}};
914  /* Get the hybrid coeffs from the Dima rate coefficient times the average Dere to Dima ratio.
915  */
916  ASSERT( nelem>=0 && nelem<LIMELM && ion>=0 && ion<=nelem );
917  double rate = coll_ion(nelem+1,nelem+1-ion,t) * DereRatio[nelem][ion];
918  ASSERT( rate >=0. );
919  return rate;
920 }
921 realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe )
922 {
923  realnum rate;
924 
925  DEBUG_ENTRY( "h_coll_str()" );
926 
927  ASSERT( ipLo < ipHi );
928 
929 # if !defined NDEBUG
930  long ipISO = ipH_LIKE;
931  long nelem = ipHYDROGEN;
932 # endif
933  ASSERT( N_(ipLo) < N_(ipHi) );
934  ASSERT( N_(ipHi) <= 5 );
935 
936  rate = (realnum)HCS[ipHi-1][ipLo][ipTe];
937 
938  return rate;
939 }
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
atmdat_adfa.h
x2
static double x2[63]
Definition: atmdat_3body.cpp:18
exp10
double exp10(double x)
Definition: cddefines.h:1368
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
x1
static double x1[83]
Definition: atmdat_3body.cpp:27
t_ADfA::h_coll_str
realnum h_coll_str(long ipLo, long ipHi, long ipTe)
Definition: atmdat_adfa.cpp:921
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat_adfa.h:61
t_ADfA::fe
realnum fe[13][3]
Definition: atmdat_adfa.h:31
NRECCOEFCNO
const int NRECCOEFCNO
Definition: atmdat_adfa.h:7
t_ADfA::PHH
realnum PHH[NHYDRO_MAX_LEVEL][5]
Definition: atmdat_adfa.h:24
t_ADfA::HRF
realnum HRF[NHYDRO_MAX_LEVEL][9]
Definition: atmdat_adfa.h:33
t_ADfA::coll_ion_wrapper
double coll_ion_wrapper(long int z, long int n, double t)
Definition: atmdat_adfa.cpp:805
t_atmdat::CIRCData
CollIonRC CIRCData
Definition: atmdat.h:443
t_atmdat::lgCollIonOn
bool lgCollIonOn
Definition: atmdat.h:349
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_ADfA::PH2
realnum PH2[30][30][7]
Definition: atmdat_adfa.h:22
t_ADfA::L
long int L[7]
Definition: atmdat_adfa.h:18
t_ADfA::NTOT
long int NTOT[30]
Definition: atmdat_adfa.h:20
t_ADfA::coll_ion
double coll_ion(long int iz, long int in, double t)
Definition: atmdat_adfa.cpp:766
t_ADfA::version
phfit_version version
Definition: atmdat_adfa.h:16
t_ADfA::ST
realnum ST[9][405]
Definition: atmdat_adfa.h:27
POW2
#define POW2
Definition: cddefines.h:979
N_
#define N_(A_)
Definition: iso.h:22
PHFIT_UNDEF
Definition: atmdat_adfa.h:8
t_ADfA::PH1
realnum PH1[7][30][30][6]
Definition: atmdat_adfa.h:21
t_ADfA::rrec
realnum rrec[30][30][2]
Definition: atmdat_adfa.h:29
realnum
float realnum
Definition: cddefines.h:124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
t_ADfA::CF
double CF[30][30][5]
Definition: atmdat_adfa.h:39
t_ADfA::H_rad_rec
double H_rad_rec(long int iz, long int n, double t)
Definition: atmdat_adfa.cpp:706
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
a1
static double a1[83]
Definition: atmdat_3body.cpp:27
powi
double powi(double, long int)
Definition: service.cpp:594
t_atmdat::DIMA
Definition: atmdat.h:442
t_ADfA::NINN
long int NINN[30]
Definition: atmdat_adfa.h:19
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
PHFIT95
Definition: atmdat_adfa.h:8
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
LIMELM
const int LIMELM
Definition: cddefines.h:308
t_ADfA::rad_rec
double rad_rec(long int iz, long int in, double t)
Definition: atmdat_adfa.cpp:632
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:630
t_ADfA::hpfit
double hpfit(long int iz, long int n, double e)
Definition: atmdat_adfa.cpp:394
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
t_atmdat::HYBRID
Definition: atmdat.h:442
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
t_ADfA::HCS
double HCS[14][10][8]
Definition: atmdat_adfa.h:45
t_ADfA::coll_ion_hybrid
double coll_ion_hybrid(long int z, long int n, double t)
Definition: atmdat_adfa.cpp:869
nint
long nint(double x)
Definition: cddefines.h:758
t_ADfA::rec_lines
void rec_lines(double t, realnum r[][NRECCOEFCNO])
Definition: atmdat_adfa.cpp:473
t_ADfA::STH
realnum STH[NHYDRO_MAX_LEVEL]
Definition: atmdat_adfa.h:37
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:349
t_ADfA::rnew
realnum rnew[30][30][4]
Definition: atmdat_adfa.h:30
t_ADfA::phfit
double phfit(long int nz, long int ne, long int is, double e)
Definition: atmdat_adfa.cpp:269
NHYDRO_MAX_LEVEL
const int NHYDRO_MAX_LEVEL
Definition: cddefines.h:316
t_ADfA::t_ADfA
t_ADfA()
Definition: atmdat_adfa.cpp:12
t_ADfA::P
realnum P[8][110]
Definition: atmdat_adfa.h:26
Generated by   doxygen 1.8.5