42 #define N_PHOTON_GFF 145L
64 fprintf(
ioQQQ,
" PROBLEM DISASTER - the kinetic temperature, %.3eK,"
65 " is above the upper limit of the code, %.3eK.\n",
67 fprintf(
ioQQQ,
" This calculation is aborting.\n Sorry.\n");
75 fprintf(
ioQQQ,
" PROBLEM DISASTER - the kinetic temperature, %.3eK,"
76 " is below the lower limit of the code, %.3eK.\n",
78 fprintf(
ioQQQ,
" This calculation is aborting.\n Sorry.\n");
116 fprintf(
ioQQQ,
" PROBLEM DISASTER - the kinetic temperature, %.3eK,"
117 " is above the upper limit of the code, %.3eK.\n",
119 fprintf(
ioQQQ,
" This calculation is aborting.\n Sorry.\n");
126 fprintf(
ioQQQ,
" PROBLEM DISASTER - the kinetic temperature, %.3eK,"
127 " is below the lower limit of the code, %.3eK.\n",
129 fprintf(
ioQQQ,
" This calculation is aborting.\n Sorry.\n");
147 static double tgffused=-1.,
149 static long int nff_defined=-1;
150 static long maxion = 0, oldmaxion = 0;
151 static double ttused = 0.;
152 static double edused = 0.;
153 static bool lgZLogSet =
false;
176 fprintf(
ioQQQ,
"tfidle called with a zero or negative electron density,%10.2e\n",
184 fprintf(
ioQQQ,
"tfidle called with a negative electron temperature,%10.2e\n",
192 for( nelem=0; nelem<
LIMELM; ++nelem )
221 for( i=1; i < 7; i++ )
285 if( fabs(1.-edused/
phycon.
te)>=0.00001 )
310 if( fabs(1.-tgffused/
phycon.
te)>=0.00001
347 for( nelem = 0; nelem <
LIMELM; nelem++ )
362 || maxion>oldmaxion )
364 static bool lgFirstRunDone =
false;
402 if( lgFirstRunDone ==
false )
405 lgFirstRunDone =
true;
414 for( ion=1; ion<=
LIMELM; ion++ )
420 fprintf(
ioQQQ,
" LinterpTable for GffTable called with te out of bounds \n");
432 fprintf(
ioQQQ,
" LinterpVector for GffTable called with photon energy out of bounds \n");
447 fprintf(
ioQQQ,
" tfidle; gaunt factors?" );
450 fprintf(
ioQQQ,
"%2f g 1 2=%10.2e%9.2ld flag%2f guv(1,n)%10.2e\n",
558 for( nelem=0; nelem <
LIMELM; nelem++ )
599 long i,i1,i2,i3,j,charge,GffMAGIC = 80321;
605 # define chLine_LENGTH 1000
615 TeGFF[N_TE_GFF-1] += 0.01f;
619 PhoGFF[i] = 0.125f*i - 8.f;
623 for( i = 1; i <=
LIMELM; i++ )
633 for( i = 1; i <=
LIMELM; i++ )
641 fprintf(
ioQQQ,
" FillGFF opening gauntff.dat:");
645 ioDATA =
open_data(
"gauntff.dat",
"r" );
649 fprintf(
ioQQQ,
" Defaulting to on-the-fly computation, ");
650 fprintf(
ioQQQ,
"but the code runs much faster if you compile gauntff.dat!\n");
657 for( charge=1; charge<=
LIMELM; charge++ )
675 fprintf(
ioQQQ,
" FillGFF could not read first line of gauntff.dat.\n");
683 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
686 " FillGFF: the version of gauntff.dat is not the current version.\n" );
688 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
693 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
695 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
700 for( charge = 1; charge <=
LIMELM; charge++ )
707 fprintf(
ioQQQ,
" FillGFF could not read first line of gauntff.dat.\n");
715 if( i1!=charge || i2!=i )
717 fprintf(
ioQQQ,
" FillGFF detected insanity in gauntff.dat.\n");
719 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
730 fprintf(
ioQQQ,
" FillGFF detected insanity in gauntff.dat.\n");
732 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
743 fprintf(
ioQQQ,
" FillGFF could not read first line of gauntff.dat.\n");
751 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
754 " FillGFF: the version of gauntff.dat is not the current version.\n" );
756 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
761 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
763 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
771 fprintf(
ioQQQ,
" - opened and read ok.\n");
781 fprintf(ioGFF,
"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
788 for( charge = 1; charge <=
LIMELM; charge++ )
792 fprintf(ioGFF,
"%li\t%li", charge, i );
800 fprintf(ioGFF,
"\t%f",
GauntFF[charge][i][j] );
802 fprintf(ioGFF,
"\n" );
807 fprintf(ioGFF,
"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
816 fprintf(
ioQQQ,
"FillGFF: compilation complete, gauntff.dat created.\n" );
825 enum {DEBUG_LOC=
false};
832 long logu, loggamma2;
834 for( logu=-4; logu<=4; logu++)
840 for(loggamma2=-4; loggamma2<=4; loggamma2++)
842 double SutherlandGff[9][9]=
843 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
844 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
845 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
846 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
847 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
848 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
849 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
850 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
851 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
856 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
860 fprintf(
ioQQQ,
" tfidle found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
861 logu, loggamma2, error );
877 double GauntAtLowerPho, GauntAtUpperPho;
878 static long int ipTe=-1, ipPho=-1;
890 fprintf(
ioQQQ,
" InterpolateGff called with te out of bounds \n");
903 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) );
934 fprintf(
ioQQQ,
" InterpolateGff called with photon energy out of bounds \n");
939 if( log10(ERyd) >
PhoGFF[i] && log10(ERyd) <=
PhoGFF[i+1] )
945 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) );
948 else if( log10(ERyd) <
PhoGFF[ipPho] )
951 while( log10(ERyd) <
PhoGFF[ipPho] && ipPho > 0)
954 else if( log10(ERyd) >
PhoGFF[ipPho+1] )
972 (
GauntFF[charge][ipPho+1][ipTe+1] -
GauntFF[charge][ipPho+1][ipTe]) +
GauntFF[charge][ipPho+1][ipTe];
975 (GauntAtUpperPho - GauntAtLowerPho) + GauntAtLowerPho;
1005 fhunt (v,ltb,x,&ipx);
1009 if( ipx == -1 || ipx == ltb )
1014 ASSERT( (ipx >=0) && (ipx < ltb-1) );
1015 ASSERT( x >= v[ipx] && x <= v[ipx+1]);
1017 frac = (x - v[ipx]) / (v[ipx+1] - v[ipx]);
1018 for( i=0; i<lta; i++ )
1020 a[i] = frac*t[i][ipx+1]+(1.f-frac)*t[i][ipx];
1035 if( yy[0] < v[0] || yy[ly-1] > v[ltb-1] )
1042 for(j = 1; j < ltb && n < ly; j++) {
1043 while (yl < v[j] && n < ly) {
1044 frac = (yl-v[j-1])/(v[j]-v[j-1]);
1045 for(i = 0; i < lta; i++)
1046 a[i][n] = frac*t[i][j]+(1.f-frac)*t[i][j-1];
1050 ASSERT( yy[n] > yy[n-1] );
1060 long int jl, jm, jh, in;
1064 up = (xx[n-1] > xx[0]);
1065 if(jl < 0 || jl >= n)
1073 if((x >= xx[jl]) == up)
1081 while ((x >= xx[jh]) == up)
1101 while ((x < xx[jl]) == up)
1117 if((x > xx[jm]) == up)