18 const long VERSION_MAGIC = 20061204
L;
20 const static char chFile[] =
"phfit.dat";
25 long i=-1, j=-1, k=-1, n;
27 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
28 if( lgErr || i != VERSION_MAGIC )
30 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile, i );
31 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
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 );
43 lgErr = lgErr || ( fscanf( io,
"%ld %ld %ld", &i, &j, &k ) != 3 );
44 if( i == -1 && j == -1 && k == -1 )
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 )
53 lgErr = lgErr || ( fscanf( io,
"%ld %ld", &i, &j ) != 2 );
54 if( i == -1 && j == -1 )
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 )
65 const static char chFile2[] =
"hpfit.dat";
69 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
70 if( lgErr || i != VERSION_MAGIC )
72 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile2, i );
73 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
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 )
89 const static char chFile3[] =
"rec_lines.dat";
93 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
94 if( lgErr || i != VERSION_MAGIC )
96 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile3, i );
97 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
101 for( i=0; i < 110; i++ )
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 )
110 for( i=0; i < 405; i++ )
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 )
122 const static char chFile4[] =
"rad_rec.dat";
126 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
127 if( lgErr || i != VERSION_MAGIC )
129 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile4, i );
130 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
136 lgErr = lgErr || ( fscanf( io,
"%ld %ld", &i, &j ) != 2 );
137 if( i == -1 && j == -1 )
139 lgErr = lgErr || ( fscanf( io,
"%le %le", &help[0], &help[1] ) != 2 );
140 for(
int l=0; l < 2; ++l )
145 lgErr = lgErr || ( fscanf( io,
"%ld %ld", &i, &j ) != 2 );
146 if( i == -1 && j == -1 )
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 )
155 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
158 lgErr = lgErr || ( fscanf( io,
"%le %le %le", &help[0], &help[1], &help[2] ) != 3 );
159 for(
int l=0; l < 3; ++l )
167 const static char chFile5[] =
"h_rad_rec.dat";
171 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
172 if( lgErr || i != VERSION_MAGIC )
174 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile5, i );
175 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
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 )
191 const static char chFile6[] =
"h_phot_cs.dat";
195 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
196 if( lgErr || i != VERSION_MAGIC )
198 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile6, i );
199 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
205 lgErr = lgErr || ( fscanf( io,
"%le", &help[0] ) != 1 );
213 const static char chFile7[] =
"coll_ion.dat";
217 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
218 if( lgErr || i != VERSION_MAGIC )
220 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile7, i );
221 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
227 lgErr = lgErr || ( fscanf( io,
"%ld %ld", &i, &j ) != 2 );
228 if( i == -1 && j == -1 )
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 );
238 const static char chFile8[] =
"h_coll_str.dat";
242 lgErr = lgErr || ( fscanf( io,
"%ld", &i ) != 1 );
243 if( lgErr || i != VERSION_MAGIC )
245 fprintf(
ioQQQ,
" File %s has incorrect version: %ld\n", chFile8, i );
246 fprintf(
ioQQQ,
" I expected to find version: %ld\n", VERSION_MAGIC );
252 lgErr = lgErr || ( fscanf( io,
"%ld %ld", &i, &j ) != 2 );
253 if( i == -1 && j == -1 )
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 );
312 if( nz < 1 || nz > 30 )
317 if( ne < 1 || ne > nz )
323 if( nz == ne && nz > 18 )
325 if( nz == (ne + 1) && ((((nz == 20 || nz == 21) || nz == 22) ||
326 nz == 25) || nz == 26) )
333 if( (is == 6 && (nz == 20 || nz == 19)) && ne >= 19 )
338 ASSERT( is >= 1 && is <= 7 );
340 if( e <
PH1[is-1][ne-1][nz-1][0] )
346 if( ((nz == 15 || nz == 17) || nz == 19) || (nz > 20 && nz != 26) )
358 einn =
PH1[nint-1][ne-1][nz-1][0];
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);
374 if( (is < nout && is > nint) && e < einn )
378 p1 = -
PH2[ne-1][nz-1][3];
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);
396 double crs , ex , eth;
414 eth =
ph1(0,0,iz-1,0)/
POW2((
double)m);
415 ex =
MAX2(1. , e/eth );
417 crs =
hpfit( iz , n , ex );
459 fprintf(
ioQQQ,
" hpfit called with too large n, =%li\n" , n );
468 q = 3.5 + l - 0.5*
PHH[n][1];
486 eth =
ph1(0,0,iz-1,0)/
POW2((
double)m);
487 ex =
MAX2(1. , e/eth );
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/
526 static long jd[6]={143,145,157,360,376,379};
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,
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};
553 for( i=0; i < 110; i++ )
558 z =
P[0][i] -
P[1][i] + 1.0;
566 a1 = p1*pow(0.004,p2)/(1.0 + p3*pow(0.004,p4));
567 a = a1/sqrt(te/0.004);
573 a1 = p1*pow(2.0,p2)/(1.0 + p3*pow(2.0,p4));
574 a = a1/pow(te/2.0,1.5);
578 a = p1*pow(te,p2)/(1.0 + p3*pow(te,p4));
581 rr[3][i] = 1.0e-13*z*a*
P[7][i];
584 for( i=0; i < 405; i++ )
598 x = p5*(1.0/te - 1.0/p6);
605 a1 = (p1/p6 + p2 + p3*p6 + p4*p6*p6)/pow(p6,1.5)/exp(p5/
614 a1 = (p1/6.0 + p2 + p3*6.0 + p4*36.0)/pow(6.0,1.5)/
616 a = a1/pow(te/6.0,1.5);
620 a = (p1/te + p2 + p3*te + p4*te*te)/pow(te,1.5)/exp(p5/
624 dr[3][i] = 1.0e-12*a;
627 for( i=0; i < 6; i++ )
631 dr[3][ipj1-1] += dr[3][ipj2-1];
635 for( i=0; i < 38; i++ )
639 rr[3][ipj1-1] += dr[3][ipj2-1];
643 for( i=0; i < 110; i++ )
652 for( i=0; i < 405; i++ )
704 if( iz < 1 || iz > 30 )
706 fprintf(
ioQQQ,
" rad_rec called with insane atomic number, =%4ld\n",
710 if( in < 1 || in > iz )
712 fprintf(
ioQQQ,
" rad_rec called with insane number elec =%4ld\n",
716 if( (((in <= 3 || in == 11) || (iz > 5 && iz < 9)) || iz == 10) ||
717 (iz == 26 && in > 11) )
719 tt = sqrt(t/
rnew[in-1][iz-1][2]);
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]));
727 if( iz == 26 && in <= 13 )
729 rate =
fe[in-1][0]/pow(tt,
fe[in-1][1] +
730 fe[in-1][2]*log10(tt));
734 rate =
rrec[in-1][iz-1][0]/pow(tt,(
double)
rrec[in-1][iz-1][1]);
779 TeScaled = t/
POW2((
double)iz);
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);
790 rate = (
HRF[n][0] +
HRF[n][2]*x +
HRF[n][4]*
792 (1.0 +
HRF[n][1]*x +
HRF[n][3]*x*x +
HRF[n][5]*
794 rate = pow(10.0,rate)/TeScaled;
829 if( iz < 1 || iz > 30 )
835 if( in < 1 || in > iz )
842 u =
CF[in-1][iz-1][0]/te;
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])*