46 long i, i1, j, nelem, ipHi, ipLo;
51 # define chLine_LENGTH 1000
58 fprintf(
ioQQQ,
" HeCollidSetup opening he1_cs.dat:");
65 fprintf(
ioQQQ,
" HeCollidSetup could not read first line of he1_cs.dat.\n");
77 " HeCollidSetup: the version of he1_cs.dat is not the current version.\n" );
79 " HeCollidSetup: I expected to find the number %i and got %li instead.\n" ,
81 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
87 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
107 fprintf(
ioQQQ,
" HeCollidSetup could not find line in CS temperatures in he1_cs.dat.\n");
132 for( ipLo=
ipHe1s1S; ipLo<ipHi; ++ipLo )
144 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
148 if( (chLine[0] ==
' ') || (chLine[0]==
'\n') )
150 if( chLine[0] !=
'#')
167 if( (chTemp = strchr( chTemp,
'\t' )) == NULL )
169 fprintf(
ioQQQ,
" HeCollidSetup no init cs\n" );
175 for( i=0; i<
iso.
nCS[ipISO]; ++i )
178 if( (chTemp = strchr( chTemp,
'\t' )) == NULL )
180 fprintf(
ioQQQ,
" HeCollidSetup not scan cs, current indices: %li %li\n", ipHi, ipLo );
184 sscanf( chTemp ,
"%le" , &a );
207 const char *where =
" ";
219 cs =
AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
224 cs =
IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
249 enum {DEBUG_LOC=
false};
254 fprintf(
ioQQQ,
"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n",
260 return MAX2(cs, 1.e-10f);
333 ASSERT( *factor1 == -1.f );
343 *factor1 = (2.f*((
realnum)ipLo-3.f)+1.f) / 9.f;
353 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
387 ASSERT( (flow >= 0.f) && (flow <= 1.f) );
412 ASSERT( *factor1 == -1.f );
464 *factor1 = (2.f*((
realnum)ipLo-3.f)+1.f) / 9.f;
470 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
480 ASSERT( *factor1 == -1.f );
528 cs = (
realnum)pow( 10. , -1.45*x + 6.75);
533 cs = (
realnum)pow( 10. , -3.33*x+15.15);
539 cs = (
realnum)pow( 10. , -2.3*x+10.3);
566 *factor1 = (2.f*((
realnum)ipLo-3.f)+1.f) / 9.f;
572 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
634 cs = 0.25f/
POW2(nelem+1.f);
637 cs = 0.4f/
POW2(nelem+1.f);
640 cs = 0.15f/
POW2(nelem+1.f);
643 cs = 0.45f/
POW2(nelem+1.f);
646 cs = 0.75f/
POW2(nelem+1.f);
649 cs = 1.3f/
POW2(nelem+1.f);
663 cs = 2.75f/
POW2(nelem+1.f);
666 cs = 60.f/
POW2(nelem+1.f);
669 cs = 180.f/
POW2(nelem+1.f);
672 cs = 300.f/
POW2(nelem+1.f);
675 cs = 5.8f/
POW2(nelem+1.f);
689 cs = 0.56f/
POW2(nelem+1.f);
692 cs = 1.74f/
POW2(nelem+1.f);
695 cs = 2.81f/
POW2(nelem+1.f);
698 cs = 190.f/
POW2(nelem+1.f);
712 cs = 8.1f/
POW2(nelem+1.f);
715 cs = 8.2f/
POW2(nelem+1.f);
718 cs = 3.9f/
POW2(nelem+1.f);
732 cs = 30.f/
POW2(nelem+1.f);
735 cs = 11.7f/
POW2(nelem+1.f);
759 ASSERT( *factor1 == -1.f );
810 *factor1 = (2.f*((
realnum)ipHi-3.f)+1.f) / 9.f;
877 double coll_str, deltaE;
881 if(
Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
915 double integrand,
cross_section, osc_strength, coll_str, zOverB2;
916 double deltaE, betaone, zeta_of_betaone, cs2, temp, stat_weight;
917 double I_energy_eV, Dubya, proj_energy;
937 proj_energy += deltaE;
938 Dubya = 0.5*(2.*proj_energy-deltaE);
943 ASSERT( I_energy_eV > 0. );
944 ASSERT( osc_strength > 0. );
947 zOverB2 = 0.5*
POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
953 betaone = sqrt( 1./zOverB2 );
955 else if( zOverB2 < 0.54 )
958 betaone = (1./3.)*(log(
PI)-log(zOverB2)+1.28);
962 betaone += 0.5*(log(
PI)-log(zOverB2));
969 double zetaOVERbeta2[101] = {
970 1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
971 7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
972 4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
973 3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
974 2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
975 1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
976 1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
977 7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
978 4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
979 3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
980 2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
981 1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
982 7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
984 ASSERT( zOverB2 >= zetaOVERbeta2[100] );
987 for( i=0; i< 100; ++i )
989 if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
997 ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
999 betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) /
1000 (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
1001 (pow(10., ((
double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((
double)ip_zOverB2)/100. - 1.)) +
1002 pow(10., (
double)ip_zOverB2/100. - 1.);
1006 zeta_of_betaone = zOverB2 *
POW2(betaone);
1012 cross_section = cs2;
1015 cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
1018 coll_str = cross_section * stat_weight * ( proj_energy/
EVRYD );
1020 integrand = exp( -1.*(proj_energy-deltaE)*
EVDEGK/temp ) * coll_str;
1035 TwoLogDebye, TwoLogRc1,
1037 bestfactor, factorpart,
1038 reduced_mass, reduced_mass_2_emass,
1055 tau =
StatesElem[ipISO][nelem][ipLo].lifetime;
1073 TwoLogRc1 = 10.95 + log10(
phycon.
te * tau * tau / reduced_mass_2_emass );
1076 Dnl =
POW2( ChargIncoming / (
double)(nelem + 1. - ipISO)) * 6. *
POW2( (
double)
StatesElem[ipISO][nelem][ipLo].n ) *
1083 factorpart = (11.54 + log10(
phycon.
te / Dnl / reduced_mass_2_emass ) );
1085 if( (factor1 = factorpart + TwoLogDebye) <= 0.)
1087 if( (factor2 = factorpart + TwoLogRc1) <= 0.)
1091 bestfactor =
MIN2(factor1,factor2);
1093 ASSERT( bestfactor > 0. );
1096 if( bestfactor > 100. )
1102 rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnl /
phycon.
sqrte * bestfactor;
1118 cs = rate / (
COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) *
1207 double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
1208 double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
1210 double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
1211 double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 ;
1213 long ipISO, nelem, n, l, lp, s, ipLo, ipHi, Collider =
global_Collider;
1229 stat_weight =
StatesElem[ipISO][nelem][ipLo].g;
1242 ASSERT( reduced_mass > 0. );
1247 ASSERT( aveRadius < 1.e-4 );
1261 if( lgParamIsRedVel )
1264 reduced_vel = velOrEner;
1273 E_Proj_Ryd = velOrEner;
1278 ASSERT( reduced_vel > 1.e-10 );
1279 ASSERT( reduced_vel < 1.e10 );
1284 ConstantFactors = 4.5*
PI*
POW2(ColliderCharge*aveRadius/reduced_vel);
1287 reduced_b_min = 1.5 * ColliderCharge / reduced_vel;
1288 ASSERT( reduced_b_min > 1.e-10 );
1289 ASSERT( reduced_b_min < 1.e10 );
1299 quantum_defect1 = (double)n- (
double)nelem/sqrt(
iso.
xIsoLevNIonRyd[ipISO][nelem][ipLo]);
1300 quantum_defect2 = (double)n- (
double)nelem/sqrt(
iso.
xIsoLevNIonRyd[ipISO][nelem][ipHi]);
1303 ASSERT( fabs(quantum_defect1) < 1.0 );
1304 ASSERT( fabs(quantum_defect1) > 0.0 );
1305 ASSERT( fabs(quantum_defect2) < 1.0 );
1306 ASSERT( fabs(quantum_defect2) > 0.0 );
1309 omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*
POW2((
double)l/(
double)n)) /
POW3( (
double)n ) / (
double)l );
1311 omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*
POW2((
double)lp/(
double)n)) /
POW3( (
double)n ) / (
double)lp );
1313 omega_qd = 0.5*( omega_qd1 + omega_qd2 );
1317 reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
1323 reduced_b_max =
MAX2( reduced_b_max, reduced_b_min );
1324 ASSERT( reduced_b_max > 0. );
1326 alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max);
1327 alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min);
1332 alphamin =
MAX2( alphamin, 1.e-30 );
1333 alphamax =
MAX2( alphamax, 1.e-20 );
1337 if( alphamax > alphamin )
1340 step = (alphamax - alphamin)/5.;
1347 cross_section = ConstantFactors * CSIntegral;
1350 coll_str = cross_section * stat_weight /
PI / BOHR_RADIUS_CM /
BOHR_RADIUS_CM;
1351 coll_str *= E_Proj_Ryd;
1360 double integrand, deltaPhi, b;
1364 ASSERT( alpha >= 1.e-30 );
1379 integrand = 1./alpha/alpha/alpha;
1388 double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
1395 cosU1 = 2.*
POW2((
double)l/(
double)n) - 1.;
1396 cosU2 = 2.*
POW2((
double)lp/(
double)n) - 1.;
1398 sinU1 = sqrt( 1. - cosU1*cosU1 );
1399 sinU2 = sqrt( 1. - cosU2*cosU2 );
1402 cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
1403 sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
1404 cosChi = 2. *
POW2( cosChiOver2 ) - 1.;
1408 if( -1.*cosU2 - cosChi < 0. )
1415 ASSERT( sinChiOver2 > 0. );
1416 ASSERT( sinChiOver2*sinChiOver2 >
POW2((
double)lp/(
double)n) );
1419 probability = (double)lp/(
POW2((
double)n)*sinChiOver2*sqrt(
POW2(sinChiOver2) -
POW2((
double)lp/(
double)n) ) );
1424 double OneMinusCosChi = 1. - cosChi;
1426 if( OneMinusCosChi == 0. )
1428 double hold = sin( deltaPhi / 2. );
1430 OneMinusCosChi = 8. * alpha * alpha *
POW2( hold );
1433 if( OneMinusCosChi == 0. )
1440 A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
1441 B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
1454 probability = 2.*lp/(
PI* n*n*
POW2( sinChiOver2 ));
1458 probability *=
ellpk( -A/(B-A) );
1459 probability /= sqrt( B - A );
1463 probability *=
ellpk( A/B );
1464 probability /= sqrt( B );
1476 STATIC void updateHeColl(
long int nelem)
1478 long int ipLo , ipHi, i;
1483 static double TeUsedForCS[
LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1484 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0};
1496 if( (TeUsedForCS[nelem] == 0.) ||
1497 ( TeUsedForCS[nelem]/
phycon.
te > 1.15 ) ||
1498 ( TeUsedForCS[nelem]/
phycon.
te < 0.85 ) )
1502 for( ipLo=
ipHe1s1S; ipLo < ipHi; ipLo++ )