00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "cddefines.h"
00018 #include "physconst.h"
00019 #include "hydro_bauman.h"
00020 #include "iso.h"
00021 #include "helike.h"
00022 #include "helike_recom.h"
00023 #include "thirdparty.h"
00024 #include "dense.h"
00025 #include "opacity.h"
00026 #include "atmdat.h"
00027 #include "taulines.h"
00028
00029
00030 STATIC double ExponentialInt( double v );
00031 STATIC double X1Int( double u );
00032 STATIC double X2Int( double u );
00033
00034 STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s);
00035 STATIC double GetHS98CrossSection( long n, long l, long s, double EgammaRyd );
00036
00037 static double Xn_S59;
00038
00039
00040
00041
00042
00043
00044 double He_cross_section( double EgammaRyd , double EthRyd, long n, long l, long S, long nelem )
00045 {
00046
00047 double cs = cross_section( EgammaRyd, EthRyd, nelem, n, l, S );
00048
00049
00050 if( nelem==ipHELIUM && n <=4 && l<=2 )
00051 {
00052 double rescaled[21] = {
00053 7.394,
00054 5.485, 9.219, 15.985, 15.985, 15.985, 13.504,
00055 8.018, 14.417, 28.501, 18.486, 18.132, 27.009,
00056 10.721, 20.235, 41.568, 36.717, 35.766, -1.000, -1.000, 41.787 };
00057 long ipLev = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[n][l][S];
00058 ASSERT( rescaled[ipLev] > 0. );
00059 cs *= rescaled[ipLev]/cross_section( EthRyd, EthRyd, nelem, n, l, S );
00060 }
00061
00062
00063 return cs * (1.e-18);
00064 }
00065
00066
00067
00068 STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long S)
00069 {
00070
00071
00072
00073 double E0[29] = {
00074 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,
00075 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,
00076 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};
00077 double sigma[29] = {
00078 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,
00079 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,
00080 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};
00081 double y_a[29] = {
00082 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,
00083 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,
00084 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};
00085 double P[29] = {
00086 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,
00087 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,
00088 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};
00089 double y_w[29] =
00090 {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};
00091 double yzero[29] =
00092 {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};
00093 double yone[29] =
00094 {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};
00095
00096 double pcs,Egamma,y,F,x;
00097 double rel_photon_energy;
00098
00099 Egamma = EgammaRyd * EVRYD;
00100
00101
00102
00103
00104 rel_photon_energy = EgammaRyd / EthRyd;
00105 rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );
00106
00107 long s=0;
00108 if( S == 1 )
00109 s=0;
00110 else if( S == 3 )
00111 s=1;
00112 else
00113 TotalInsanity();
00114
00115 if( nelem==ipHELIUM && n<=25 && l<=4 )
00116 {
00117
00118 pcs = GetHS98CrossSection( n, l, s, EgammaRyd );
00119 }
00120 else if( nelem==ipHELIUM && n>25 && l<=2 )
00121 {
00122 double scale[3][2] = {
00123 {1.4673,1.3671},
00124 {1.5458,1.5011},
00125 {1.4912,1.5144}};
00126
00127 long ipLev = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[25][l][S];
00128 double EthRyd_25 = iso_sp[ipHE_LIKE][nelem].fb[ipLev].xIsoLevNIonRyd;
00129 pcs = GetHS98CrossSection( 25, l, s, EthRyd_25 * rel_photon_energy );
00130 pcs *= pow((double)n/25., scale[l][s]);
00131 }
00132 else if( n==1 )
00133 {
00134
00135
00136
00137 x = Egamma/E0[nelem-1] - yzero[nelem-1];
00138 y = sqrt(x*x + yone[nelem-1]*yone[nelem-1]);
00139 F = ((x-1)*(x-1)+y_w[nelem-1]*y_w[nelem-1])
00140 * pow(y,0.5*P[nelem-1]-5.5) * pow((1+sqrt(y/y_a[nelem-1])),-P[nelem-1]);
00141 pcs = sigma[nelem-1]*F;
00142 }
00143 else if( nelem>=ipLITHIUM && nelem<=ipCALCIUM && n<11 && OP_Helike_NumPts[nelem][n][l][s]>0 )
00144 {
00145
00146 long numDataPoints = OP_Helike_NumPts[nelem][n][l][s];
00147 ASSERT( numDataPoints > 0 );
00148 ASSERT( OP_Helike_Xsectn[nelem][n][l][s] != NULL );
00149
00150 if( EgammaRyd < OP_Helike_Energy[nelem][n][l][s][numDataPoints-1] )
00151 {
00152 pcs = linint( OP_Helike_Energy[nelem][n][l][s], OP_Helike_Xsectn[nelem][n][l][s], numDataPoints, EgammaRyd );
00153 }
00154 else
00155 {
00156
00157 pcs = OP_Helike_Xsectn[nelem][n][l][s][numDataPoints-1] * POW3( OP_Helike_Energy[nelem][n][l][s][numDataPoints-1]/EgammaRyd );
00158 }
00159 }
00160 else
00161 {
00162
00163 pcs = (1.e18)*H_photo_cs(rel_photon_energy , n, l, nelem);
00164 }
00165
00166 ASSERT( pcs > 0. && pcs < 1.E10 );
00167
00168 return pcs;
00169 }
00170
00171 STATIC double GetHS98CrossSection( long n, long l, long s, double EgammaRyd )
00172 {
00173 double pcs;
00174 ASSERT( n<=25 );
00175 ASSERT( l<=4 );
00176 ASSERT( s==0 || s==1 );
00177
00178
00179 if( EgammaRyd < HS_He1_Energy[n][l][s][NUM_HS98_DATA_POINTS-1] )
00180 {
00181 pcs = linint( HS_He1_Energy[n][l][s], HS_He1_Xsectn[n][l][s], NUM_HS98_DATA_POINTS, EgammaRyd );
00182 }
00183 else
00184 {
00185
00186 pcs = HS_He1_Xsectn[n][l][s][NUM_HS98_DATA_POINTS-1] * POW3( HS_He1_Energy[n][l][s][NUM_HS98_DATA_POINTS-1]/EgammaRyd );
00187 }
00188
00189 return pcs;
00190 }
00191
00192
00193 double Recomb_Seaton59( long nelem, double temp, long n)
00194 {
00195 double lambda = TE1RYD * nelem * nelem / temp;
00196
00197
00198 double x = lambda / n / n;
00199 double AlphaN;
00200 double SzeroOfX = 0.;
00201 double SoneOfX = 0.;
00202 double StwoOfX = 0.;
00203 double SnOfLambda = 0.;
00204 double lowerlimit, upperlimit, step;
00205
00206 Xn_S59 = x;
00207
00208
00209 lowerlimit = x;
00210 step = 3. * x;
00211 upperlimit = lowerlimit + step;
00212 SzeroOfX = qg32( lowerlimit, upperlimit, ExponentialInt);
00213
00214 do
00215 {
00216 lowerlimit = upperlimit;
00217 step *= 2;
00218 upperlimit = lowerlimit + step;
00219 SzeroOfX += qg32( lowerlimit, upperlimit, ExponentialInt);
00220 } while ( upperlimit < 20. );
00221
00222
00223
00224
00225
00226
00227 lowerlimit = 0.;
00228 step = 0.5;
00229 upperlimit = lowerlimit + step;
00230 SoneOfX = qg32( lowerlimit, upperlimit, X1Int);
00231 StwoOfX = qg32( lowerlimit, upperlimit, X2Int);
00232
00233 do
00234 {
00235 lowerlimit = upperlimit;
00236 step *= 2;
00237 upperlimit = lowerlimit + step;
00238 SoneOfX += qg32( lowerlimit, upperlimit, X1Int);
00239 StwoOfX += qg32( lowerlimit, upperlimit, X2Int);
00240 } while ( upperlimit < 200. );
00241
00242 SoneOfX *= 0.1728 * pow( x, 1./3. );
00243 StwoOfX *= -0.0496 * pow( x, 2./3. );
00244
00245
00246 SnOfLambda = SzeroOfX + pow(1./lambda, 1./3.)*SoneOfX + pow(1./lambda, 2./3.)*StwoOfX;
00247
00248 AlphaN = 5.197E-14 * nelem * pow(x, 1.5) * SnOfLambda;
00249
00250 return AlphaN;
00251
00252 }
00253
00254 STATIC double ExponentialInt( double v )
00255 {
00256 double Integrand;
00257
00258 Integrand = exp( -1. * v + Xn_S59) / v;
00259
00260 return Integrand;
00261 }
00262
00263 STATIC double X1Int( double u )
00264 {
00265 double Integrand;
00266
00267 Integrand = pow(1./(u + 1.), 5./3.) * (u - 1.) * exp( -1. * Xn_S59 * u );
00268
00269 return Integrand;
00270 }
00271
00272 STATIC double X2Int( double u )
00273 {
00274 double Integrand;
00275
00276 Integrand = pow(1./(u + 1.), 7./3.) * (u*u + 4./3.*u + 1.) * exp( -1. * Xn_S59 * u );
00277
00278 return Integrand;
00279 }
00280