/* Copyright (c) 2014, Peter A.M. van Hoof. This program is provided 'as-is', without any express or implied warranty. In no event will the author be held liable for any damages arising from the use of this program. Permission is granted to anyone to use this program for any purpose, including commercial applications, and to alter it and redistribute it freely, subject to the following restrictions: 1. The origin of this program must not be misrepresented; you must not claim that you created the original program. If you use this program in a product, an acknowledgment in the product documentation would be appreciated but is not required. 2. Altered program versions must be plainly marked as such, and must not be misrepresented as being the original program. 3. This notice may not be removed or altered from any further distribution. Peter A.M. van Hoof Royal Observatory of Belgium Ringlaan 3 B-1180 Brussels Belgium p.vanhoof@oma.be */ /* If you use any of these data in a scientific publication, please refer to van Hoof P.A.M., Williams R.J.R., Volk K., Chatzikos M., Ferland G.J., Lykins M., Porter R.L., Wang Y. Accurate determination of the free-free Gaunt factor, I -- non-relativistic Gaunt factors 2014, MNRAS, 444, 420 van Hoof P.A.M., Ferland G.J., Williams R.J.R., Volk K., Chatzikos M., Lykins M., Porter R.L. Accurate determination of the free-free Gaunt factor, II -- relativistic Gaunt factors 2015, MNRAS, 449, 2112 */ #include #include #include #include static const long gaunt_magic = 20141008L; static const char* gauntff_file = "gauntff_freqint_Z01.dat"; #define NP_GAM2 161 double log_g2[NP_GAM2], gff[NP_GAM2]; double lg_gam2_min, step, lg_gam2_max; size_t min(size_t x, size_t y) { return ( x < y ) ? x : y; } size_t max(size_t x, size_t y) { return ( x > y ) ? x : y; } void read_table(const char* fnam) { char* ptr; long magic; size_t np_gam2, i, ipg2, len; FILE* io = fopen( fnam, "r" ); if( io == NULL ) { fprintf( stderr, "failed to open file %s\n", fnam ); exit(1); } ptr = NULL; /* skip copyright statement */ while( 1 ) { getline( &ptr, &len, io ); if( ptr[0] != '#' ) break; free( ptr ); ptr = NULL; } /* read magic number */ sscanf( ptr, "%ld", &magic ); free( ptr ); ptr = NULL; if( magic != gaunt_magic ) { fprintf( stderr, "read_table() found wrong magic number in file %s.\n", fnam ); exit(1); } /* read dimensions of the table */ getline( &ptr, &len, io ); sscanf( ptr, "%ld", &np_gam2 ); free( ptr ); ptr = NULL; assert( np_gam2 == NP_GAM2 ); /* read start value for log(gamma^2) */ getline( &ptr, &len, io ); sscanf( ptr, "%le", &lg_gam2_min ); free( ptr ); ptr = NULL; /* read step size in dex */ getline( &ptr, &len, io ); sscanf( ptr, "%le", &step ); free( ptr ); ptr = NULL; lg_gam2_max = lg_gam2_min + (double)(np_gam2-1)*step; /* read atomic number */ getline( &ptr, &len, io ); free( ptr ); ptr = NULL; /* next lines are comments */ for( i=0; i < 2; ++i ) { getline( &ptr, &len, io ); free( ptr ); ptr = NULL; } for( ipg2=0; ipg2 < np_gam2; ++ipg2 ) { fscanf( io, "%le %le", &log_g2[ipg2], &gff[ipg2] ); gff[ipg2] = log(gff[ipg2]); } fclose(io); } double lagrange(const double x[], /* x[n] */ const double y[], /* y[n] */ long n, double xval) { long i, j; double yval = 0.; for( i=0; i < n; i++ ) { double l = 1.; for( j=0; j < n; j++ ) { if( i != j ) l *= (xval-x[j])/(x[i]-x[j]); } yval += y[i]*l; } return yval; } double gauntff(double gam2) { /* use 3rd-order interpolation */ static const size_t NINTERP = 4; long ipg2; double lg_gam2; assert( gam2 > 0. ); lg_gam2 = log10(gam2); assert( lg_gam2_min <= lg_gam2 && lg_gam2 <= lg_gam2_max ); ipg2 = min(max((size_t)((lg_gam2-lg_gam2_min)/step),1)-1,NP_GAM2-NINTERP); return lagrange(&log_g2[ipg2], &gff[ipg2], NINTERP, lg_gam2); } int main() { double gam2; read_table( gauntff_file ); printf( "enter gam2: " ); scanf( "%le", &gam2 ); printf( "gam2 = %e, gauntff = %e\n", gam2, exp(gauntff(gam2)) ); return 0; }