135 double rel_photon_energy,
168 double photon_energy,
529 inline double log10_prodxx(
long int lp,
double Ksqrd );
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)(l > lp ? 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)(l > lp ? 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" );
1334 long int rindx = 2*
q;
1336 if( rcsvV[rindx] == 0. )
1341 double Ksqrd = K * K;
1342 double n2 = (double)(n*n);
1344 double dd1 = (double)(2 * n);
1345 double dd2 = ( 1. + ( n2 * Ksqrd));
1350 double G1 = ((dd2 * GK) / dd1);
1362 else if( q == (n - 2) )
1364 double Ksqrd = (K*K);
1365 double n2 = (double)(n*n);
1370 double dd1 = (double) (2 * n);
1371 double dd2 = ( 1. + ( n2 * Ksqrd));
1372 double G1 = ((dd2 * GK) / dd1);
1377 double dd3 = (double)((2 * n) - 1);
1378 double dd4 = (double)(n - 1);
1379 double dd5 = (4. + (dd4 * dd2));
1380 double G2 = (dd3 * dd5 * G1);
1411 long int lp1 = (q + 1);
1412 long int lp2 = (q + 2);
1414 double Ksqrd = (K*K);
1415 double n2 = (double)(n * n);
1416 double lp1s = (double)(lp1 * lp1);
1417 double lp2s = (double)(lp2 * lp2);
1419 double d1 = (4. * n2);
1420 double d2 = (4. * lp1s);
1421 double d3 = (double)((lp1)*((2 *
q) + 3));
1422 double d4 = (1. + (n2 * Ksqrd));
1423 double d5 = (d1 - d2 + (d3 * d4));
1424 double d5_1 = d5 *
bhGp( (q+1), K, n, l, lp, rcsvV, GK );
1431 double d6 = (n2 - lp2s);
1432 double d7 = (1. + (lp1s * Ksqrd));
1433 double d8 = (d1 * d6 * d7);
1434 double d8_1 = d8 *
bhGp( (q+2), K, n, l, lp, rcsvV, GK );
1435 double d9 = (d5_1 - d8_1);
1460 ASSERT( rcsvV[rindx] != 0. );
1461 return rcsvV[rindx];
1485 long int rindx = 2*
q;
1487 if( rcsvV_mxq[rindx].q == 0 )
1498 double Ksqrd = (K * K);
1499 double n2 = (double)(n*n);
1501 double dd1 = (double) (2 * n);
1502 double dd2 = ( 1. + ( n2 * Ksqrd));
1503 double dd3 = dd2/dd1;
1516 rcsvV_mxq[rindx].
q = 1;
1517 rcsvV_mxq[rindx].
mx = G1_mx;
1521 else if( q == (n - 2) )
1533 double Ksqrd = (K*K);
1534 double n2 = (double)(n*n);
1535 double dd1 = (double)(2 * n);
1536 double dd2 = ( 1. + ( n2 * Ksqrd) );
1537 double dd3 = (dd2/dd1);
1538 double dd4 = (double) ((2 * n) - 1);
1539 double dd5 = (double) (n - 1);
1540 double dd6 = (4. + (dd5 * dd2));
1541 double dd7 = dd4 * dd6;
1567 rcsvV_mxq[rindx].
q = 1;
1568 rcsvV_mxq[rindx].
mx = G2_mx;
1586 long int lp1 = (q + 1);
1587 long int lp2 = (q + 2);
1589 double Ksqrd = (K * K);
1590 double n2 = (double)(n * n);
1591 double lp1s = (double)(lp1 * lp1);
1592 double lp2s = (double)(lp2 * lp2);
1594 double d1 = (4. * n2);
1595 double d2 = (4. * lp1s);
1596 double d3 = (double)((lp1)*((2 *
q) + 3));
1597 double d4 = (1. + (n2 * Ksqrd));
1599 double d5 = d1 - d2 + (d3 * d4);
1602 double d6 = (n2 - lp2s);
1605 double d7 = (1. + (lp1s * Ksqrd));
1608 double d8 = (d1 * d6 * d7);
1619 mx t0_mx =
bhGp_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1620 mx t1_mx =
bhGp_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1625 mx result_mx =
sub_mx( d9_mx, d10_mx );
1643 rcsvV_mxq[rindx].
q = 1;
1644 rcsvV_mxq[rindx].
mx = result_mx;
1650 ASSERT( rcsvV_mxq[rindx].q != 0 );
1651 rcsvV_mxq[rindx].
q = 1;
1652 return rcsvV_mxq[rindx].
mx;
1691 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1692 #pragma optimize("", off)
1709 long int rindx = 2*q+1;
1711 if( rcsvV[rindx] == 0. )
1721 else if( q == n - 2 )
1749 dd1 = (double) ((2 * n) - 1);
1752 dd2 = (1. + (n2 * Ksqrd));
1755 G2 = dd1 * dd2 * n1 * GK;
1764 long int lp2 = (q + 2);
1765 long int lp3 = (q + 3);
1767 double lp2s = (double)(lp2 * lp2);
1768 double lp3s = (double)(lp3 * lp3);
1783 double Ksqrd = (K*K);
1784 double n2 = (double)(n*n);
1785 double d1 = (4. * n2);
1786 double d2 = (4. * lp2s);
1787 double d3 = (double)(lp2)*((2*
q)+3);
1788 double d4 = (1. + (n2 * Ksqrd));
1790 double d5 = d1 - d2 + (d3 * d4);
1798 double d6 = (n2 - lp2s);
1800 double d7 = (1. + (lp3s * Ksqrd));
1802 double d8 = d1 * d6 * d7;
1803 double d9 = d5 *
bhGm( (q+1), K, n, l, lp, rcsvV, GK );
1804 double d10 = d8 *
bhGm( (q+2), K, n, l, lp, rcsvV, GK );
1805 double d11 = d9 - d10;
1829 ASSERT( rcsvV[rindx] != 0. );
1830 return rcsvV[rindx];
1833 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1834 #pragma optimize("", on)
1856 long int rindx = 2*q+1;
1858 if( rcsvV_mxq[rindx].q == 0 )
1863 mx result_mx = GK_mx;
1866 rcsvV_mxq[rindx].
q = 1;
1867 rcsvV_mxq[rindx].
mx = result_mx;
1873 else if( q == n - 2 )
1875 double Ksqrd = (K * K);
1876 double n1 = (double)n;
1877 double n2 = (double) (n*n);
1878 double dd1 = (double)((2 * n) - 1);
1879 double dd2 = (1. + (n2 * Ksqrd));
1881 double dd3 = (dd1*dd2*n1);
1901 rcsvV_mxq[rindx].
q = 1;
1902 rcsvV_mxq[rindx].
mx = G2_mx;
1920 long int lp2 = (q + 2);
1921 long int lp3 = (q + 3);
1923 double lp2s = (double)(lp2 * lp2);
1924 double lp3s = (double)(lp3 * lp3);
1925 double n2 = (double)(n*n);
1926 double Ksqrd = (K * K);
1932 double d1 = (4. * n2);
1933 double d2 = (4. * lp2s);
1934 double d3 = (double)(lp2)*((2*
q)+3);
1935 double d4 = (1. + (n2 * Ksqrd));
1937 double d5 = d1 - d2 + (d3 * d4);
1945 double d6 = (n2 - lp2s);
1946 double d7 = (1. + (lp3s * Ksqrd));
1948 double d8 = d1 * d6 * d7;
1958 mx d9_mx =
bhGm_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1959 mx d10_mx =
bhGm_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1962 mx result_mx =
sub_mx( d11_mx , d12_mx );
1963 rcsvV_mxq[rindx].
q = 1;
1964 rcsvV_mxq[rindx].
mx = result_mx;
1985 ASSERT( rcsvV_mxq[rindx].q != 0 );
1986 return rcsvV_mxq[rindx].
mx;
2065 double ld3 = (ld1 / ld2);
2086 double d2 = sqrt( ld3 * partprod );
2087 double d3 =
powi( (2 * n) , (l - n) );
2088 double d4 =
bhG( K, n, l, lp, rcsvV );
2089 double d5 = (d2 * d3);
2090 double d6 = (d5 * d4);
2094 ASSERT( ((n-l)-1) >= 0 );
2096 ASSERT( partprod != 0. );
2132 double d1 = (double)(2*n);
2133 double d2 = (double)(l-n);
2134 double Ksqrd = (K*K);
2181 double ld4 = (1./2.) * ( ld3 + ld1 - ld2 );
2189 double ld5 = d2 * log10( d1 );
2209 double ld6 = (ld5+ld4);
2212 mx dd1_mx =
bhG_mx( K, n, l, lp, rcsvV_mxq );
2215 double result =
unmxify( dd2_mx );
2231 double Ksqrd =(K*K);
2232 double partprod = 1.;
2234 for( s = 0; s <= lp; s = s + 1 )
2236 double s2 = (double)(s*s);
2246 partprod *= ( 1. + ( s2 * Ksqrd ) );
2338 if( n > 60 || np > 60 )
2372 double d1 =
hv( n, np, iz );
2374 double d3 =
pow3(d2);
2375 double lg = (double)(l > lp ? l : lp);
2376 double Two_L_Plus_One = (double)(2*l + 1);
2377 double d6 = lg / Two_L_Plus_One;
2378 double d7 =
hri( n, l, np, lp , iz );
2379 double d8 = d7 * d7;
2380 double result =
CONST_ONE * d3 * d6 * d8;
2385 fprintf(
ioQQQ,
"Principal Quantum Number `n' too large.\n");
2390 fprintf(
ioQQQ,
" The charge is impossible.\n");
2393 if( n < 1 || np < 1 || l >= n || lp >= np )
2395 fprintf(
ioQQQ,
" The quantum numbers are impossible.\n");
2400 fprintf(
ioQQQ,
" The principal quantum numbers are such that n <= n'.\n");
2418 double d1 =
hv( n, np , iz );
2420 double d3 =
pow3(d2);
2421 double lg = (double)(l > lp ? l : lp);
2422 double Two_L_Plus_One = (double)(2*l + 1);
2423 double d6 = lg / Two_L_Plus_One;
2424 double d7 =
hri_log10( n, l, np, lp , iz );
2425 double d8 = d7 * d7;
2426 double result =
CONST_ONE * d3 * d6 * d8;
2431 fprintf(
ioQQQ,
" The charge is impossible.\n");
2434 if( n < 1 || np < 1 || l >= n || lp >= np )
2436 fprintf(
ioQQQ,
" The quantum numbers are impossible.\n");
2441 fprintf(
ioQQQ,
" The principal quantum numbers are such that n <= n'.\n");
2497 inline double hv(
long int n,
long int nprime,
long int iz )
2501 double n1 = (double)n;
2503 double np1 = (double)nprime;
2504 double np2 = np1*np1;
2506 double izsqrd = (double)(iz*iz);
2508 double d1 = 1. / n2;
2509 double d2 = 1. / np2;
2510 double d3 = izsqrd * rmr *
EN1RYD;
2511 double d4 = d2 - d1;
2512 double result = d3 * d4;
2522 fprintf(
ioQQQ,
" The principal quantum numbers are such that n <= n'.\n");
2579 double Z = (double)iz;
2593 ASSERT( n > np || ( n == np && l == lp + 1 ));
2595 ASSERT( lp == l + 1 || lp == l - 1 );
2605 else if( l == lp - 1 )
2615 printf(
"BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2622 ld1 =
hrii(a, b, c, d ) /
Z;
2649 inline double hri_log10(
long int n,
long int l,
long int np,
long int lp ,
long int iz )
2664 double Z = (double)iz;
2672 ASSERT( n > np || ( n == np && l == lp + 1 ));
2674 ASSERT( lp == l + 1 || lp == l - 1 );
2684 else if( l == lp - 1 )
2694 printf(
"BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2706 STATIC double hrii(
long int n,
long int l,
long int np,
long int lp)
2719 long int a = 0, b = 0, c = 0;
2720 long int i1 = 0, i2 = 0, i3 = 0, i4 = 0;
2726 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
2727 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
2728 double d00 = 0., d01 = 0.;
2742 printf(
"BadMagic: Energy requirements not met.\n\n" );
2749 d5 = (double)(i1 - i2);
2751 d7 = (double)n * d6;
2755 else if( l == np && lp == (l - 1) )
2766 d1 = (double)( 2*n - 2 );
2767 d2 = (double)( 2*n - 1 );
2771 d5 = (double)(4 * n * (n - 1));
2773 d6 = (double)(i1 * i1);
2783 d12 = d4 * d8 * d11;
2797 for( i1 = -l; i1 <= l; i1 = i1 + 1 )
2799 d1 = (double)(n - i1);
2808 d5 = (double)( 4. * n * l );
2810 d6 = (double)( i3 * i3 );
2812 d8 =
powi( d7, l+1 );
2816 d9 = (double)( i3 ) / (double)( i4 );
2817 d10 =
powi( d9 , i4 );
2824 d14 = d4 * d8 * d10 * d13;
2859 else if( lp == l + 1 )
2865 printf(
" BadMagic: Don't know what to do here.\n\n");
2882 fsf =
fsff( n, l, np );
2926 d2 = (double)(n - np);
2929 d5 = (double)(n * np);
2934 d00 =
F21( a, b, c, y, A );
2987 d01 =
F21(a, b, c, y, A );
2996 d1 =
pow2( (
double)i1 );
2998 d2 =
pow2( (
double)i2 );
3035 double log10_fsf = 0.;
3049 printf(
"BadMagic: l'= l+1 for n'= n.\n\n" );
3060 long int i1 = n * n;
3061 long int i2 = l * l;
3063 double d1 = 3. / 2.;
3064 double d2 = (double)n;
3065 double d3 = (double)(i1 - i2);
3066 double d4 = sqrt(d3);
3067 double result = d1 * d2 * d4;
3073 else if( l == np && lp == l - 1 )
3084 double d1 = (double)((2*n-2)*(2*n-1));
3085 double d2 = sqrt( d1 );
3086 double d3 = (double)(4*n*(n-1));
3087 double d4 = (double)(2*n-1);
3090 double d8 =
powi(d7,n);
3092 double d10 = d4 - d9;
3093 double d11 = 0.25*d10;
3094 double result = (d2 * d8 * d11);
3103 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0., ld7 = 0.;
3124 for(
long int i1 = (n-l); i1 <= (n+l); i1++ )
3126 double d1 = (double)(i1);
3146 ld3 = 0.5 * (ld1 - ld2);
3158 double d1 = (double)(l+1);
3159 double d2 = (double)(4*n*l);
3160 double d3 = (double)(n-l);
3161 double d4 = log10(d2);
3162 double d5 = log10(d3);
3164 ld4 = d1 * (d4 - 2.*d5);
3180 ld5 = d2 * (d3 - d4);
3193 double d6 = 0.25*d5;
3204 ld7 = ld3 + ld4 + ld5 + ld6;
3206 result = pow( 10., ld7 );
3215 long int a = 0, b = 0, c = 0;
3216 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
3217 mx d00={0.0,0}, d01={0.0,0}, d02={0.0,0}, d03={0.0,0};
3223 else if( lp == l + 1 )
3229 printf(
" BadMagic: Don't know what to do here.\n\n");
3310 d2 = (double)(n - np);
3314 d5 = (double)(n * np);
3347 d00 =
F21_mx( a, b, c, y, A );
3400 d01 =
F21_mx(a, b, c, y, A );
3425 d1 = (double)((n - np)*(n -np));
3426 d2 = (double)((n + np)*(n + np));
3432 while ( fabs(d02.m) > 1.0e+25 )
3439 d03.m = d00.
m * (1. - (d02.m/d00.
m) *
powi( 10. , (d02.x - d00.
x) ) );
3441 result = pow( 10., (log10_fsf + d03.x) ) * d03.m;
3446 return fabs(result);
3465 long int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
3466 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
3492 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3528 d2 =
powi( (
double)i0 , i1 );
3534 i1 = n + np - 2*l - 2;
3535 d2 =
powi( (
double)i0 , i1 );
3542 d2 =
powi( (
double)i0 , i1 );
3557 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3565 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3573 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3581 printf(
"BadMagic: Relational error amongst n, l, n' and l'\n" );
3591 d5 = sqrt(d1)*sqrt(d2);
3615 double d0 = 0., d1 = 0.;
3616 double ld0 = 0., ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0.;
3632 d0 = (double)(2*l - 1);
3646 result = -(ld0 + ld1);
3655 d0 = (double)(4 * n * np);
3656 d1 = (double)(l + 1);
3657 result += d1 * log10(d0);
3667 d0 = (double)(n - np);
3668 d1 = (double)(n + np - 2*l - 2);
3669 result += d1 * log10(fabs(d0));
3674 d0 = (double)(n + np);
3675 d1 = (double)(-n - np);
3676 result += d1 * log10(d0);
3700 ld4 = 0.5*((ld0+ld1)-(ld2+ld3));
3795 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3796 long int i1 = 0, i2 = 0;
3805 d1 = (double)i1 / (
double)i2;
3807 d2 =
hv( n, np, iz );
3809 d4 =
hri( n, l, np, lp ,iz );
3812 d6 = d0 * d1 * d3 * d5;
3818 inline double OscStr_f_log10(
long int n ,
long int l ,
long int np ,
long int lp ,
long int iz )
3820 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3821 long int i1 = 0, i2 = 0;
3830 d1 = (double)i1 / (
double)i2;
3832 d2 =
hv( n, np, iz );
3837 d6 = d0 * d1 * d3 * d5;
3842 STATIC double F21(
long int a ,
long int b,
long int c,
double y,
char A )
3856 ASSERT( A ==
'a' || A ==
'b' );
3875 yV = (
double*)
CALLOC(
sizeof(
double), (size_t)(-a + 5) );
3950 d1 =
F21i(a, b, c, y, yV );
3959 mx result_mx = {0.0,0};
3968 ASSERT( A ==
'a' || A ==
'b' );
4062 result_mx =
F21i_log(a, b, c, y, yV );
4067 STATIC double F21i(
long int a,
long int b,
long int c,
double y,
double *yV )
4071 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
4072 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
4073 long int i1 = 0, i2 = 0;
4089 else if( yV[-a] != 0. )
4148 d2= (double)i1 * d1;
4153 d8=
F21i( (
long int)(a + 1), b, c, y, yV );
4154 d9=
F21i( (
long int)(a + 2), b, c, y, yV );
4172 mx result_mx = {0.0,0};
4174 if( yV[-a].
q != 0. )
4190 yV[-a].
mx = result_mx;
4195 double d1 = (double)b;
4196 double d2 = (double)c;
4197 double d3 = y * (d1/d2);
4203 result_mx.
m = 1. - d3;
4206 while ( fabs(result_mx.
m) > 1.0e+25 )
4208 result_mx.
m /= 1.0e+25;
4216 yV[-a].
mx = result_mx;
4265 mx d8 = {0.0,0}, d9 = {0.0,0}, d10 = {0.0,0}, d11 = {0.0,0};
4267 double db = (double)b;
4268 double d00 = (double)(a + 1 - c);
4269 double d0 = (double)(a + 1);
4271 double d2 = d0 * d1;
4272 double d3 = d2 / d00;
4274 double d5 = d00 + d4;
4275 double d6 = d5 / d00;
4287 d8=
F21i_log( (a + 1), b, c, y, yV );
4288 d9=
F21i_log( (a + 2), b, c, y, yV );
4293 d10.m = d8.
m * (1. - (d9.m/d8.
m) *
powi( 10., (d9.x - d8.
x)));
4308 result_mx.
x = d11.x;
4309 result_mx.
m = d11.m * (1. + (d10.m/d11.m) *
powi( 10. , (d10.x - d11.x) ) );
4316 while ( fabs(result_mx.
m) > 1.0e+25 )
4318 result_mx.
m /= 1.0e+25;
4324 yV[-a].
mx = result_mx;
4333 while( fabs(target.
m) > 1.0e+25 )
4335 target.
m /= 1.0e+25;
4338 while( fabs(target.
m) < 1.0e-25 )
4340 target.
m *= 1.0e+25;
4348 mx result = {0.0,0};
4353 result.
m = a.
m * (1. + (b.
m/a.
m) *
powi( 10. , (b.
x - a.
x) ) );
4365 mx result = {0.0,0};
4367 minusb.
m = -minusb.
m;
4369 result =
add_mx( a, minusb );
4377 mx result_mx = {0.0,0};
4388 return a_mx.
m *
powi( 10., a_mx.
x );
4393 mx result_mx = {0.0,0};
4395 while ( log10_a > 25.0 )
4401 while ( log10_a < -25.0 )
4407 result_mx.
m = pow(10., log10_a);
4414 mx result = {0.0,0};
4416 result.
m = (a.
m * b.
m);
4417 result.
x = (a.
x + b.
x);
4436 double partsum = 0.;
4437 for(
long int s = 1; s <= lp; s++ )
4439 double s2 =
pow2((
double)s);
4440 partsum += log10( 1. + s2*Ksqrd );