135 double rel_photon_energy,
168 double photon_energy,
566 double rel_photon_energy,
596 double rel_photon_energy,
612 double electron_energy;
614 double xn_sqrd = (double)(n*n);
615 double z_sqrd = (double)(iz*iz);
616 double Z = (double)iz;
621 if( rel_photon_energy < 1.+FLT_EPSILON )
627 if( n < 1 || l >= n )
629 fprintf(
ioQQQ,
" The quantum numbers are impossible.\n");
635 fprintf(
ioQQQ,
" This value of n is too large.\n");
642 electron_energy = (rel_photon_energy-1.) * (z_sqrd/xn_sqrd);
643 k = sqrt( ( electron_energy ) );
647 dim_rcsvV = (((n * 2) - 1) + 1);
649 for( i=0; i<dim_rcsvV; ++i )
667 double rel_photon_energy,
678 long int dim_rcsvV_mxq;
680 mxq *rcsvV_mxq = NULL;
682 double electron_energy;
685 double xn_sqrd = (double)(n*n);
686 double z_sqrd = (double)(iz*iz);
687 double Z = (double)iz;
692 if( rel_photon_energy < 1.+FLT_EPSILON )
695 fprintf(
ioQQQ,
"PROBLEM IN HYDRO_BAUMAN: rel_photon_energy, n, l, iz: %e\t%li\t%li\t%li\n",
702 if( n < 1 || l >= n )
704 fprintf(
ioQQQ,
" The quantum numbers are impossible.\n");
710 electron_energy = (rel_photon_energy-1.) * (z_sqrd/xn_sqrd);
712 k = sqrt( ( electron_energy ) );
718 dim_rcsvV_mxq = (((n * 2) - 1) + 1);
721 rcsvV_mxq = (
mxq*)
CALLOC( (
size_t)dim_rcsvV_mxq,
sizeof(
mxq) );
723 t1 =
bh_log( K, n, l, rcsvV_mxq );
727 t1 =
MAX2( t1, 1.0e-250 );
734 fprintf(
ioQQQ,
"PROBLEM: Hydro_Bauman...t1\t%e\n", t1 );
763 for( lp = l - 1; lp <= l + 1; lp = lp + 2 )
800 for( lp = l - 1; lp <= l + 1; lp = lp + 2 )
843 double Two_L_Plus_One = (double)(2*l + 1);
844 double lg = (double)
max(l,lp);
846 double n2 = (double)(n*n);
848 double Ksqrd = (K*K);
862 double d2 = 1. + n2*Ksqrd;
863 double d5 =
bhg( K, n, l, lp, rcsvV );
864 double Theta = d2 * d5 * d5;
865 double d7 = (lg/Two_L_Plus_One) * Theta;
867 ASSERT( Two_L_Plus_One != 0. );
925 double n2 = (double)(n*n);
926 double Ksqrd = (K*K);
927 double Two_L_Plus_One = (double)(2*l + 1);
928 double lg = (double)
max(l,lp);
934 ASSERT( Two_L_Plus_One != 0. );
952 d2 = ( 1. + n2 * (Ksqrd) );
956 d5 =
bhg_log( K, n, l, lp, rcsvV_mxq );
958 d5 =
MAX2( d5, 1.0E-150 );
961 Theta = d2 * d5 * d5;
964 d7 = (lg/Two_L_Plus_One) * Theta;
1018 double n1 = (double)n;
1019 double n2 = (double)(n * n);
1020 double Ksqrd = K * K;
1023 double ld2 =
powi((
double)(4*n), n);
1024 double ld3 = exp(-(
double)(2 * n));
1034 double G0 =
SQRTPIBY2 * (8. * n1 * ld2 * ld3) / ld1;
1036 double d1 = sqrt( 1. - exp(( -2. *
PI )/ K ));
1037 double d2 =
powi(( 1. + (n2 * Ksqrd)), ( n + 2 ));
1038 double d3 = atan( n1 * K );
1039 double d4 = ((2. / K) * d3);
1040 double d5 = (double)( 2 * n );
1041 double d6 = exp( d5 - d4 );
1042 double GK = ( d6 /( d1 * d2 ) ) * G0;
1045 ASSERT( (l == lp - 1) || (l == lp + 1) );
1050 ASSERT( ((2*n) - 1) < 1755 );
1051 ASSERT( ((2*n) - 1) >= 0 );
1053 ASSERT( (1.0 / ld1) != 0. );
1078 double result =
bhGm( l, K, n, l, lp, rcsvV, GK );
1085 else if( l == lp + 1 )
1087 double result =
bhGp( l, K, n, l, lp, rcsvV, GK );
1094 printf(
"BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
1112 double log10_GK = 0.;
1113 double log10_G0 = 0.;
1115 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
1116 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0.;
1118 double n1 = (double)n;
1119 double n2 = n1 * n1;
1120 double Ksqrd = K * K;
1125 ASSERT( (l == lp - 1) || (l == lp + 1) );
1130 ASSERT( ((2*n) - 1) >= 0 );
1163 ld2 = n1 * log10( 4. * n1 );
1171 ld3 = (-(2. * n1)) * (
LOG10_E);
1182 log10_G0 = log10(
SQRTPIBY2 * 8. * n1) + ( (ld2 + ld3) - ld1);
1201 d1 = (1. - exp(-(2. *
PI )/ K ));
1202 ld4 = (1./2.) * log10( d1 );
1213 d2 = ( 1. + (n2 * Ksqrd));
1214 ld5 = (n1 + 2.) * log10( d2 );
1227 d3 = atan( n1 * K );
1235 d5 = (double) ( 2. * n1 );
1254 log10_GK = (ld6 -(ld4 + ld5)) + log10_G0;
1255 ASSERT( log10_GK != 0. );
1262 mx result_mx =
bhGm_mx( l, K, n, l, lp, rcsvV_mxq , GK_mx );
1268 else if( l == lp + 1 )
1270 mx result_mx =
bhGp_mx( l, K, n, l, lp, rcsvV_mxq , GK_mx );
1277 printf(
"BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
1280 printf(
"This code should be inaccessible\n\n" );
1336 long int rindx = 2*
q;
1338 if( rcsvV[rindx] == 0. )
1343 double Ksqrd = K * K;
1344 double n2 = (double)(n*n);
1346 double dd1 = (double)(2 * n);
1347 double dd2 = ( 1. + ( n2 * Ksqrd));
1352 double G1 = ((dd2 * GK) / dd1);
1364 else if(
q == (n - 2) )
1366 double Ksqrd = (K*K);
1367 double n2 = (double)(n*n);
1372 double dd1 = (double) (2 * n);
1373 double dd2 = ( 1. + ( n2 * Ksqrd));
1374 double G1 = ((dd2 * GK) / dd1);
1379 double dd3 = (double)((2 * n) - 1);
1380 double dd4 = (double)(n - 1);
1381 double dd5 = (4. + (dd4 * dd2));
1382 double G2 = (dd3 * dd5 * G1);
1413 long int lp1 = (
q + 1);
1414 long int lp2 = (
q + 2);
1416 double Ksqrd = (K*K);
1417 double n2 = (double)(n * n);
1418 double lp1s = (double)(lp1 * lp1);
1419 double lp2s = (double)(lp2 * lp2);
1421 double d1 = (4. * n2);
1422 double d2 = (4. * lp1s);
1423 double d3 = (double)((lp1)*((2 *
q) + 3));
1424 double d4 = (1. + (n2 * Ksqrd));
1425 double d5 = (d1 - d2 + (d3 * d4));
1426 double d5_1 = d5 *
bhGp( (
q+1), K, n, l, lp, rcsvV, GK );
1433 double d6 = (n2 - lp2s);
1434 double d7 = (1. + (lp1s * Ksqrd));
1435 double d8 = (d1 * d6 * d7);
1436 double d8_1 = d8 *
bhGp( (
q+2), K, n, l, lp, rcsvV, GK );
1437 double d9 = (d5_1 - d8_1);
1462 ASSERT( rcsvV[rindx] != 0. );
1463 return rcsvV[rindx];
1487 long int rindx = 2*
q;
1489 if( rcsvV_mxq[rindx].
q == 0 )
1500 double Ksqrd = (K * K);
1501 double n2 = (double)(n*n);
1503 double dd1 = (double) (2 * n);
1504 double dd2 = ( 1. + ( n2 * Ksqrd));
1505 double dd3 = dd2/dd1;
1518 rcsvV_mxq[rindx].
q = 1;
1519 rcsvV_mxq[rindx].
mx = G1_mx;
1523 else if(
q == (n - 2) )
1535 double Ksqrd = (K*K);
1536 double n2 = (double)(n*n);
1537 double dd1 = (double)(2 * n);
1538 double dd2 = ( 1. + ( n2 * Ksqrd) );
1539 double dd3 = (dd2/dd1);
1540 double dd4 = (double) ((2 * n) - 1);
1541 double dd5 = (double) (n - 1);
1542 double dd6 = (4. + (dd5 * dd2));
1543 double dd7 = dd4 * dd6;
1569 rcsvV_mxq[rindx].
q = 1;
1570 rcsvV_mxq[rindx].
mx = G2_mx;
1588 long int lp1 = (
q + 1);
1589 long int lp2 = (
q + 2);
1591 double Ksqrd = (K * K);
1592 double n2 = (double)(n * n);
1593 double lp1s = (double)(lp1 * lp1);
1594 double lp2s = (double)(lp2 * lp2);
1596 double d1 = (4. * n2);
1597 double d2 = (4. * lp1s);
1598 double d3 = (double)((lp1)*((2 *
q) + 3));
1599 double d4 = (1. + (n2 * Ksqrd));
1601 double d5 = d1 - d2 + (d3 * d4);
1604 double d6 = (n2 - lp2s);
1607 double d7 = (1. + (lp1s * Ksqrd));
1610 double d8 = (d1 * d6 * d7);
1621 mx t0_mx =
bhGp_mx( (
q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1622 mx t1_mx =
bhGp_mx( (
q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1627 mx result_mx =
sub_mx( d9_mx, d10_mx );
1645 rcsvV_mxq[rindx].
q = 1;
1646 rcsvV_mxq[rindx].
mx = result_mx;
1652 ASSERT( rcsvV_mxq[rindx].
q != 0 );
1653 rcsvV_mxq[rindx].
q = 1;
1654 return rcsvV_mxq[rindx].
mx;
1693#if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1694#pragma optimize("", off)
1711 long int rindx = 2*
q+1;
1713 if( rcsvV[rindx] == 0. )
1723 else if(
q == n - 2 )
1751 dd1 = (double) ((2 * n) - 1);
1754 dd2 = (1. + (n2 * Ksqrd));
1757 G2 = dd1 * dd2 * n1 * GK;
1766 long int lp2 = (
q + 2);
1767 long int lp3 = (
q + 3);
1769 double lp2s = (double)(lp2 * lp2);
1770 double lp3s = (double)(lp3 * lp3);
1785 double Ksqrd = (K*K);
1786 double n2 = (double)(n*n);
1787 double d1 = (4. * n2);
1788 double d2 = (4. * lp2s);
1789 double d3 = (double)(lp2)*((2*
q)+3);
1790 double d4 = (1. + (n2 * Ksqrd));
1792 double d5 = d1 - d2 + (d3 * d4);
1800 double d6 = (n2 - lp2s);
1802 double d7 = (1. + (lp3s * Ksqrd));
1804 double d8 = d1 * d6 * d7;
1805 double d9 = d5 *
bhGm( (
q+1), K, n, l, lp, rcsvV, GK );
1806 double d10 = d8 *
bhGm( (
q+2), K, n, l, lp, rcsvV, GK );
1807 double d11 = d9 - d10;
1831 ASSERT( rcsvV[rindx] != 0. );
1832 return rcsvV[rindx];
1835#if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1836#pragma optimize("", on)
1858 long int rindx = 2*
q+1;
1860 if( rcsvV_mxq[rindx].
q == 0 )
1865 mx result_mx = GK_mx;
1868 rcsvV_mxq[rindx].
q = 1;
1869 rcsvV_mxq[rindx].
mx = result_mx;
1875 else if(
q == n - 2 )
1877 double Ksqrd = (K * K);
1878 double n1 = (double)n;
1879 double n2 = (double) (n*n);
1880 double dd1 = (double)((2 * n) - 1);
1881 double dd2 = (1. + (n2 * Ksqrd));
1883 double dd3 = (dd1*dd2*n1);
1903 rcsvV_mxq[rindx].
q = 1;
1904 rcsvV_mxq[rindx].
mx = G2_mx;
1922 long int lp2 = (
q + 2);
1923 long int lp3 = (
q + 3);
1925 double lp2s = (double)(lp2 * lp2);
1926 double lp3s = (double)(lp3 * lp3);
1927 double n2 = (double)(n*n);
1928 double Ksqrd = (K * K);
1934 double d1 = (4. * n2);
1935 double d2 = (4. * lp2s);
1936 double d3 = (double)(lp2)*((2*
q)+3);
1937 double d4 = (1. + (n2 * Ksqrd));
1939 double d5 = d1 - d2 + (d3 * d4);
1947 double d6 = (n2 - lp2s);
1948 double d7 = (1. + (lp3s * Ksqrd));
1950 double d8 = d1 * d6 * d7;
1960 mx d9_mx =
bhGm_mx( (
q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1961 mx d10_mx =
bhGm_mx( (
q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1964 mx result_mx =
sub_mx( d11_mx , d12_mx );
1965 rcsvV_mxq[rindx].
q = 1;
1966 rcsvV_mxq[rindx].
mx = result_mx;
1987 ASSERT( rcsvV_mxq[rindx].
q != 0 );
1988 return rcsvV_mxq[rindx].
mx;
2067 double ld3 = (ld1 / ld2);
2088 double d2 = sqrt( ld3 * partprod );
2089 double d3 =
powi( (2 * n) , (l - n) );
2090 double d4 =
bhG( K, n, l, lp, rcsvV );
2091 double d5 = (d2 * d3);
2092 double d6 = (d5 * d4);
2096 ASSERT( ((n-l)-1) >= 0 );
2098 ASSERT( partprod != 0. );
2134 double d1 = (double)(2*n);
2135 double d2 = (double)(l-n);
2136 double Ksqrd = (K*K);
2183 double ld4 = (1./2.) * ( ld3 + ld1 - ld2 );
2191 double ld5 = d2 * log10( d1 );
2211 double ld6 = (ld5+ld4);
2214 mx dd1_mx =
bhG_mx( K, n, l, lp, rcsvV_mxq );
2217 double result =
unmxify( dd2_mx );
2233 double Ksqrd =(K*K);
2234 double partprod = 1.;
2236 for( s = 0; s <= lp; s = s + 1 )
2238 double s2 = (double)(s*s);
2248 partprod *= ( 1. + ( s2 * Ksqrd ) );
2340 if( n > 60 || np > 60 )
2374 double d1 =
hv( n, np, iz );
2376 double d3 =
pow3(d2);
2377 double lg = (double)(l > lp ? l : lp);
2378 double Two_L_Plus_One = (double)(2*l + 1);
2379 double d6 = lg / Two_L_Plus_One;
2380 double d7 =
hri( n, l, np, lp , iz );
2381 double d8 = d7 * d7;
2382 double result =
CONST_ONE * d3 * d6 * d8;
2387 fprintf(
ioQQQ,
"Principal Quantum Number `n' too large.\n");
2392 fprintf(
ioQQQ,
" The charge is impossible.\n");
2395 if( n < 1 || np < 1 || l >= n || lp >= np )
2397 fprintf(
ioQQQ,
" The quantum numbers are impossible.\n");
2402 fprintf(
ioQQQ,
" The principal quantum numbers are such that n <= n'.\n");
2420 double d1 =
hv( n, np , iz );
2422 double d3 =
pow3(d2);
2423 double lg = (double)(l > lp ? l : lp);
2424 double Two_L_Plus_One = (double)(2*l + 1);
2425 double d6 = lg / Two_L_Plus_One;
2426 double d7 =
hri_log10( n, l, np, lp , iz );
2427 double d8 = d7 * d7;
2428 double result =
CONST_ONE * d3 * d6 * d8;
2433 fprintf(
ioQQQ,
" The charge is impossible.\n");
2436 if( n < 1 || np < 1 || l >= n || lp >= np )
2438 fprintf(
ioQQQ,
" The quantum numbers are impossible.\n");
2443 fprintf(
ioQQQ,
" The principal quantum numbers are such that n <= n'.\n");
2499inline double hv(
long int n,
long int nprime,
long int iz )
2503 double n1 = (double)n;
2505 double np1 = (double)nprime;
2506 double np2 = np1*np1;
2508 double izsqrd = (double)(iz*iz);
2510 double d1 = 1. / n2;
2511 double d2 = 1. / np2;
2512 double d3 = izsqrd * rmr *
EN1RYD;
2513 double d4 = d2 - d1;
2514 double result = d3 * d4;
2524 fprintf(
ioQQQ,
" The principal quantum numbers are such that n <= n'.\n");
2581 double Z = (double)iz;
2595 ASSERT( n > np || ( n == np && l == lp + 1 ));
2597 ASSERT( lp == l + 1 || lp == l - 1 );
2607 else if( l == lp - 1 )
2617 printf(
"BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2624 ld1 =
hrii(a, b, c, d ) / Z;
2651inline double hri_log10(
long int n,
long int l,
long int np,
long int lp ,
long int iz )
2666 double Z = (double)iz;
2674 ASSERT( n > np || ( n == np && l == lp + 1 ));
2676 ASSERT( lp == l + 1 || lp == l - 1 );
2686 else if( l == lp - 1 )
2696 printf(
"BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2708STATIC double hrii(
long int n,
long int l,
long int np,
long int lp)
2721 long int a = 0, b = 0, c = 0;
2722 long int i1 = 0, i2 = 0, i3 = 0, i4 = 0;
2728 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
2729 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
2730 double d00 = 0., d01 = 0.;
2744 printf(
"BadMagic: Energy requirements not met.\n\n" );
2751 d5 = (double)(i1 - i2);
2753 d7 = (double)n * d6;
2757 else if( l == np && lp == (l - 1) )
2768 d1 = (double)( 2*n - 2 );
2769 d2 = (double)( 2*n - 1 );
2773 d5 = (double)(4 * n * (n - 1));
2775 d6 = (double)(i1 * i1);
2785 d12 = d4 * d8 * d11;
2799 for( i1 = -l; i1 <= l; i1 = i1 + 1 )
2801 d1 = (double)(n - i1);
2810 d5 = (double)( 4. * n * l );
2812 d6 = (double)( i3 * i3 );
2814 d8 =
powi( d7, l+1 );
2818 d9 = (double)( i3 ) / (double)( i4 );
2819 d10 =
powi( d9 , i4 );
2826 d14 = d4 * d8 * d10 * d13;
2861 else if( lp == l + 1 )
2867 printf(
" BadMagic: Don't know what to do here.\n\n");
2884 fsf =
fsff( n, l, np );
2928 d2 = (double)(n - np);
2931 d5 = (double)(n * np);
2936 d00 =
F21( a, b, c, y, A );
2989 d01 =
F21(a, b, c, y, A );
2998 d1 =
pow2( (
double)i1 );
3000 d2 =
pow2( (
double)i2 );
3037 double log10_fsf = 0.;
3051 printf(
"BadMagic: l'= l+1 for n'= n.\n\n" );
3062 long int i1 = n * n;
3063 long int i2 = l * l;
3065 double d1 = 3. / 2.;
3066 double d2 = (double)n;
3067 double d3 = (double)(i1 - i2);
3068 double d4 = sqrt(d3);
3069 double result = d1 * d2 * d4;
3075 else if( l == np && lp == l - 1 )
3086 double d1 = (double)((2*n-2)*(2*n-1));
3087 double d2 = sqrt( d1 );
3088 double d3 = (double)(4*n*(n-1));
3089 double d4 = (double)(2*n-1);
3092 double d8 =
powi(d7,n);
3094 double d10 = d4 - d9;
3095 double d11 = 0.25*d10;
3096 double result = (d2 * d8 * d11);
3105 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0., ld7 = 0.;
3126 for(
long int i1 = (n-l); i1 <= (n+l); i1++ )
3128 double d1 = (double)(i1);
3148 ld3 = 0.5 * (ld1 - ld2);
3160 double d1 = (double)(l+1);
3161 double d2 = (double)(4*n*l);
3162 double d3 = (double)(n-l);
3163 double d4 = log10(d2);
3164 double d5 = log10(d3);
3166 ld4 = d1 * (d4 - 2.*d5);
3182 ld5 = d2 * (d3 - d4);
3195 double d6 = 0.25*d5;
3206 ld7 = ld3 + ld4 + ld5 + ld6;
3208 result = pow( 10., ld7 );
3217 long int a = 0, b = 0, c = 0;
3218 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
3219 mx d00={0.0,0}, d01={0.0,0}, d02={0.0,0}, d03={0.0,0};
3225 else if( lp == l + 1 )
3231 printf(
" BadMagic: Don't know what to do here.\n\n");
3312 d2 = (double)(n - np);
3316 d5 = (double)(n * np);
3349 d00 =
F21_mx( a, b, c, y, A );
3402 d01 =
F21_mx(a, b, c, y, A );
3427 d1 = (double)((n - np)*(n -np));
3428 d2 = (double)((n + np)*(n + np));
3434 while ( fabs(d02.m) > 1.0e+25 )
3441 d03.m = d00.
m * (1. - (d02.m/d00.
m) *
powi( 10. , (d02.x - d00.
x) ) );
3443 result = pow( 10., (log10_fsf + d03.x) ) * d03.m;
3448 return fabs(result);
3451 printf(
" This code should be inaccessible\n\n");
3470 long int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
3471 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
3497 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3533 d2 =
powi( (
double)i0 , i1 );
3539 i1 = n + np - 2*l - 2;
3540 d2 =
powi( (
double)i0 , i1 );
3547 d2 =
powi( (
double)i0 , i1 );
3562 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3570 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3578 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3586 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3596 d5 = sqrt(d1)*sqrt(d2);
3620 double d0 = 0., d1 = 0.;
3621 double ld0 = 0., ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0.;
3637 d0 = (double)(2*l - 1);
3651 result = -(ld0 + ld1);
3660 d0 = (double)(4 * n * np);
3661 d1 = (double)(l + 1);
3662 result += d1 * log10(d0);
3672 d0 = (double)(n - np);
3673 d1 = (double)(n + np - 2*l - 2);
3674 result += d1 * log10(fabs(d0));
3679 d0 = (double)(n + np);
3680 d1 = (double)(-n - np);
3681 result += d1 * log10(d0);
3705 ld4 = 0.5*((ld0+ld1)-(ld2+ld3));
3800 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3801 long int i1 = 0, i2 = 0;
3810 d1 = (double)i1 / (
double)i2;
3812 d2 =
hv( n, np, iz );
3814 d4 =
hri( n, l, np, lp ,iz );
3817 d6 = d0 * d1 * d3 * d5;
3823inline double OscStr_f_log10(
long int n ,
long int l ,
long int np ,
long int lp ,
long int iz )
3825 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3826 long int i1 = 0, i2 = 0;
3835 d1 = (double)i1 / (
double)i2;
3837 d2 =
hv( n, np, iz );
3842 d6 = d0 * d1 * d3 * d5;
3847STATIC double F21(
long int a ,
long int b,
long int c,
double y,
char A )
3861 ASSERT( A ==
'a' || A ==
'b' );
3880 yV = (
double*)
CALLOC(
sizeof(
double), (
size_t)(-a + 5) );
3955 d1 =
F21i(a, b, c, y, yV );
3964 mx result_mx = {0.0,0};
3973 ASSERT( A ==
'a' || A ==
'b' );
4067 result_mx =
F21i_log(a, b, c, y, yV );
4072STATIC double F21i(
long int a,
long int b,
long int c,
double y,
double *yV )
4076 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
4077 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
4078 long int i1 = 0, i2 = 0;
4094 else if( yV[-a] != 0. )
4153 d2= (double)i1 * d1;
4158 d8=
F21i( (
long int)(a + 1), b, c, y, yV );
4159 d9=
F21i( (
long int)(a + 2), b, c, y, yV );
4177 mx result_mx = {0.0,0};
4179 if( yV[-a].
q != 0. )
4195 yV[-a].
mx = result_mx;
4200 double d1 = (double)b;
4201 double d2 = (double)c;
4202 double d3 = y * (d1/d2);
4208 result_mx.
m = 1. - d3;
4211 while ( fabs(result_mx.
m) > 1.0e+25 )
4213 result_mx.
m /= 1.0e+25;
4221 yV[-a].
mx = result_mx;
4270 mx d8 = {0.0,0}, d9 = {0.0,0}, d10 = {0.0,0}, d11 = {0.0,0};
4272 double db = (double)b;
4273 double d00 = (double)(a + 1 - c);
4274 double d0 = (double)(a + 1);
4276 double d2 = d0 * d1;
4277 double d3 = d2 / d00;
4279 double d5 = d00 + d4;
4280 double d6 = d5 / d00;
4292 d8=
F21i_log( (a + 1), b, c, y, yV );
4293 d9=
F21i_log( (a + 2), b, c, y, yV );
4298 d10.m = d8.
m * (1. - (d9.m/d8.
m) *
powi( 10., (d9.x - d8.
x)));
4313 result_mx.
x = d11.x;
4314 result_mx.
m = d11.m * (1. + (d10.m/d11.m) *
powi( 10. , (d10.x - d11.x) ) );
4321 while ( fabs(result_mx.
m) > 1.0e+25 )
4323 result_mx.
m /= 1.0e+25;
4329 yV[-a].
mx = result_mx;
4338 while( fabs(target.
m) > 1.0e+25 )
4340 target.
m /= 1.0e+25;
4343 while( fabs(target.
m) < 1.0e-25 )
4345 target.
m *= 1.0e+25;
4353 mx result = {0.0,0};
4358 result.
m = a.
m * (1. + (b.
m/a.
m) *
powi( 10. , (b.
x - a.
x) ) );
4370 mx result = {0.0,0};
4372 minusb.
m = -minusb.
m;
4374 result =
add_mx( a, minusb );
4382 mx result_mx = {0.0,0};
4393 return a_mx.
m *
powi( 10., a_mx.
x );
4398 mx result_mx = {0.0,0};
4400 while ( log10_a > 25.0 )
4406 while ( log10_a < -25.0 )
4412 result_mx.
m = pow(10., log10_a);
4419 mx result = {0.0,0};
4421 result.
m = (a.
m * b.
m);
4422 result.
x = (a.
x + b.
x);
4441 double partsum = 0.;
4442 for(
long int s = 1; s <= lp; s++ )
4444 double s2 =
pow2((
double)s);
4445 partsum += log10( 1. + s2*Ksqrd );
#define DEBUG_ENTRY(funcname)
double powi(double, long int)
double hri(long int n, long int l, long int np, long int lp, long int iz)
mx mxify_log10(double log10_a)
mx add_mx(const mx &a, const mx &b)
double H_photo_cs_log10(double photon_energy, long int n, long int l, long int iz)
STATIC double H_Einstein_A_lin(long int n, long int l, long int np, long int lp, long int iz)
double OscStr_f(long int n, long int l, long int np, long int lp, long int iz)
STATIC mx F21i_log(long int a, long int b, long int c, double y, mxq *yV)
static const double PHYSICAL_CONSTANT_TWO
STATIC double fsff(long int n, long int l, long int np)
double local_product(double K, long int lp)
STATIC mx bhGm_mx(long int q, double K, long int n, long int l, long int lp, mxq *rcsvV_mxq, const mx &GK_mx)
double unmxify(const mx &a_mx)
double hv(long int n, long int nprime, long int iz)
void normalize_mx(mx &target)
STATIC double log10_fsff(long int n, long int l, long int np)
STATIC double hrii(long int n, long int l, long int np, long int lp)
double OscStr_f_log10(long int n, long int l, long int np, long int lp, long int iz)
STATIC double bh(double k, long int n, long int l, double *rcsvV)
double H_Einstein_A_log10(long int n, long int l, long int np, long int lp, long int iz)
STATIC double bhintegrand_log(double k, long int n, long int l, long int lp, mxq *rcsvV_mxq)
mx sub_mx(const mx &a, const mx &b)
STATIC double bhg(double K, long int n, long int l, long int lp, double *rcsvV)
STATIC mx F21_mx(long int a, long int b, long int c, double y, char A)
double hri_log10(long int n, long int l, long int np, long int lp, long int iz)
STATIC double bhGp(long int q, double K, long int n, long int l, long int lp, double *rcsvV, double GK)
STATIC mx bhG_mx(double K, long int n, long int l, long int lp, mxq *rcsvV_mxq)
mx mult_mx(const mx &a, const mx &b)
STATIC double bh_log(double k, long int n, long int l, mxq *rcsvV_mxq)
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
double log10_prodxx(long int lp, double Ksqrd)
STATIC double bhg_log(double K, long int n, long int l, long int lp, mxq *rcsvV_mxq)
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
static const double CONST_ONE
STATIC double bhG(double K, long int n, long int l, long int lp, double *rcsvV)
STATIC double F21i(long int a, long int b, long int c, double y, double *yV)
STATIC double bhintegrand(double k, long int n, long int l, long int lp, double *rcsvV)
STATIC double F21(long int a, long int b, long int c, double y, char A)
STATIC double H_photo_cs_lin(double rel_photon_energy, long int n, long int l, long int iz)
STATIC mx bhGp_mx(long int q, double K, long int n, long int l, long int lp, mxq *rcsvV_mxq, const mx &GK_mx)
STATIC double bhGm(long int q, double K, long int n, long int l, long int lp, double *rcsvV, double GK)
STATIC double hrii_log(long int n, long int l, long int np, long int lp)
UNUSED const double BOHR_RADIUS_CM
UNUSED const double LOG10_E
UNUSED const double HPLANCK
UNUSED const double FINE_STRUCTURE
UNUSED const double SPEEDLIGHT
UNUSED const double ELECTRON_MASS
UNUSED const double SQRTPIBY2
UNUSED const double EN1RYD
UNUSED const double PROTON_MASS
double lfactorial(long n)
static const int NPRE_FACTORIAL