24 STATIC double Hydcs123(
long int ilow,
long int ihigh,
long int iz,
long int chType);
37 static const realnum HCSTE[
NHCSTE] = {5802.f,11604.f,34812.f,58020.f,116040.f,174060.f,232080.f,290100.f};
48 if( ipLo==1 && ipHi==2 )
50 fprintf(
ioQQQ,
"HCSAR_interp was called for the 2s-2p transition, which it cannot do\n");
67 for( ip=1; ip<
NHCSTE; ++ip )
79 fprintf(
ioQQQ,
" insane cs returned by HCSAR_interp, values are\n");
127 static const double ap[5] = {-2113.113,729.0084,1055.397,854.632,938.9912};
128 static const double bp[5] = {-6783.515,-377.7190,724.1936,493.1107,735.7466};
129 static const double cp[5] = {-3049.719,226.2320,637.8630,388.5465,554.6369};
130 static const double dp[5] = {3514.5153,88.60169,-470.4055,-329.4914,-450.8459};
131 static const double ep[5] = {0.005251557,0.009059154,0.008725781,0.009952418,0.01098687};
132 static const double ae[5] = {-767.5859,-643.1189,-461.6836,-429.0543,-406.5285};
133 static const double be[5] = {-1731.9178,-1442.548,-1055.364,-980.3079,-930.9266};
134 static const double ce[5] = {-939.1834,-789.9569,-569.1451,-530.1974,-502.0939};
135 static const double de[5] = {927.4773,773.2008,564.3272,524.2944,497.7763};
136 static const double ee[5] = {-0.002528027,-0.003793665,-0.002122103,-0.002234207,-0.002317720};
137 static const double A[2] = {4.4394,0.0};
138 static const double B[2] = {0.8949,0.8879};
139 static const double C0[2] = {-0.6012,-0.2474};
140 static const double C1[2] = {-3.9710,-3.7562};
141 static const double C2[2] = {-4.2176,2.0491};
142 static const double D0[2] = {2.930,0.0539};
143 static const double D1[2] = {1.7990,3.4009};
144 static const double D2[2] = {4.9347,-1.7770};
232 fprintf(
ioQQQ,
" Insane levels in Hydcs123\n" );
243 else if( ipLow == 2 )
249 else if( ipLow == 3 )
256 fprintf(
ioQQQ,
" Insane levels in Hydcs123\n" );
261 Charge = (double)(nelem + 1);
263 ChargeSquared = Charge*Charge;
265 if( ipLow == 2 && ipHi == 3 )
273 fprintf(
ioQQQ,
" insane charge given to Hydcs123\n" );
276 else if( nelem == 1 )
284 else if( nelem <= 5 )
291 else if( nelem <= 11 )
298 else if( nelem <= 15 )
305 else if( nelem <= 17 )
326 TeUse =
MIN2(x,0.80);
327 x =
MAX2(0.025,TeUse);
333 Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*
pow2(x)*log(x) + de[j-1]*
334 exp(x) + ee[j-1]*log(x)/
pow2(x);
336 else if( chType ==
'p' )
338 Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*
pow2(x)*log(x) + dp[j-1]*
339 exp(x) + ep[j-1]*log(x)/
pow2(x);
344 fprintf(
ioQQQ,
" insane collision species given to Hydcs123\n" );
350 TeUse =
MIN2(x,0.80);
351 x =
MAX2(0.025,TeUse);
354 Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*
pow2(x)*log(x) +
355 de[k-1]*exp(x) + ee[k-1]*log(x)/
pow2(x);
359 Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*
pow2(x)*log(x) +
360 dp[k-1]*exp(x) + ep[k-1]*log(x)/
pow2(x);
369 slope = (Ratehigh - Ratelow)/(zhigh - zlow);
370 rate = slope*(Charge - zlow) + Ratelow;
372 rate = rate/ChargeSquared/Charge*1.0e-7;
374 templow = 0.025*27.211396*
TE1RYD/
EVRYD*ChargeSquared;
375 temphigh = 0.80*27.211396*
TE1RYD/
EVRYD*ChargeSquared;
377 temp =
MAX2(TeUse,templow);
386 else if( ipHi == 4 || ipHi == 5 || ipHi == 6 )
392 fprintf(
ioQQQ,
" insane charge given to Hydcs123\n" );
395 else if( nelem == 1 )
402 else if( nelem > 1 && nelem <= 5 )
407 Ratehigh =
C6cs123(ipLow,ipHi);
409 else if( nelem > 5 && nelem <= 9 )
416 else if( nelem > 9 && nelem <= 19 )
423 else if( nelem > 19 && nelem <= 25 )
448 slope = (Ratehigh - Ratelow)/(zhigh - zlow);
449 rate = slope*(Charge - zlow) + Ratelow;
466 return( Hydcs123_v );
472 EE = ChargeSquared*
EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp);
482 C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared;
483 D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared;
509 rate = (B[i] + D*
q)/expq + (A[i] + C*q - D*q*q)*
511 rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);
513 rate *= expq*gLo/gHi;
518 return( Hydcs123_v );
529 static const double a[3] = {-92.23774,-1631.3878,-6326.4947};
530 static const double b[3] = {-11.93818,-218.3341,-849.8927};
531 static const double c[3] = {0.07762914,1.50127,5.847452};
532 static const double d[3] = {78.401154,1404.8475,5457.9291};
533 static const double e[3] = {332.9531,5887.4263,22815.211};
550 t =
MIN2(TeUse,1.6e6);
553 if( i == 1 && j == 2 )
556 fprintf(
ioQQQ,
" Carbon VI 2s-1s not done in C6cs123\n" );
560 else if( i == 1 && j == 3 )
563 fprintf(
ioQQQ,
" Carbon VI 2p-1s not done in C6cs123\n" );
567 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) )
570 C6cs123_v = a[0] + b[0]*x + c[0]*
pow2(x)*sqrt(x) + d[0]*log(x) +
573 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) )
576 C6cs123_v = a[1] + b[1]*x + c[1]*
pow2(x)*sqrt(x) + d[1]*log(x) +
579 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) )
582 C6cs123_v = a[2] + b[2]*x + c[2]*
pow2(x)*sqrt(x) + d[2]*log(x) +
587 fprintf(
ioQQQ,
" insane levels for C VI n=1,2,3 !!!\n" );
601 static const double a[3] = {-12.5007,-187.2303,-880.18896};
602 static const double b[3] = {-1.438749,-22.17467,-103.1259};
603 static const double c[3] = {0.008219688,0.1318711,0.6043752};
604 static const double d[3] = {10.116516,153.2650,717.4036};
605 static const double e[3] = {45.905343,685.7049,3227.2836};
624 t =
MIN2(TeUse,1.585e7);
627 if( i == 1 && j == 2 )
630 fprintf(
ioQQQ,
" Ca XX 2s-1s not done in Ca20cs123\n" );
634 else if( i == 1 && j == 3 )
637 fprintf(
ioQQQ,
" Ca XX 2p-1s not done in Ca20cs123\n" );
641 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ))
644 Ca20cs123_v = a[0] + b[0]*x + c[0]*
pow2(x)*sqrt(x) + d[0]*log(x) +
647 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ))
650 Ca20cs123_v = a[1] + b[1]*x + c[1]*
pow2(x)*sqrt(x) + d[1]*log(x) +
653 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ))
656 Ca20cs123_v = a[2] + b[2]*x + c[2]*
pow2(x)*sqrt(x) + d[2]*log(x) +
661 fprintf(
ioQQQ,
" insane levels for Ca XX n=1,2,3 !!!\n" );
664 return( Ca20cs123_v );
675 static const double a[3] = {3.346644,151.2435,71.7095};
676 static const double b[3] = {0.5176036,20.05133,13.1543};
677 static const double c[3] = {-0.00408072,-0.1311591,-0.1099238};
678 static const double d[3] = {-3.064742,-129.8303,-71.0617};
679 static const double e[3] = {-11.87587,-541.8599,-241.2520};
696 t =
MIN2(TeUse,1.6e6);
699 if( i == 1 && j == 2 )
702 fprintf(
ioQQQ,
" Neon X 2s-1s not done in Ne10cs123\n" );
706 else if( i == 1 && j == 3 )
709 fprintf(
ioQQQ,
" Neon X 2p-1s not done in Ne10cs123\n" );
713 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) )
716 Ne10cs123_v = a[0] + b[0]*x + c[0]*
pow2(x)*sqrt(x) + d[0]*log(x) +
719 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) )
722 Ne10cs123_v = a[1] + b[1]*x + c[1]*
pow2(x)*sqrt(x) + d[1]*log(x) +
725 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) )
728 Ne10cs123_v = a[2] + b[2]*x + c[2]*
pow2(x)*sqrt(x) + d[2]*log(x) +
733 fprintf(
ioQQQ,
" insane levels for Ne X n=1,2,3 !!!\n" );
736 return( Ne10cs123_v );
745 static const double a[11]={0.12176209,0.32916723,0.46546497,0.044501688,
746 0.040523277,0.5234889,1.4903214,1.4215094,1.0295881,4.769306,9.7226127};
747 static const double b[11]={0.039936166,2.9711166e-05,-0.020835863,3.0508137e-04,
748 -2.004485e-15,4.41475e-06,1.0622666e-05,2.0538877e-06,0.80638448,2.0967075e-06,
750 static const double c[11]={143284.77,0.73158545,-2.159172,0.43254802,2.1338557,
751 8.9899702e-06,-2.9001451e-12,1.762076e-05,52741.735,-2153.1219,-3.3996921e-11};
773 else if( t > 5.0e05 )
780 if( i == 1 && j == 2 )
783 He2cs123_v = a[0] + b[0]*exp(-t/c[0]);
785 else if( i == 1 && j == 3 )
788 He2cs123_v = a[1] + b[1]*pow(t,c[1]);
790 else if( i == 1 && j == 4 )
793 He2cs123_v = a[2] + b[2]*log(t) + c[2]/log(t);
795 else if( i == 1 && j == 5 )
798 He2cs123_v = a[3] + b[3]*pow(t,c[3]);
800 else if( i == 1 && j == 6 )
803 He2cs123_v = a[4] + b[4]*pow(t,c[4]);
805 else if( i == 2 && j == 4 )
808 He2cs123_v = (a[5] + c[5]*t)/(1 + b[5]*t);
810 else if( i == 2 && j == 5 )
813 He2cs123_v = a[6] + b[6]*t + c[6]*t*t;
815 else if( i == 2 && j == 6 )
818 He2cs123_v = (a[7] + c[7]*t)/(1 + b[7]*t);
820 else if( i == 3 && j == 4 )
823 He2cs123_v = a[8] + b[8]*exp(-t/c[8]);
825 else if( i == 3 && j == 5 )
828 He2cs123_v = a[9] + b[9]*t + c[9]/t;
830 else if( i == 3 && j == 6 )
833 He2cs123_v = a[10] + b[10]*t + c[10]*t*t;
838 fprintf(
ioQQQ,
" insane levels for He II n=1,2,3 !!!\n" );
841 return( He2cs123_v );
852 static const double a[3] = {-4.238398,-238.2599,-1211.5237};
853 static const double b[3] = {-0.4448177,-27.06869,-136.7659};
854 static const double c[3] = {0.0022861,0.153273,0.7677703};
855 static const double d[3] = {3.303775,191.7165,972.3731};
856 static const double e[3] = {15.82689,878.1333,4468.696};
873 t =
MIN2(TeUse,1.585e7);
876 if( i == 1 && j == 2 )
879 fprintf(
ioQQQ,
" Fe XXVI 2s-1s not done in Fe26cs123\n" );
883 else if( i == 1 && j == 3 )
886 fprintf(
ioQQQ,
" Fe XXVI 2p-1s not done in Fe26cs123\n" );
890 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) )
893 Fe26cs123_v = a[0] + b[0]*x + c[0]*
pow2(x)*sqrt(x) + d[0]*log(x) +
896 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) )
899 Fe26cs123_v = a[1] + b[1]*x + c[1]*
pow2(x)*sqrt(x) + d[1]*log(x) +
902 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) )
905 Fe26cs123_v = a[2] + b[2]*x + c[2]*
pow2(x)*sqrt(x) + d[2]*log(x) +
910 fprintf(
ioQQQ,
" insane levels for Ca XX n=1,2,3 !!!\n" );
913 return( Fe26cs123_v );
974 double A, D, L, F, G, H;
975 double np, n, s, Z_3, xPlus, xMinus, y;
993 A = (8./3./s) * pow(np/s/n, 3.) * (0.184 - 0.04 * pow( s, -0.66)) * pow( 1. - 0.2*s/n/np, 1.+2.*s);
995 Z_3 = (double)(nelem + 1. - ipISO);
997 D = exp( - Z_3 * Z_3 / n / np / Ebar / Ebar );
999 L = log( (1. + 0.53 * Ebar * Ebar * n * np / Z_3 / Z_3) / (1. + 0.4*Ebar) );
1001 F = pow( 1. - 0.3 * s * D / n /np, 1. + 2.*s );
1003 G = 0.5*
POW3( Ebar * n * n / Z_3 / np );
1005 xPlus = 2. * Z_3 / ( n * n * Ebar * ( sqrt( 2. - n*n/np/np ) + 1. ) );
1006 xMinus = 2. * Z_3 / ( n * n * Ebar * ( sqrt( 2. - n*n/np/np ) - 1. ) );
1008 y = 1. / (1. - D * log ( 18. * s )/ 4. / s);
1010 H =
POW2( xMinus) * log( 1. + 2.*xMinus/3. ) / ( 2.*y + 1.5*xMinus );
1011 H -=
POW2( xPlus ) * log( 1. + 2.* xPlus/3. ) / ( 2.*y + 1.5*xPlus );
1014 cross_section = (A*D*L + F*G*H);
1019 stat_weight = 2. * n * n;
1021 stat_weight = 4. * n * n;
1031 STATIC void TestPercivalRichards(
void )
1036 for(
long i=0; i<5; i++ )
1038 double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.};
1044 for(
long i=0; i<5; i++ )
1046 double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.};
1058 long int ipCollider )
1069 if(
N_(ipLo) ==
N_(ipHi) )
1072 abs(
L_(ipLo) -
L_(ipHi) ) != 1 )
1088 abs(
L_(ipLo) -
L_(ipHi) ) != 1 )
1121 CStemp =
Hydcs123(ipLo + 1,ipHi + 1,nelem,
'e');
1125 CStemp =
Hydcs123(ipLo + 1,ipHi + 1,nelem,
'p');
1127 else if(
N_(ipLo) ==
N_(ipHi) )