24STATIC double Hydcs123(
long int ilow,
long int ihigh,
long int iz,
long int chType);
37static 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 );
932 if( !
iso_ctrl.lgCS_therm_ave[ipISO] )
938 if( (
dense.eden > 10000.) && (
dense.eden < 1E10 ) )
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 );
1019 stat_weight = 2. * n * n;
1021 stat_weight = 4. * n * n;
1031STATIC 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 )
1065 if( !
iso_ctrl.lgColl_excite[ipISO] )
1070 if(
opac.lgCaseB_HummerStorey )
1072 if(
N_(ipLo) ==
N_(ipHi) )
1075 abs(
L_(ipLo) -
L_(ipHi) ) != 1 )
1096 abs(
L_(ipLo) -
L_(ipHi) ) != 1 )
1124 CStemp =
Hydcs123(ipLo + 1,ipHi + 1,nelem,
'e');
1128 CStemp =
Hydcs123(ipLo + 1,ipHi + 1,nelem,
'p');
1130 else if(
N_(ipLo) ==
N_(ipHi) )
1142 else if (abs(
L_(ipLo) -
L_(ipHi) ) == 1)
double qg32(double, double, double(*)(double))
bool fp_equal(sys_float x, sys_float y, int n=3)
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
realnum h_coll_str(long ipLo, long ipHi, long ipTe)
double CS_l_mixing_VF01(long int ipISO, long int nelem, long int n, long int l, long int lp, long int s, double temp, long int Collider)
double CS_l_mixing_PS64(long nelem, double tau, double target_charge, long n, long l, double gHi, long Collider)
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
double hydro_vs_deexcit(long ipISO, long nelem, long ipHi, long ipLo, double Aul)
double CS_VS80(long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider)
realnum HydroCSInterp(long int nelem, long int ipHi, long int ipLo, long int ipCollider)
static const realnum HCSTE[NHCSTE]
static double global_deltaE
STATIC double He2cs123(long int i, long int j)
STATIC double Therm_ave_coll_str_int_PR78(double EOverKT)
STATIC double Ca20cs123(long int i, long int j)
STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp)
STATIC double CS_PercivalRichards78(double Ebar)
STATIC double Hydcs123(long int ilow, long int ihigh, long int iz, long int chType)
STATIC realnum HCSAR_interp(int ipLo, int ipHi)
STATIC double Fe26cs123(long int i, long int j)
STATIC double C6cs123(long int i, long int j)
STATIC double Ne10cs123(long int i, long int j)
t_iso_sp iso_sp[NISO][LIMELM]
UNUSED const double BOHR_RADIUS_CM
UNUSED const double ELECTRON_MASS
UNUSED const double EN1RYD
UNUSED const double EVRYD
UNUSED const double PROTON_MASS
UNUSED const double EVDEGK
UNUSED const double COLL_CONST
UNUSED const double TE1RYD