cloudy trunk
Loading...
Searching...
No Matches
atmdat_adfa.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/*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#include "taulines.h"
9#include "thirdparty.h"
10
13{
14 DEBUG_ENTRY( "t_ADfA()" );
15
16 /* this option, use the new atmdat_rad_rec recombination rates */
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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
269double 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
394double 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 cs = 0.0;
425
426 ASSERT( iz > 0 && e>0. );
427
428 if( n >= NHYDRO_MAX_LEVEL )
429 {
430 fprintf( ioQQQ, " hpfit called with too large n, =%li\n" , n );
432 }
433
434 l = 0;
435 if( n == 2 )
436 {
437 l = 1;
438 }
439 q = 3.5 + l - 0.5*PHH[n][1];
440
441 if( n == 0 )
442 {
443 m = 1;
444 }
445 else
446 {
447 if( n == 1 )
448 {
449 m = 2;
450 }
451 else
452 {
453 m = n;
454 }
455 }
456
457 eth = ph1(0,0,iz-1,0)/POW2((double)m);
458 ex = MAX2(1. , e/eth );
459
460 /* Don't just force to be at least one...make sure e/eth is close to one or greater. */
461 ASSERT( e/eth > 0.95 );
462
463 if( ex < 1.0 )
464 {
465 return(0.);
466 }
467
468 x = ex/PHH[n][0];
469 cs = (PHH[n][4]*pow(1.0 + ((double)PHH[n][2]/x),(double)PHH[n][3])/
470 pow(x,q)/pow(1.0 + sqrt(x),(double)PHH[n][1])*8.79737e-17/
471 POW2((double)iz));
472 return(cs);
473}
474
475void t_ADfA::rec_lines(double t,
476 realnum r[][471])
477{
478 long int i,
479 j,
480 ipj1,
481 ipj2;
482
483 double a,
484 a1,
485 dr[4][405],
486 p1,
487 p2,
488 p3,
489 p4,
490 p5,
491 p6,
492 rr[4][110],
493 te,
494 x,
495 z;
496
497 static long jd[6]={143,145,157,360,376,379};
498
499 static long ip[38]={7,9,12,13,14,16,18,19,20,21,22,44,45,49,50,
500 52,53,54,55,56,57,58,59,60,66,67,78,83,84,87,88,95,96,97,100,
501 101,103,104};
502
503 static long id[38]={7,3,10,27,23,49,51,64,38,47,60,103,101,112,
504 120,114,143,145,157,152,169,183,200,163,225,223,237,232,235,
505 249,247,300,276,278,376,360,379,384};
506
507 DEBUG_ENTRY( "rec_lines()" );
508
509 /*effective recombination coefficients for lines of C, N, O, by D. Verner
510 * Version 2, April 30, 1997
511 ******************************************************************************
512 *** This subroutine calculates effective recombination coefficients
513 *** for 110 permitted recombination lines of C, N, O (Pequignot, Petitjean,
514 *** & Boisson, 1991, A&A, 251, 680) and 405 permitted dielectronic
515 *** recombination lines (Nussbaumer & Storey, 1984, A&AS, 56, 293)
516 *** Input parameter: t - temperature, K
517 *** Output parameters: r(i,j), i=1,471
518 *** r(i,1) - atomic number
519 *** r(i,2) - number of electrons
520 *** r(i,3) - wavelength, angstrom
521 *** r(i,4) - rate coefficient, cm^3 s^(-1)
522 ****************************************************************************** */
523
524 for( i=0; i < 110; i++ )
525 {
526 rr[0][i] = P[0][i];
527 rr[1][i] = P[1][i];
528 rr[2][i] = P[2][i];
529 z = P[0][i] - P[1][i] + 1.0;
530 te = 1.0e-04*t/z/z;
531 p1 = P[3][i];
532 p2 = P[4][i];
533 p3 = P[5][i];
534 p4 = P[6][i];
535 if( te < 0.004 )
536 {
537 a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
538 a = a1/sqrt(te/0.004);
539 }
540 else
541 {
542 if( te > 2.0 )
543 {
544 a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
545 a = a1/pow(te/2.0,1.5);
546 }
547 else
548 {
549 a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
550 }
551 }
552 rr[3][i] = 1.0e-13*z*a*P[7][i];
553 }
554
555 for( i=0; i < 405; i++ )
556 {
557 dr[0][i] = ST[0][i];
558 dr[1][i] = ST[1][i];
559 dr[2][i] = ST[2][i];
560 te = 1.0e-04*t;
561 p1 = ST[3][i];
562 p2 = ST[4][i];
563 p3 = ST[5][i];
564 p4 = ST[6][i];
565 p5 = ST[7][i];
566 p6 = ST[8][i];
567 if( te < p6 )
568 {
569 x = p5*(1.0/te - 1.0/p6);
570 if( x > 80.0 )
571 {
572 a = 0.0;
573 }
574 else
575 {
576 a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/
577 p6);
578 a = a1/exp(x);
579 }
580 }
581 else
582 {
583 if( te > 6.0 )
584 {
585 a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/
586 exp(p5/6.0);
587 a = a1/pow(te/6.0,1.5);
588 }
589 else
590 {
591 a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/
592 te);
593 }
594 }
595 dr[3][i] = 1.0e-12*a;
596 }
597
598 for( i=0; i < 6; i++ )
599 {
600 ipj1 = jd[i];
601 ipj2 = ipj1 + 1;
602 dr[3][ipj1-1] += dr[3][ipj2-1];
603 dr[0][ipj2-1] = 0.0;
604 }
605
606 for( i=0; i < 38; i++ )
607 {
608 ipj1 = ip[i];
609 ipj2 = id[i];
610 rr[3][ipj1-1] += dr[3][ipj2-1];
611 dr[0][ipj2-1] = 0.0;
612 }
613
614 for( i=0; i < 110; i++ )
615 {
616 r[0][i] = (realnum)rr[0][i];
617 r[1][i] = (realnum)rr[1][i];
618 r[2][i] = (realnum)rr[2][i];
619 r[3][i] = (realnum)rr[3][i];
620 }
621
622 j = 110;
623 for( i=0; i < 405; i++ )
624 {
625 if( dr[0][i] > 1.0 )
626 {
627 j += 1;
628 r[0][j-1] = (realnum)dr[0][i];
629 r[1][j-1] = (realnum)dr[1][i];
630 r[2][j-1] = (realnum)dr[2][i];
631 r[3][j-1] = (realnum)dr[3][i];
632 }
633 }
634 return;
635}
636
637double t_ADfA::rad_rec(long int iz,
638 long int in,
639 double t)
640{
641 /*
642 *** Version 4. June 29, 1999.
643 *** Written by D. A. Verner, verner@pa.uky.edu
644 ******************************************************************************
645 *** This subroutine calculates rates of radiative recombination for all ions
646 *** of all elements from H through Zn by use of the following fits:
647 *** H-like, He-like, Li-like, Na-like -
648 *** >>refer all rec_coef Verner, D. A., & Ferland, G. J. 1996, ApJS, 103, 467
649 *** Other ions of C, N, O, Ne - Pequignot et al. 1991, A&A, 251, 680,
650 *** refitted by Verner & Ferland formula to ensure correct asymptotes
651 *** Fe XVII-XXIII -
652 *** >>refer Fe17-23 recom Arnaud, M., & Raymond, J. 1992, ApJ, 398, 394
653 *** Fe I-XV - refitted by Verner & Ferland formula to ensure correct asymptotes
654 *** Other ions of Mg, Si, S, Ar, Ca, Fe, Ni -
655 *** -
656 *** >>refer all recom Shull, M. J., & Van Steenberg, M. 1982, ApJS, 48, 95
657 *** Other ions of Na, Al -
658 *** >>refer Na, Al recom Landini, M., & Monsignori-Fossi, B.C. 1990, A&AS, 82, 229
659 *** Other ions of F, P, Cl, K, Ti, Cr, Mn, Co (excluding Ti I-II, Cr I-IV,
660 *** Mn I-V, Co I) -
661 *** >>refer many recom Landini, M., & Monsignori-Fossi, B.C. 1991, A&AS, 91, 183
662 *** All other species - interpolations of the power-law fits
663 *** Input parameters: iz - atomic number
664 *** in - number of electrons from 1 to iz
665 *** t - temperature, K
666 *** return result: - rate coefficient, cm^3 s^(-1)
667 ******************************************************************************
668 */
669 double tt;
670 double rate;
671
672 DEBUG_ENTRY( "rad_rec()" );
673
674 rate = 0.0;
675 if( iz < 1 || iz > 30 )
676 {
677 fprintf( ioQQQ, " rad_rec called with insane atomic number, =%4ld\n",
678 iz );
680 }
681 if( in < 1 || in > iz )
682 {
683 fprintf( ioQQQ, " rad_rec called with insane number elec =%4ld\n",
684 in );
686 }
687 if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) ||
688 (iz == 26 && in > 11) )
689 {
690 tt = sqrt(t/rnew[in-1][iz-1][2]);
691 rate =
692 rnew[in-1][iz-1][0]/(tt*pow(tt + 1.0,1.0 - rnew[in-1][iz-1][1])*
693 pow(1.0 + sqrt(t/rnew[in-1][iz-1][3]),1.0 + rnew[in-1][iz-1][1]));
694 }
695 else
696 {
697 tt = t*1.0e-04;
698 if( iz == 26 && in <= 13 )
699 {
700 rate = fe[in-1][0]/pow(tt,fe[in-1][1] +
701 fe[in-1][2]*log10(tt));
702 }
703 else
704 {
705 rate = rrec[in-1][iz-1][0]/pow(tt,(double)rrec[in-1][iz-1][1]);
706 }
707 }
708
709 return rate;
710}
711
712double t_ADfA::H_rad_rec(long int iz,
713 long int n,
714 double t)
715{
716 /*
717 * Version 4, October 9, 1997
718 ******************************************************************************
719 *** This subroutine calculates state-specific recombination rates
720 *** for hydrogen and hydrogen-like ions.
721 *** Input parameters: iz - atomic number
722 *** n - shell number, from 0 to 400:
723 *** 0 - 1s
724 *** 1 - 2s
725 *** 2 - 2p
726 *** 3 - 3
727 *** ......
728 *** t - temperature, K
729 *** Output parameter: r - rate coefficient, cm^3 s^(-1)
730 *** If n is negative, the subroutine returns the total recombination
731 *** rate coefficient
732 ******************************************************************************
733 */
734 double rate,
735 TeScaled,
736 x,
737 x1,
738 x2;
739
740 DEBUG_ENTRY( "H_rad_rec()" );
741
742 rate = 0.0;
743
744 /* iz is charge, must be 1 or greater */
745 ASSERT( iz > 0 );
746
747 /* n is level number, must be less than dim or hydro vectors */
749
750 TeScaled = t/POW2((double)iz);
751
752 if( n < 0 )
753 {
754 x1 = sqrt(TeScaled/3.148);
755 x2 = sqrt(TeScaled/7.036e05);
756 rate = 7.982e-11/x1/pow(1.0 + x1,0.252)/pow(1.0 + x2,1.748);
757 }
758 else
759 {
760 x = log10(TeScaled);
761 rate = (HRF[n][0] + HRF[n][2]*x + HRF[n][4]*
762 x*x + HRF[n][6]*powi(x,3) + HRF[n][8]*powi(x,4))/
763 (1.0 + HRF[n][1]*x + HRF[n][3]*x*x + HRF[n][5]*
764 powi(x,3) + HRF[n][7]*powi(x,4));
765 rate = pow(10.0,rate)/TeScaled;
766 }
767 rate *= iz;
768
769 return rate;
770}
771
772/*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
773 * returns collisional ionization rate coefficient cm^3 s^-1*/
775 /* atomic number, 1 for hydrogen */
776 long int iz,
777 /* number of bound electrons before ionization*/
778 long int in,
779 /* temperature */
780 double t)
781{
782 double rate, te, u;
783
784 DEBUG_ENTRY( "coll_ion()" );
785 /*D Verner's routine to compute collisional ionization rate coefficients
786 * Version 3, April 21, 1997
787 * Cu (Z=29) and Zn (Z=30) are added (fits from Ni, correct thresholds).
788 ******************************************************************************
789 *** This subroutine calculates rates of direct collisional ionization
790 *** for all ionization stages of all elements from H to Ni (Z=28)
791 *** by use of the fits from
792 *>>refer all coll_ion Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1
793 *** Input parameters: iz - atomic number on pphysical scale, H is 1
794 *** in - number of electrons from 1 to iz
795 *** t - temperature, K
796 *** Output parameter: c - rate coefficient, cm^3 s^(-1)
797 ****************************************************************************** */
798 rate = 0.0;
799
800 te = t*EVRYD/TE1RYD;
801 u = CF[in-1][iz-1][0]/te;
802 if( u > 80.0 )
803 {
804 return( 0. );
805 }
806
807 rate = (CF[in-1][iz-1][2]*(1.0 + CF[in-1][iz-1][1]*
808 sqrt(u))/(CF[in-1][iz-1][3] + u)*pow(u,CF[in-1][iz-1][4])*
809 exp(-u));
810
811 return(rate);
812}
813
815 /* (atomic number - 1), 0 for hydrogen */
816 long int z,
817 /* stage of ionization, 0 for atom */
818 long int n,
819 /* temperature K */
820 double t)
821{
822 double rate = 0.0;
823
824 DEBUG_ENTRY( "coll_ion_wrapper()" );
825
826 if( z < 0 || z > LIMELM-1 )
827 {
828 /* return zero rate is atomic number outside range of code */
829 return( 0. );
830 }
831
832 if( n < 0 || n > z )
833 {
834 /* return zero rate is ion stage is impossible */
835 return( 0. );
836 }
837
838 /*Get the rate coefficients from either Dima or Hybrid.
839 * The enum variable CIRCData is set to a default in init_defaults_preparse
840 * It can be changed parse_set via the input command set collisional ionization XXX */
841 if( atmdat.CIRCData == t_atmdat::DIMA)
842 {
843 /* collisional ionization by thermal electrons
844 * >>chng 97 mar 19, to Dima's new routine using
845 * >>refer all coll_ion Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1
846 *coll_ion(atomic number, number of electrons, temperature)*/
847 rate = coll_ion( z+1, z+1-n, t );
848 }
849 else if( atmdat.CIRCData == t_atmdat::HYBRID)
850 {
851 // Get the rate coefficient by scaling Dima to Dere.
852 rate = coll_ion_hybrid( z, n, t );
853 }
854 else
855 {
856 //CIRCData is an enum so it must be equal to one of the enum values.
858 }
859
860 ASSERT(rate >= 0.0);
861 return(rate);
862}
863
864/*coll_ion D Verner's routine to compute collisional ionization rate coefficients,
865 * returns collisional ionization rate coefficient cm^3 s^-1*/
867 /* (atomic number - 1) 0 for hydrogen */
868 long int nelem,
869 /* stage of ionization, 0 for atom */
870 long int ion,
871 /* temperature */
872 double t)
873{
874
875 DEBUG_ENTRY( "coll_ion_hybrid()" );
876
880 static const double DereRatio[30][30]=
881 {{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},
882 {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},
883 {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},
884 {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},
885 {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},
886 {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},
887 {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},
888 {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},
889 {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},
890 {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},
891 {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},
892 {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},
893 {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},
894 {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},
895 {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},
896 {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},
897 {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},
898 {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},
899 {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},
900 {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},
901 {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},
902 {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},
903 {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},
904 {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},
905 {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},
906 {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},
907 {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},
908 {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},
909 {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},
910 {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}};
911 /* Get the hybrid coeffs from the Dima rate coefficient times the average Dere to Dima ratio.
912 */
913 ASSERT( nelem>=0 && nelem<LIMELM && ion>=0 && ion<=nelem );
914 double rate = coll_ion(nelem+1,nelem+1-ion,t) * DereRatio[nelem][ion];
915 ASSERT( rate >=0. );
916 return(rate);
917}
918realnum t_ADfA::h_coll_str( long ipLo, long ipHi, long ipTe )
919{
920 realnum rate;
921
922 DEBUG_ENTRY( "h_coll_str()" );
923
924 ASSERT( ipLo < ipHi );
925
926# if !defined NDEBUG
927 long ipISO = ipH_LIKE;
928 long nelem = ipHYDROGEN;
929# endif
930 ASSERT( N_(ipLo) < N_(ipHi) );
931 ASSERT( N_(ipHi) <= 5 );
932
933 rate = (realnum)HCS[ipHi-1][ipLo][ipTe];
934
935 return rate;
936}
@ PHFIT95
Definition atmdat.h:276
@ PHFIT_UNDEF
Definition atmdat.h:276
t_atmdat atmdat
Definition atmdat.cpp:6
static double x2[63]
static double x1[83]
static double a1[83]
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int NHYDRO_MAX_LEVEL
Definition cddefines.h:266
const int LIMELM
Definition cddefines.h:258
#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
long nint(double x)
Definition cddefines.h:719
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double powi(double, long int)
Definition service.cpp:604
realnum fe[13][3]
Definition atmdat.h:299
long int NTOT[30]
Definition atmdat.h:288
realnum ST[9][405]
Definition atmdat.h:295
realnum PH2[30][30][7]
Definition atmdat.h:290
realnum PHH[NHYDRO_MAX_LEVEL][5]
Definition atmdat.h:292
double phfit(long int nz, long int ne, long int is, double e)
void rec_lines(double t, realnum r[][471])
realnum ph1(int i, int j, int k, int l) const
Definition atmdat.h:329
double coll_ion_wrapper(long int z, long int n, double t)
double HCS[14][10][8]
Definition atmdat.h:313
double rad_rec(long int iz, long int in, double t)
double coll_ion_hybrid(long int z, long int n, double t)
realnum STH[NHYDRO_MAX_LEVEL]
Definition atmdat.h:305
realnum h_coll_str(long ipLo, long ipHi, long ipTe)
long int L[7]
Definition atmdat.h:286
double CF[30][30][5]
Definition atmdat.h:307
double H_rad_rec(long int iz, long int n, double t)
phfit_version version
Definition atmdat.h:284
realnum P[8][110]
Definition atmdat.h:294
realnum HRF[NHYDRO_MAX_LEVEL][9]
Definition atmdat.h:301
realnum rnew[30][30][4]
Definition atmdat.h:298
double hpfit(long int iz, long int n, double e)
long int NINN[30]
Definition atmdat.h:287
double coll_ion(long int iz, long int in, double t)
realnum PH1[7][30][30][6]
Definition atmdat.h:289
realnum rrec[30][30][2]
Definition atmdat.h:297
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
#define N_(A_)
Definition iso.h:20
const int ipH_LIKE
Definition iso.h:62
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double TE1RYD
Definition physconst.h:183
static double * ex
Definition species2.cpp:28
@ HYBRID
Definition atmdat.h:259