29 STATIC double ritoa(
long li,
long lf,
long nelem,
double k,
double RI2 );
49 d2 =
zJint * sin(theta),
60 long int rep = 0, ddiv, divsor;
69 if( (fabs(vv)) - (
int)(fabs(vv)) > 0.5 )
70 ddiv = (int)(fabs(vv)) + 1;
72 ddiv = (int)(fabs(vv));
74 divsor = ((ddiv == 0) ? 1 : ddiv);
78 for( rep = 0; rep < divsor; rep++ )
81 rl = (((double) rep)/((double) divsor)),
82 ru = (((
double) (rep+1))/((double) divsor)),
147 double nstar,
long int l,
148 double npstar,
long int lp,
152 double n_c = ((2.0 * nstar * npstar ) / ( nstar + npstar ));
154 double D_n = (nstar - npstar);
155 double D_l = (double) ( l - lp );
156 double lg = (double) ( (lp > l) ? lp : l);
160 double f = ( 1.0 - g );
161 double e = (( f >= 0.0) ? sqrt( f ) : 0.0 );
163 double x = (e * D_n);
164 double z = (-1.0 * x);
165 double v1 = (D_n + 1.0);
166 double v2 = (D_n - 1.0);
168 double d1,d2,d7,d8,d9,d34,d56,d6_1;
196 d2 = (n_c * n_c)/(2.0 * D_n);
198 d34 = (1.0 - ((D_l * lg)/n_c)) *
AngerJ( v1, z );
200 d56 = (1.0 + ((D_l * lg)/n_c)) *
AngerJ( v2, z );
204 d7 = (2./
PI) * sin( d6_1 ) * (1.0 - e);
206 d8 = d1 * d2 * ( (d34) - (d56) + d7 );
213 ASSERT( (l == lp + 1) || ( l == lp - 1) );
227 double As2nuFrom1S[28] = {1940.,1.82E+04,9.21E+04,3.30E+05,9.44E+05,2.31E+06,5.03E+06,1.00E+07,
228 1.86E+07,3.25E+07,5.42E+07,8.69E+07,1.34E+08,2.02E+08,2.96E+08,4.23E+08,5.93E+08,8.16E+08,
229 1.08E+09,1.43E+09,1.88E+09,2.43E+09,3.25E+09,3.95E+09,4.96E+09,6.52E+09,7.62E+09,9.94E+09};
235 double As2nuFrom3S[28] = {1.25E-06,5.53E-05,8.93E-04,8.05E-03,4.95E-02,2.33E-01,8.94E-01,2.95E+00,
236 8.59E+00,2.26E+01,5.49E+01,1.24E+02,2.64E+02,5.33E+02,1.03E+03,1.91E+03,3.41E+03,5.91E+03,
237 9.20E+03,1.50E+04,2.39E+04,3.72E+04,6.27E+04,8.57E+04,1.27E+05,2.04E+05,2.66E+05,4.17E+05};
251 double ForbiddenHe[5] = { 1.272E-4, 51.02, 1E-20, 177.58, 0.327 };
254 A = ForbiddenHe[ipHi - 1];
264 A = (3.9061E-7) * pow( (
double)nelem+1., 10.419 ) + As2nuFrom3S[nelem-2];
267 A = As2nuFrom1S[nelem-2];
277 A = ( 11.431 * pow((
double)nelem, 9.091) );
282 A = ( 383.42 * pow((
double)nelem, 7.8901) );
290 A = ( 0.11012 * pow((
double)nelem, 7.6954) );
310 else if( nelem>
ipHELIUM &&
L_(ipHi)==1 &&
S_(ipHi)==3 &&
311 L_(ipLo)==0 &&
S_(ipLo)==1 &&
N_(ipLo) <
N_(ipHi) )
313 A = 8.0E-3 * exp(9.283/sqrt((
double)
N_(ipLo))) * pow((
double)nelem,9.091) /
314 pow((
double)
N_(ipHi),2.877);
319 else if( nelem >
ipHELIUM &&
L_(ipHi)==0 &&
S_(ipHi)==1 &&
320 L_(ipLo)==1 &&
S_(ipLo)==3 &&
N_(ipLo) <
N_(ipHi) )
322 A = 2.416 * exp(-0.256*
N_(ipLo)) * pow((
double)nelem,9.159) / pow((
double)
N_(ipHi),3.111);
328 A *= (2.*(ipLo-3)+1.0)/3.0;
335 double As_3TripS_to_2TripS[9] = {
336 7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
337 1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
343 A = As_3TripS_to_2TripS[nelem-1];
352 A= 7.22E-9*pow((
double)nelem, 9.33);
370 A= 0.1834*pow((
double)nelem, 6.5735);
375 else if( nelem==
ipHELIUM && ipHi==4 && ipLo==3 )
389 else if( nelem==
ipHELIUM && ipHi==5 && ipLo==4 )
419 long nelem ,
double Enerwn ,
421 double Eff_nupper,
long lHi,
long sHi,
long jHi,
423 double Eff_nlower,
long lLo,
long sLo,
long jLo )
429 long nHi, nLo, ipHi, ipLo;
437 nHi = (int)(Eff_nupper + 0.4);
438 nLo = (int)(Eff_nlower + 0.4);
441 ASSERT( fabs(Eff_nupper-(
double)nHi) < 0.4 );
442 ASSERT( fabs(Eff_nlower-(
double)nLo) < 0.4 );
445 if( (nHi==2) && (lHi==1) && (sHi==3) )
447 ASSERT( (jHi>=0) && (jHi<=2) );
452 if( (nLo==2) && (lLo==1) && (sLo==3) )
454 ASSERT( (jLo>=0) && (jLo<=2) );
472 if( (sHi == sLo) && (abs((
int)(lHi - lLo)) == 1) )
498 ASSERT( (lHi == 1) && (sHi == 1) );
502 Aul = (1.59208e10) / pow(Eff_nupper,3.0);
509 else if( lHi>=2 && lLo>=2 && nHi>nLo )
529 else if(
N_(ipHi)>10 &&
N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
532 double emisOscStr, x, a, b, c;
533 double extrapol_Params[2][4][4][3] = {
537 { 0.8267396 , 1.4837624 , -0.4615955 },
538 { 1.2738405 , 1.5841806 , -0.3022984 },
539 { 1.6128996 , 1.6842538 , -0.2393057 },
540 { 1.8855491 , 1.7709125 , -0.2115213 }},
542 { -1.4293664 , 2.3294080 , -0.0890470 },
543 { -0.3608082 , 2.3337636 , -0.0712380 },
544 { 0.3027974 , 2.3326252 , -0.0579008 },
545 { 0.7841193 , 2.3320138 , -0.0497094 }},
547 { 1.1341403 , 3.1702435 , -0.2085843 },
548 { 1.7915926 , 2.4942946 , -0.2266493 },
549 { 2.1979400 , 2.2785377 , -0.1518743 },
550 { 2.5018229 , 2.1925720 , -0.1081966 }},
552 { 0.0000000 , 0.0000000 , 0.0000000 },
553 { -2.6737396 , 2.9379143 , -0.3805367 },
554 { -1.4380124 , 2.7756396 , -0.2754625 },
555 { -0.6630196 , 2.6887253 , -0.2216493 }},
560 { 0.3075287 , 0.9087130 , -1.0387207 },
561 { 0.687069 , 1.1485864 , -0.6627317 },
562 { 0.9776064 , 1.3382024 , -0.5331906 },
563 { 1.2107725 , 1.4943721 , -0.4779232 }},
565 { -1.3659605 , 2.3262253 , -0.0306439 },
566 { -0.2899490 , 2.3279391 , -0.0298695 },
567 { 0.3678878 , 2.3266603 , -0.0240021 },
568 { 0.8427457 , 2.3249540 , -0.0194091 }},
570 { 1.3108281 , 2.8446367 , -0.1649923 },
571 { 1.8437692 , 2.2399326 , -0.2583398 },
572 { 2.1820792 , 2.0693762 , -0.1864091 },
573 { 2.4414052 , 2.0168255 , -0.1426083 }},
575 { 0.0000000 , 0.0000000 , 0.0000000 },
576 { -1.9219877 , 2.7689624 , -0.2536072 },
577 { -0.7818065 , 2.6595150 , -0.1895313 },
578 { -0.0665624 , 2.5955623 , -0.1522616 }},
586 else if( lLo==1 && lHi==0 )
590 else if( lLo==1 && lHi==2 )
600 ASSERT( (
int)((sHi-1)/2) == 0 || (
int)((sHi-1)/2) == 1 );
601 a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
602 b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
603 c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
607 emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
608 (2.*lLo+1)/(2.*lHi+1);
614 Aul *= (2.*(ipLo-3)+1.0)/9.0;
623 RI2 =
scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(
double)(ipHELIUM));
626 Aul =
ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
631 Aul *= (2.*(ipLo-3)+1.0)/9.0;
660 if( nLo==nHi && lHi<=2 && lLo<=2 )
667 Aul = 3.31E7 + 1.13E6 * pow((
double)nelem+1.,1.76);
669 Aul = 2.73E7 + 1.31E6 * pow((
double)nelem+1.,1.76);
671 Aul = 3.68E7 + 1.04E7 * exp(((
double)nelem+1.)/5.29);
674 fprintf(
ioQQQ,
" PROBLEM Impossible value for ipHi in he_1trans.\n");
682 Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
691 if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
693 Aul = 0.4 * 3.85E8 * pow((
double)nelem,1.6)/pow((
double)nHi,5.328);
696 else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
698 Aul = 1.95E4 * pow((
double)nelem,1.6) / pow((
double)nHi, 4.269);
701 else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
703 Aul = 6.646E7 * pow((
double)nelem,1.5) / pow((
double)nHi, 5.077);
707 ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
708 Aul = 3.9E6 * pow((
double)nelem,1.6) / pow((
double)nHi, 4.9);
709 if( (lHi >2) || (lLo > 2) )
719 else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
735 ASSERT( (lHi == 1) && (sHi == 1) );
742 Aul = 1.375E10 * pow((
double)nelem, 3.9) / pow((
double)nHi,3.1);
748 ASSERT( (lHi == 1) && (sHi == 1) );
751 Aul = 5.0e8 * pow((
double)nelem,4.) / pow((
double)nHi, 2.889);
757 ASSERT( (lHi == 1) && (sHi == 3) );
761 Aul = 1.5 * 3.405E8 * pow((
double)nelem,4.) / pow((
double)nHi, 2.883);
763 Aul = 2.5 * 4.613E7 * pow((
double)nelem,4.) / pow((
double)nHi, 2.672);
765 Aul = 3.0 * 1.436E7 * pow((
double)nelem,4.) / pow((
double)nHi, 2.617);
775 RI2 =
scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(
double)(nelem));
777 Aul =
ritoa(lHi,lLo,nelem,Enerwn,RI2);
782 Aul *= (2.*(ipLo-3)+1.0)/9.0;
803 ASSERT( (sHi != sLo) || (abs((
int)(lHi - lLo)) != 1) );
815 fprintf(
ioQQQ,
" he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
830 long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
831 double Ass, Att, Ast, Ats;
832 double SinHi, SinLo, CosHi, CosLo;
833 double HiMixingAngle, LoMixingAngle , error;
834 double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
845 if( ( sHi == 3 || sLo == 3 ) ||
846 ( abs(lHi - lLo) != 1 ) ||
848 ( lHi <= 1 || lLo <= 1 ) ||
849 ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
850 ( nHi > nLo && lHi != 1 && lLo != 1 ) )
863 HiMixingAngle = 0.01;
871 HiMixingAngle =
PI/4.;
876 LoMixingAngle = 0.01;
884 LoMixingAngle =
PI/4.;
888 ASSERT( ipHiTrip > ipLoTrip );
889 ASSERT( ipHiTrip > ipLoSing );
890 ASSERT( ipHiSing > ipLoTrip );
891 ASSERT( ipHiSing > ipLoSing );
893 SinHi = sin( HiMixingAngle );
894 SinLo = sin( LoMixingAngle );
895 CosHi = cos( HiMixingAngle );
896 CosLo = cos( LoMixingAngle );
909 temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
910 fssNew = Kss*
POW2( temp );
911 temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
912 fttNew = Ktt*
POW2( temp );
913 temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
914 fstNew = Kst*
POW2( temp );
915 temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
916 ftsNew = Kts*
POW2( temp );
932 (Ass+Ast+Ats+Att) - 1.f );
936 fprintf(
ioQQQ,
"FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
937 ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
954 STATIC double ritoa(
long li,
long lf,
long nelem,
double k,
double RI2)
965 double fmean,mu,EinsteinA,w,RI2_cm;
975 lg = (lf > li) ? lf : li;
977 fmean = 2.0*mu*w*lg*RI2_cm/((3.0*
H_BAR) * (2.0*li+1.0));
989 double Enerwn, n_eff_hi, n_eff_lo;
992 double z4 =
POW2((
double)nelem);
997 Enerwn =
Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN;
1015 Aul =
he_1trans( nelem, Enerwn, n_eff_hi,
L_(ipLo)+1,
1016 S_(ipLo), -1, n_eff_lo,
L_(ipLo),
S_(ipLo), ipLo-3);
1020 Aul *= (2.*
L_(ipLo)+3.) *
S_(ipLo) / (4.*(double)
N_(ipHi)*(double)
N_(ipHi));
1025 Aul1 =
he_1trans( nelem, Enerwn, n_eff_hi,
L_(ipLo)-1,
1026 S_(ipLo), -1, n_eff_lo,
L_(ipLo),
S_(ipLo), ipLo-3);
1031 Aul += Aul1*(2.*
L_(ipLo)-1.) *
S_(ipLo) / (4.*(double)
N_(ipHi)*(double)
N_(ipHi));
1045 Aul =
he_1trans( nelem, -1.*Enerwn, n_eff_lo,
1046 L_(ipLo),
S_(ipLo), ipLo-3, n_eff_hi,
L_(ipHi),
S_(ipHi), ipHi-3);
1050 Aul =
he_1trans( nelem, Enerwn, n_eff_hi,
1051 L_(ipHi),
S_(ipHi), ipHi-3, n_eff_lo,
L_(ipLo),
S_(ipLo), ipLo-3);
1063 char chLine[chLine_LENGTH];
1068 long nelem, ipLo, ipHi, i, i1, i2, i3;
1088 fprintf(
ioQQQ,
" HelikeTransProbSetup opening he_transprob.dat:");
1090 ioDATA =
open_data(
"he_transprob.dat",
"r" );
1093 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1095 fprintf(
ioQQQ,
" HelikeTransProbSetup could not read first line of he_transprob.dat.\n");
1104 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1106 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1108 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
1129 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1132 while( chLine[0]==
'#' )
1134 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1142 if( i1<0 || i2<=i1 )
1144 fprintf(
ioQQQ,
" HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1151 for( i=0; i<1; ++i )
1153 if( (chTemp = strchr( chTemp,
'\t' )) == NULL )
1155 fprintf(
ioQQQ,
" HelikeTransProbSetup could not init he_transprob\n" );
1164 if( (chTemp = strchr( chTemp,
'\t' )) == NULL )
1166 fprintf(
ioQQQ,
" HelikeTransProbSetup could not scan he_transprob\n" );
1171 sscanf( chTemp ,
"%le" , &
TransProbs[nelem][i2][i1] );
1176 fprintf(
ioQQQ,
" HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1183 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1185 fprintf(
ioQQQ,
" HelikeTransProbSetup could not read last line of he_transprob.dat.\n");
1194 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1196 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1198 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );