44double He_cross_section(
double EgammaRyd ,
double EthRyd,
long n,
long l,
long S,
long nelem )
50 if( nelem==
ipHELIUM && n <=5 && l<=2 )
52 double rescaled[31] = {
54 5.485, 9.219, 15.985, 15.985, 15.985, 13.504,
55 8.018, 14.417, 28.501, 18.486, 18.132, 27.009,
56 10.721, 20.235, 41.568, 36.717, 35.766, -1.000, -1.000, 41.787,
57 13.527, 26.539, 55.692, 55.010, 53.514, -1.000, -1.000, -1.000, -1.000, 58.120 };
59 ASSERT( rescaled[ipLev] > 0. );
60 cs *= rescaled[ipLev]/
cross_section( EthRyd, EthRyd, nelem, n, l,
S );
75 1.36E+01,2.01E+01,1.76E+01,3.34E+01,4.62E+01,6.94E+01,8.71E+01,1.13E+02,1.59E+02,2.27E+02,
76 2.04E+02,2.74E+02,2.75E+02,3.38E+02,4.39E+02,4.17E+02,4.47E+02,5.18E+02,6.30E+02,6.27E+02,
77 8.66E+02,7.67E+02,9.70E+02,9.66E+02,1.06E+03,1.25E+03,1.35E+03,1.43E+03,1.56E+03};
79 9.49E+02,3.20E+02,5.46E+02,2.85E+02,2.34E+02,1.52E+02,1.33E+02,1.04E+02,6.70E+01,4.00E+01,
80 6.14E+01,4.04E+01,4.75E+01,3.65E+01,2.45E+01,3.14E+01,3.11E+01,2.59E+01,1.94E+01,2.18E+01,
81 1.23E+01,1.76E+01,1.19E+01,1.31E+01,1.20E+01,9.05E+00,8.38E+00,8.06E+00,7.17E+00};
83 1.47E+00,7.39E+00,1.72E+01,2.16E+01,2.18E+01,2.63E+01,2.54E+01,2.66E+01,3.35E+01,5.32E+01,
84 2.78E+01,3.57E+01,2.85E+01,3.25E+01,4.41E+01,3.16E+01,3.04E+01,3.28E+01,3.92E+01,3.45E+01,
85 5.89E+01,3.88E+01,5.35E+01,4.83E+01,5.77E+01,6.79E+01,7.43E+01,7.91E+01,9.10E+01};
87 3.19E+00,2.92E+00,3.16E+00,2.62E+00,2.58E+00,2.32E+00,2.34E+00,2.26E+00,2.00E+00,1.68E+00,
88 2.16E+00,1.92E+00,2.14E+00,2.00E+00,1.77E+00,2.04E+00,2.09E+00,2.02E+00,1.86E+00,2.00E+00,
89 1.62E+00,1.93E+00,1.70E+00,1.79E+00,1.72E+00,1.61E+00,1.59E+00,1.58E+00,1.54E+00};
91 {2.039,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};
93 {0.4434,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};
95 {2.136,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};
97 double pcs,Egamma,y,F,x;
98 double rel_photon_energy;
100 Egamma = EgammaRyd *
EVRYD;
105 rel_photon_energy = EgammaRyd / EthRyd;
106 rel_photon_energy =
MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );
116 if( nelem==
ipHELIUM && n<=25 && l<=4 )
121 else if( nelem==
ipHELIUM && n>25 && l<=2 )
123 double scale[3][2] = {
131 pcs *= pow((
double)n/25., scale[l][s]);
138 x = Egamma/E0[nelem-1] - yzero[nelem-1];
139 y = sqrt(x*x + yone[nelem-1]*yone[nelem-1]);
140 F = ((x-1)*(x-1)+y_w[nelem-1]*y_w[nelem-1])
141 * pow(y,0.5*P[nelem-1]-5.5) * pow((1+sqrt(y/y_a[nelem-1])),-P[nelem-1]);
142 pcs = sigma[nelem-1]*F;
148 ASSERT( numDataPoints > 0 );
164 pcs = (1.e18)*
H_photo_cs(rel_photon_energy , n, l, nelem);
167 ASSERT( pcs > 0. && pcs < 1.E10 );
196 double lambda =
TE1RYD * nelem * nelem / temp;
199 double x = lambda / n / n;
201 double SzeroOfX = 0.;
204 double SnOfLambda = 0.;
205 double lowerlimit, upperlimit, step;
212 upperlimit = lowerlimit + step;
217 lowerlimit = upperlimit;
219 upperlimit = lowerlimit + step;
221 }
while ( upperlimit < 20. );
230 upperlimit = lowerlimit + step;
231 SoneOfX =
qg32( lowerlimit, upperlimit,
X1Int);
232 StwoOfX =
qg32( lowerlimit, upperlimit,
X2Int);
236 lowerlimit = upperlimit;
238 upperlimit = lowerlimit + step;
239 SoneOfX +=
qg32( lowerlimit, upperlimit,
X1Int);
240 StwoOfX +=
qg32( lowerlimit, upperlimit,
X2Int);
241 }
while ( upperlimit < 200. );
243 SoneOfX *= 0.1728 * pow( x, 1./3. );
244 StwoOfX *= -0.0496 * pow( x, 2./3. );
247 SnOfLambda = SzeroOfX + pow(1./lambda, 1./3.)*SoneOfX + pow(1./lambda, 2./3.)*StwoOfX;
249 AlphaN = 5.197E-14 * nelem * pow(x, 1.5) * SnOfLambda;