236 double ForbiddenHe[5] = { 1.272E-4, 1E-20, 1E-20, 177.58, 0.327 };
238 A = ForbiddenHe[ipHi - 1];
253 A = (3.9061E-7) * pow( (
double)nelem+1., 10.419 );
266 A = ( 11.431 * pow((
double)nelem, 9.091) );
271 A = ( 383.42 * pow((
double)nelem, 7.8901) );
279 A = ( 0.11012 * pow((
double)nelem, 7.6954) );
298 else if( nelem>
ipHELIUM &&
L_(ipHi)==1 &&
S_(ipHi)==3 &&
299 L_(ipLo)==0 &&
S_(ipLo)==1 &&
N_(ipLo) <
N_(ipHi) )
301 A = 8.0E-3 * exp(9.283/sqrt((
double)
N_(ipLo))) * pow((
double)nelem,9.091) /
302 pow((
double)
N_(ipHi),2.877);
307 else if( nelem >
ipHELIUM &&
L_(ipHi)==0 &&
S_(ipHi)==1 &&
308 L_(ipLo)==1 &&
S_(ipLo)==3 &&
N_(ipLo) <
N_(ipHi) )
310 A = 2.416 * exp(-0.256*
N_(ipLo)) * pow((
double)nelem,9.159) / pow((
double)
N_(ipHi),3.111);
316 A *= (2.*(ipLo-3)+1.0)/3.0;
323 double As_3TripS_to_2TripS[9] = {
324 7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
325 1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
331 A = As_3TripS_to_2TripS[nelem-1];
340 A= 7.22E-9*pow((
double)nelem, 9.33);
358 A= 0.1834*pow((
double)nelem, 6.5735);
363 else if( nelem==
ipHELIUM && ipHi==4 && ipLo==3 )
377 else if( nelem==
ipHELIUM && ipHi==5 && ipLo==4 )
392 A = 44.326 * (1./3.) + 0.1146547 * (5./9.);
399 A = 1.459495 * (1./3.) + 3.6558e-5 * (5./9.);
406 A = 2.2297e-3 * (1./3.);
424 else if( ( ipLo==0 &&
L_(ipHi)==2 &&
S_(ipHi)==1 ) ||
427 (
N_(ipLo)==2 &&
L_(ipLo)==1 &&
S_(ipLo)==3 &&
N_(ipHi)>=3 &&
L_(ipHi)==1 &&
S_(ipHi)==3 ) ||
432 double f_params[5][4][3] = {
434 {9.360591E-007, -3.1937E-006, 3.5186E-006},
435 {4.631435E-007, -1.4973E-006, 1.4848E-006},
436 {2.493912E-007, -7.8658E-007, 7.3994E-007},
437 {1.476742E-007, -4.5953E-007, 4.1932E-007}},
439 {1.646733E-006, -2.0028E-006, -1.3552E-006},
440 {9.120593E-008, 3.1301E-007, -3.2190E-007},
441 {1.360965E-008, 1.1467E-007, 8.6977E-008},
442 {3.199421E-009, 4.5485E-008, 1.1016E-007}},
444 {1.646733E-006, -2.9720E-006, 1.5367E-006},
445 {9.120593E-008, -3.9037E-008, 3.9156E-008},
446 {1.360965E-008, 1.4671E-008, 1.5626E-008},
447 {3.199421E-009, 8.9806E-009, 1.2436E-008}},
449 {1.543812E-007, -2.8220E-007, 3.0261E-008},
450 {3.648237E-008, -6.6824E-008, 4.5879E-009},
451 {1.488556E-008, -2.7304E-008, 1.6628E-009},
452 {7.678610E-009, -1.4112E-008, 6.8453E-010}},
454 {1.543812E-007, -2.8546E-007, 4.6605E-008},
455 {3.648237E-008, -6.8422E-008, 1.7038E-008},
456 {1.488556E-008, -2.8170E-008, 8.5012E-009},
457 {7.678610E-009, -1.4578E-008, 4.6686E-009}}
470 long index_hi =
MIN2(
N_(ipHi), 6 ) - 3;
471 double f_lu =
POW2(nelem+1) * (
472 f_params[index_lo][index_hi][0] +
473 f_params[index_lo][index_hi][1]/(nelem+1) +
474 f_params[index_lo][index_hi][2]/
POW2(nelem+1) );
478 f_lu *= pow( 6. /
N_(ipHi), 3. );
484 A =
eina( gLo * f_lu, eWN, gHi );
503 long nelem ,
double Enerwn ,
505 double Eff_nupper,
long lHi,
long sHi,
long jHi,
507 double Eff_nlower,
long lLo,
long sLo,
long jLo )
513 long nHi, nLo, ipHi, ipLo;
521 nHi = (int)(Eff_nupper + 0.4);
522 nLo = (int)(Eff_nlower + 0.4);
525 ASSERT( fabs(Eff_nupper-(
double)nHi) < 0.4 );
526 ASSERT( fabs(Eff_nlower-(
double)nLo) < 0.4 );
529 if( (nHi==2) && (lHi==1) && (sHi==3) )
531 ASSERT( (jHi>=0) && (jHi<=2) );
536 if( (nLo==2) && (lLo==1) && (sLo==3) )
538 ASSERT( (jLo>=0) && (jLo<=2) );
556 if( (sHi == sLo) && (abs((
int)(lHi - lLo)) == 1) )
582 ASSERT( (lHi == 1) && (sHi == 1) );
586 Aul = (1.59208e10) / pow(Eff_nupper,3.0);
593 else if( lHi>=2 && lLo>=2 && nHi>nLo )
605 else if(
N_(ipHi)>10 &&
N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
608 double emisOscStr, x, a, b, c;
609 double extrapol_Params[2][4][4][3] = {
613 { 0.8267396 , 1.4837624 , -0.4615955 },
614 { 1.2738405 , 1.5841806 , -0.3022984 },
615 { 1.6128996 , 1.6842538 , -0.2393057 },
616 { 1.8855491 , 1.7709125 , -0.2115213 }},
618 { -1.4293664 , 2.3294080 , -0.0890470 },
619 { -0.3608082 , 2.3337636 , -0.0712380 },
620 { 0.3027974 , 2.3326252 , -0.0579008 },
621 { 0.7841193 , 2.3320138 , -0.0497094 }},
623 { 1.1341403 , 3.1702435 , -0.2085843 },
624 { 1.7915926 , 2.4942946 , -0.2266493 },
625 { 2.1979400 , 2.2785377 , -0.1518743 },
626 { 2.5018229 , 2.1925720 , -0.1081966 }},
628 { 0.0000000 , 0.0000000 , 0.0000000 },
629 { -2.6737396 , 2.9379143 , -0.3805367 },
630 { -1.4380124 , 2.7756396 , -0.2754625 },
631 { -0.6630196 , 2.6887253 , -0.2216493 }},
636 { 0.3075287 , 0.9087130 , -1.0387207 },
637 { 0.687069 , 1.1485864 , -0.6627317 },
638 { 0.9776064 , 1.3382024 , -0.5331906 },
639 { 1.2107725 , 1.4943721 , -0.4779232 }},
641 { -1.3659605 , 2.3262253 , -0.0306439 },
642 { -0.2899490 , 2.3279391 , -0.0298695 },
643 { 0.3678878 , 2.3266603 , -0.0240021 },
644 { 0.8427457 , 2.3249540 , -0.0194091 }},
646 { 1.3108281 , 2.8446367 , -0.1649923 },
647 { 1.8437692 , 2.2399326 , -0.2583398 },
648 { 2.1820792 , 2.0693762 , -0.1864091 },
649 { 2.4414052 , 2.0168255 , -0.1426083 }},
651 { 0.0000000 , 0.0000000 , 0.0000000 },
652 { -1.9219877 , 2.7689624 , -0.2536072 },
653 { -0.7818065 , 2.6595150 , -0.1895313 },
654 { -0.0665624 , 2.5955623 , -0.1522616 }},
662 else if( lLo==1 && lHi==0 )
666 else if( lLo==1 && lHi==2 )
676 ASSERT( (
int)((sHi-1)/2) == 0 || (
int)((sHi-1)/2) == 1 );
677 a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
678 b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
679 c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
683 emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
684 (2.*lLo+1)/(2.*lHi+1);
690 Aul *= (2.*(ipLo-3)+1.0)/9.0;
707 Aul *= (2.*(ipLo-3)+1.0)/9.0;
736 if( nLo==nHi && lHi<=2 && lLo<=2 )
743 Aul = 3.31E7 + 1.13E6 * pow((
double)nelem+1.,1.76);
745 Aul = 2.73E7 + 1.31E6 * pow((
double)nelem+1.,1.76);
747 Aul = 3.68E7 + 1.04E7 * exp(((
double)nelem+1.)/5.29);
750 fprintf(
ioQQQ,
" PROBLEM Impossible value for ipHi in he_1trans.\n");
758 Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
767 if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
769 Aul = 0.4 * 3.85E8 * pow((
double)nelem,1.6)/pow((
double)nHi,5.328);
772 else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
774 Aul = 1.95E4 * pow((
double)nelem,1.6) / pow((
double)nHi, 4.269);
777 else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
779 Aul = 6.646E7 * pow((
double)nelem,1.5) / pow((
double)nHi, 5.077);
783 ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
784 Aul = 3.9E6 * pow((
double)nelem,1.6) / pow((
double)nHi, 4.9);
785 if( (lHi >2) || (lLo > 2) )
795 else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
806 || ( ipLo ==
ipHe3s3S ) || ( nLo==4 && lLo==0 && sLo==3 ) )
811 ASSERT( (lHi == 1) && (sHi == 1) );
818 Aul = 1.375E10 * pow((
double)nelem, 3.9) / pow((
double)nHi,3.1);
824 ASSERT( (lHi == 1) && (sHi == 1) );
827 Aul = 5.0e8 * pow((
double)nelem,4.) / pow((
double)nHi, 2.889);
833 ASSERT( (lHi == 1) && (sHi == 3) );
837 Aul = 1.5 * 3.405E8 * pow((
double)nelem,4.) / pow((
double)nHi, 2.883);
839 Aul = 2.5 * 4.613E7 * pow((
double)nelem,4.) / pow((
double)nHi, 2.672);
841 Aul = 3.0 * 1.436E7 * pow((
double)nelem,4.) / pow((
double)nHi, 2.617);
851 RI2 =
scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(
double)(nelem));
853 Aul =
ritoa(lHi,lLo,nelem,Enerwn,RI2);
858 Aul *= (2.*(ipLo-3)+1.0)/9.0;
876 ASSERT( (sHi != sLo) || (abs((
int)(lHi - lLo)) != 1) );
886 if( Enerwn < 0. && Aul >
iso_ctrl.SmallA )
888 fprintf(
ioQQQ,
" he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
903 long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
904 double Ass, Att, Ast, Ats;
905 double SinHi, SinLo, CosHi, CosLo;
906 double HiMixingAngle, LoMixingAngle , error;
907 double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
918 if( ( sHi == 3 || sLo == 3 ) ||
919 ( abs(lHi - lLo) != 1 ) ||
921 ( lHi <= 1 || lLo <= 1 ) ||
922 ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
923 ( nHi > nLo && lHi != 1 && lLo != 1 ) )
936 HiMixingAngle = 0.01;
944 HiMixingAngle =
PI/4.;
949 LoMixingAngle = 0.01;
957 LoMixingAngle =
PI/4.;
961 ASSERT( ipHiTrip > ipLoTrip );
962 ASSERT( ipHiTrip > ipLoSing );
963 ASSERT( ipHiSing > ipLoTrip );
964 ASSERT( ipHiSing > ipLoSing );
966 SinHi = sin( HiMixingAngle );
967 SinLo = sin( LoMixingAngle );
968 CosHi = cos( HiMixingAngle );
969 CosLo = cos( LoMixingAngle );
982 temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
983 fssNew = Kss*
POW2( temp );
984 temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
985 fttNew = Ktt*
POW2( temp );
986 temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
987 fstNew = Kst*
POW2( temp );
988 temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
989 ftsNew = Kts*
POW2( temp );
1003 error = fabs( (
iso_sp[
ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()+
1005 (Ass+Ast+Ats+Att) - 1.f );
1009 fprintf(
ioQQQ,
"FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
1010 ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,