359 ASSERT( *factor1 == -1.f );
369 *factor1 = (2.f*((
realnum)ipLo-3.f)+1.f) / 9.f;
379 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
396 cs =
HeCS[ipHi][ipLo][0];
413 ASSERT( (flow >= 0.f) && (flow <= 1.f) );
415 cs =
HeCS[ipHi][ipLo][ipArray] * (1.f-flow) +
416 HeCS[ipHi][ipLo][ipArray+1] * flow;
438 ASSERT( *factor1 == -1.f );
497 *factor1 = (2.f*((
realnum)ipLo-3.f)+1.f) / 9.f;
503 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
513 ASSERT( *factor1 == -1.f );
561 cs = (
realnum)pow( 10. , -1.45*x + 6.75);
566 cs = (
realnum)pow( 10. , -3.33*x+15.15);
572 cs = (
realnum)pow( 10. , -2.3*x+10.3);
599 *factor1 = (2.f*((
realnum)ipLo-3.f)+1.f) / 9.f;
605 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
950 double stat_weight,
double I_energy_eV )
954 double betaone, zeta_of_betaone, cs2;
955 double Dubya, proj_energy;
971 proj_energy += deltaE;
972 Dubya = 0.5*(2.*proj_energy-deltaE);
977 ASSERT( I_energy_eV > 0. );
978 ASSERT( osc_strength > 0. );
981 zOverB2 = 0.5*
POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
987 betaone = sqrt( 1./zOverB2 );
989 else if( zOverB2 < 0.54 )
992 betaone = (1./3.)*(log(
PI)-log(zOverB2)+1.28);
996 betaone += 0.5*(log(
PI)-log(zOverB2));
1002 long ip_zOverB2 = 0;
1003 double zetaOVERbeta2[101] = {
1004 1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
1005 7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
1006 4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
1007 3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
1008 2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
1009 1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
1010 1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
1011 7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
1012 4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
1013 3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
1014 2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
1015 1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
1016 7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
1018 ASSERT( zOverB2 >= zetaOVERbeta2[100] );
1021 for(
long i=0; i< 100; ++i )
1023 if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
1031 ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
1033 betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) /
1034 (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
1035 (pow(10., ((
double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((
double)ip_zOverB2)/100. - 1.)) +
1036 pow(10., (
double)ip_zOverB2/100. - 1.);
1040 zeta_of_betaone = zOverB2 *
POW2(betaone);
1049 cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
1054 integrand = exp( -1.*(proj_energy-deltaE)*
EVDEGK/temp ) * coll_str;
1062 double target_charge,
1071 TwoLogDebye, TwoLogRc1,
1073 bestfactor, factorpart,
1074 reduced_mass, reduced_mass_2_emass,
1100 TwoLogRc1 = 10.95 + log10(
phycon.te * tau * tau / reduced_mass_2_emass );
1106 TwoLogRc2 = 2. * log10( 1.12 *
H_BAR * v / deltaE );
1112 Dnl =
POW2( ChargIncoming / target_charge) * 6. *
POW2( (
double)n) *
1113 (
POW2((
double)n) -
POW2((
double)l) - l - 1);
1118 factorpart = (11.54 + log10(
phycon.te / Dnl / reduced_mass_2_emass ) );
1120 if( (factor1 = factorpart + TwoLogDebye) <= 0.)
1122 if( (factor2 = factorpart + TwoLogRc1) <= 0.)
1126 bestfactor =
MIN2(factor1,factor2);
1128 ASSERT( bestfactor > 0. );
1131 if( bestfactor > 100. )
1137 rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnl /
phycon.sqrte * bestfactor;
1153 cs = rate / (
COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) *
1229STATIC double collision_strength_VF01(
long ipISO,
long nelem,
long n,
long l,
long lp,
long s,
long Collider,
double ColliderCharge,
double temp,
double velOrEner,
bool lgParamIsRedVel )
1231 double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
1232 double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
1233 double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
1234 double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 ;
1242 ipLo =
iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][s];
1243 ipHi =
iso_sp[ipISO][nelem].QuantumNumbers2Index[n][lp][s];
1244 stat_weight =
iso_sp[ipISO][nelem].st[ipLo].g();
1257 ASSERT( reduced_mass > 0. );
1262 ASSERT( aveRadius < 1.e-4 );
1275 if( lgParamIsRedVel )
1278 reduced_vel = velOrEner;
1287 E_Proj_Ryd = velOrEner;
1292 ASSERT( reduced_vel > 1.e-10 );
1293 ASSERT( reduced_vel < 1.e10 );
1300 ASSERT( reduced_b_min > 1.e-10 );
1301 ASSERT( reduced_b_min < 1.e10 );
1311 quantum_defect1 = (double)n- (
double)nelem/sqrt(
iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd);
1312 quantum_defect2 = (double)n- (
double)nelem/sqrt(
iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
1315 ASSERT( fabs(quantum_defect1) < 1.0 );
1316 ASSERT( fabs(quantum_defect1) > 0.0 );
1317 ASSERT( fabs(quantum_defect2) < 1.0 );
1318 ASSERT( fabs(quantum_defect2) > 0.0 );
1321 omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*
POW2((
double)l/(
double)n)) /
POW3( (
double)n ) / (
double)l );
1323 omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*
POW2((
double)lp/(
double)n)) /
POW3( (
double)n ) / (
double)lp );
1325 omega_qd = 0.5*( omega_qd1 + omega_qd2 );
1329 reduced_b_max = sqrt( 1.5 *
ColliderCharge * n / omega_qd )/aveRadius;
1335 reduced_b_max =
MAX2( reduced_b_max, reduced_b_min );
1336 ASSERT( reduced_b_max > 0. );
1344 func.
an = aveRadius;
1346 func.
bmax = reduced_b_max*aveRadius;
1355 alphamin =
MAX2( alphamin, 1.e-30 );
1356 alphamax =
MAX2( alphamax, 1.e-20 );
1360 if( alphamax > alphamin )
1363 step = (alphamax - alphamin)/5.;
1365 CSIntegral += VF01_alpha.
sum( alpha1, alpha1+step, func );
1366 CSIntegral += VF01_alpha.
sum( alpha1+step, alpha1+4.*step, func );
1410 double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
1417 cosU1 = 2.*
POW2((
double)l/(
double)n) - 1.;
1418 cosU2 = 2.*
POW2((
double)lp/(
double)n) - 1.;
1420 sinU1 = sqrt( 1. - cosU1*cosU1 );
1421 sinU2 = sqrt( 1. - cosU2*cosU2 );
1424 cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
1425 sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
1426 cosChi = 2. *
POW2( cosChiOver2 ) - 1.;
1430 if( -1.*cosU2 - cosChi < 0. )
1437 ASSERT( sinChiOver2 > 0. );
1438 ASSERT( sinChiOver2*sinChiOver2 >
POW2((
double)lp/(
double)n) );
1441 probability = (double)lp/(
POW2((
double)n)*sinChiOver2*sqrt(
POW2(sinChiOver2) -
POW2((
double)lp/(
double)n) ) );
1446 double OneMinusCosChi = 1. - cosChi;
1448 if( OneMinusCosChi == 0. )
1450 double hold = sin( deltaPhi / 2. );
1452 OneMinusCosChi = 8. * alpha * alpha *
POW2( hold );
1455 if( OneMinusCosChi == 0. )
1462 A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
1463 B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
1476 probability = 2.*lp/(
PI* n*n*
POW2( sinChiOver2 ));
1480 probability *=
ellpk( -A/(B-A) );
1481 probability /= sqrt( B - A );
1485 probability *=
ellpk( A/B );
1486 probability /= sqrt( B );