51 double xsqrt , x15 ,
x2;
52 double energy = energy_ryd *
EVRYD;
58 energy_keV = energy/1000.0;
68 else if(energy >= 15.4 && energy < 18.)
70 cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/
x2);
72 cross_section =
MAX2( 0. , cross_section );
75 else if(energy >= 18. && energy <= 30.)
77 cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5);
80 else if(energy > 30. && energy <= 85.)
82 cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5);
88 cross_section = 45.57*(1 - 2.003/xsqrt - 4.806/x + 50.577/x15 - 171.044/x2 + 231.608/(xsqrt*
x2) - 81.885/(x*x2))/pow(energy_keV,3.5);
91 return( cross_section * 1e-24 );
144 long int ion , nshell , low , ihi , ipop , ip;
169 for( ip=low-1; ip < ihi; ip++ )
198 static const realnum csh2p[
NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f,
218 fprintf(
ioQQQ,
" OpacityCreateAll called but NOT evaluated since already done.\n" );
230 fprintf(
ioQQQ,
" OpacityCreateAll called, evaluating.\n" );
245 for( nelem=ipISO; nelem <
LIMELM; nelem++ )
325 dx)/
POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+
326 2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0*
327 dx)/
POW3(1.e0 + 2.e0*dx));
351 POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. +
352 x*(0.130471301 + x*0.000524906)));
437 for( nelem=2; nelem <
LIMELM; nelem++ )
472 crs =
ofit(eps,opart);
476 for( n=1; n < 6; n++ )
512 crs = 3.85e-18*(4.4*pow(
rfield.
AnuOrg[i]/thres,-1.5) - 3.38*
526 (0.2602325880970085 +
527 445.8558249365131*exp(-
rfield.
AnuOrg[i]/0.1009243952792674))*
541 " OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n",
549 fprintf(
ioQQQ,
"The COMPILE OPACITIES command is currently not supported\n" );
555 " Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" ,
opac.
nOpacTot);
588 for( i=ilo-1; i < ihi; i++ )
627 fprintf(
ioQQQ,
" Too many opacities were entered into OpacityCreateReilMan. Increase the value of NOP.\n" );
628 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl );
635 for( i=0; i < ncr; i++ )
637 en[i] = cross[i*2]/13.6f;
638 cs[i] = cross[(i+1)*2-1]*1e-18f;
645 " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" );
647 " The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n",
650 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl );
651 fprintf(
ioQQQ,
" The original energy (eV) and cross section (mb) arrays follow:\n" );
652 fprintf(
ioQQQ,
" " );
654 for( i=0; i < ncross; i++ )
656 fprintf(
ioQQQ,
"%11.4e", cross[i] );
659 fprintf(
ioQQQ,
"\n" );
663 slope = (cs[1] - cs[0])/(en[1] - en[0]);
670 for( i=low-1; i < ihi; i++ )
683 fprintf(
ioQQQ,
" OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" );
684 fprintf(
ioQQQ,
" The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n",
686 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl
688 fprintf(
ioQQQ,
" The lowest energy enterd in the array was%10.2e eV\n",
690 fprintf(
ioQQQ,
" The highest energy ever needed would be%10.2eeV\n",
692 fprintf(
ioQQQ,
" The lowest energy needed was%10.2eeV\n",
697 slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
706 fprintf(
ioQQQ,
" Internal logical error in OpacityCreateReilMan.\n" );
707 fprintf(
ioQQQ,
" The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n",
708 rfield.
anu[i]*13.6, i, en[ics-1], en[ics] );
710 fprintf(
ioQQQ,
" The previous energy (eV) was%10.2e\n",
713 fprintf(
ioQQQ,
" Here comes the energy array. ICS=%4ld\n",
716 for( j=0; j < ncr; j++ )
718 fprintf(
ioQQQ,
"%10.2e", en[j] );
720 fprintf(
ioQQQ,
"\n" );
722 fprintf(
ioQQQ,
" chLabl was %4.4s\n", chLabl );
742 static const double y[7][5] = {
743 {8.915,3995.,3.242,10.44,0.0},
744 {11.31,1498.,5.27,7.319,0.0},
745 {10.5,1.059e05,1.263,13.04,0.0},
746 {19.49,48.47,8.806,5.983,0.0},
747 {50.,4.244e04,0.1913,7.012,4.454e-02},
748 {110.5,0.1588,148.3,-3.38,3.589e-02},
749 {177.4,32.37,381.2,1.083,0.0}
751 static const double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
752 static const long l[7]={1,1,1,0,1,1,0};
769 for(
int i=0; i < 7; i++ )
774 for(
int i=0; i < 7; i++ )
778 q = 5.5 - 0.5*y[i][3] + l[i];
783 pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));
818 for( ion=0; ion < nelem; ion++ )
828 nelec = nelem+1 - ion;
831 for( nshell=0; nshell <
Heavy.
nsShells[nelem][ion]; nshell++ )
853 ASSERT( low <= ihi || low<5 );
856 for( ip=low-1; ip < ihi; ip++ )
863 cs =
t_ADfA::Inst().phfit(nelem+1,nelec,nshell+1,energy);
878 fprintf(
punch.
ipPoint,
"%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n",
905 fprintf(
ioQQQ,
" NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack, please consider increasing it.\n" );
906 fprintf(
ioQQQ,
" NOTE OpacityCreate1Element doubled memory allocation to %li.\n",
ndimOpacityStack );
939 ASSERT( crs > 0. && crs < 1e-10 );
948 ASSERT( crs > 0. && crs < 1e-10 );
962 ASSERT( crs > 0. && crs < 1e-10 );
975 ASSERT( crs > 0. && crs < 1e-10 );
993 (crs > 0. && crs < 1e-10) );
1006 ASSERT( crs > 0. && crs < 1e-10 );
1023 static double y2[
NCRS];
1024 static double crs[
NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
1025 2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
1026 3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
1027 1.567,1.233,.912,.629,.39,.19};
1028 static double ener[
NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
1029 0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
1030 0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
1031 0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
1032 0.5520,0.8557,1.7669};
1033 static bool lgFirst =
true;
1050 energy = freq - 0.05552;
1051 if( energy < ener[0] || energy > ener[
NCRS-1] )
1084 rayleh_v = (8.41e-25*
powi(ener,4) + 3.37e-24*
powi(ener,6))*
1088 else if( ener < 0.646 )
1090 rayleh_v = (8.41e-25*
powi(ener,4) + 3.37e-24*
powi(ener,6) +
1094 else if( ener >= 0.646 && ener < 1.0 )
1096 rayleh_v = fabs(0.74959-ener);