562 long i1,i2,i3,j,charge,GffMAGIC = 100804;
573 for(
long i = 0; i <
N_TE_GFF; i++ )
591 for(
long i = 1; i <=
LIMELM; i++ )
601 for(
long i = 1; i <=
LIMELM; i++ )
606 if( !
rfield.lgCompileGauntFF )
609 fprintf(
ioQQQ,
" FillGFF opening gauntff.dat:");
613 ioDATA =
open_data(
"gauntff.dat",
"r" );
617 fprintf(
ioQQQ,
" Defaulting to on-the-fly computation, ");
618 fprintf(
ioQQQ,
"but the code runs much faster if you compile gauntff.dat!\n");
625 for( charge=1; charge<=
LIMELM; charge++ )
643 fprintf(
ioQQQ,
" FillGFF could not read first line of gauntff.dat.\n");
647 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
648 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
649 i3 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
654 " FillGFF: the version of gauntff.dat is not the current version.\n" );
656 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
661 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
663 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
668 for( charge = 1; charge <=
LIMELM; charge++ )
675 fprintf(
ioQQQ,
" FillGFF could not read first line of gauntff.dat.\n");
680 i1 = (long)
FFmtRead(chLine,&i3,
sizeof(chLine),&lgEOL);
681 i2 = (long)
FFmtRead(chLine,&i3,
sizeof(chLine),&lgEOL);
683 if( i1!=charge || i2!=i )
685 fprintf(
ioQQQ,
" FillGFF detected insanity in gauntff.dat.\n");
687 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
698 fprintf(
ioQQQ,
" FillGFF detected insanity in gauntff.dat.\n");
700 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
711 fprintf(
ioQQQ,
" FillGFF could not read first line of gauntff.dat.\n");
715 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
716 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
717 i3 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
722 " FillGFF: the version of gauntff.dat is not the current version.\n" );
724 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
729 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
731 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
739 fprintf(
ioQQQ,
" - opened and read ok.\n");
749 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",
756 for( charge = 1; charge <=
LIMELM; charge++ )
760 fprintf(ioGFF,
"%li\t%li", charge, i );
768 fprintf(ioGFF,
"\t%f",
GauntFF[charge][i][j] );
770 fprintf(ioGFF,
"\n" );
775 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",
784 fprintf(
ioQQQ,
"FillGFF: compilation complete, gauntff.dat created.\n" );
793 double tempsave =
phycon.te;
794 long logu, loggamma2;
796 for( logu=-4; logu<=4; logu++)
798 for(loggamma2=-4; loggamma2<=4; loggamma2++)
800 double SutherlandGff[9][9]=
801 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
802 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
803 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
804 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
805 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
806 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
807 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
808 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
809 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
813 double ERyd = pow(10.,(
double)(logu-loggamma2));
817 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
820 fprintf(
ioQQQ,
" PROBLEM DISASTER tfidle found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
821 logu, loggamma2, error );
837 double GauntAtLowerPho, GauntAtUpperPho;
838 static long int ipTe=-1, ipPho=-1;
850 fprintf(
ioQQQ,
" InterpolateGff called with te out of bounds \n");
894 fprintf(
ioQQQ,
" InterpolateGff called with photon energy out of bounds \n");
899 if( log10(ERyd) >
PhoGFF[i] && log10(ERyd) <=
PhoGFF[i+1] )
908 else if( log10(ERyd) <
PhoGFF[ipPho] )
911 while( log10(ERyd) <
PhoGFF[ipPho] && ipPho > 0)
914 else if( log10(ERyd) >
PhoGFF[ipPho+1] )
932 (
GauntFF[charge][ipPho+1][ipTe+1] -
GauntFF[charge][ipPho+1][ipTe]) +
GauntFF[charge][ipPho+1][ipTe];
935 (GauntAtUpperPho - GauntAtLowerPho) + GauntAtLowerPho;