c c Copyright (c) 2014, Peter A.M. van Hoof. c c This program is provided 'as-is', without any express or implied warranty. In c no event will the author be held liable for any damages arising from the use c of this program. c c Permission is granted to anyone to use this program for any purpose, including c commercial applications, and to alter it and redistribute it freely, subject c to the following restrictions: c c 1. The origin of this program must not be misrepresented; you must not claim c that you created the original program. If you use this program in a product, c an acknowledgment in the product documentation would be appreciated but c is not required. c 2. Altered program versions must be plainly marked as such, and must not be c misrepresented as being the original program. c 3. This notice may not be removed or altered from any further distribution. c c Peter A.M. van Hoof c Royal Observatory of Belgium c Ringlaan 3 c B-1180 Brussels c Belgium c p.vanhoof@oma.be c c c If you use any of these data in a scientific publication, please refer to c c van Hoof P.A.M., Williams R.J.R., Volk K., Chatzikos M., Ferland G.J., Lykins M., Porter R.L., Wang Y. c Accurate determination of the free-free Gaunt factor, I -- non-relativistic Gaunt factors c 2014, MNRAS, 444, 420 c c van Hoof P.A.M., Ferland G.J., Williams R.J.R., Volk K., Chatzikos M., Lykins M., Porter R.L. c Accurate determination of the free-free Gaunt factor, II -- relativistic Gaunt factors c 2015, MNRAS, 449, 2112 c #define gaunt_magic 20141008 c use 3rd-order interpolation #define NINTERP 4 #define gauntff_file "gauntff_freqint_Z01.dat" #define NP_GAM2 161 program interp implicit none double precision gam2, gauntff, res call read_table( gauntff_file ) print*, "enter gam2: " read*, gam2 res = exp( gauntff( gam2 ) ) print*, "gam2 = ", gam2, ", gauntff = ", res end subroutine read_table(fnam) implicit none character fnam*(*) common/gffdat/ log_g2, gff, lg_gam2_min, step, lg_gam2_max double precision log_g2(NP_GAM2), gff(NP_GAM2), lg_gam2_min, 1 step, lg_gam2_max integer magic, np_gam2, i, ipg2 character line*100 open(unit=10, file=fnam, status="old", err=999) c skip copyright statement do i=1,28 read(10,*) line enddo c read magic number read(10,*) magic if( magic.ne.gaunt_magic ) then print*, "read_table found wrong magic number in file", fnam stop endif c read dimensions of the table read(10,*) np_gam2 c read start value for log(gamma^2) read(10,*) lg_gam2_min c read step size in dex read(10,*) step lg_gam2_max = lg_gam2_min + dble(np_gam2-1)*step c read atomic number when present read(10,*) line c next lines are comments do i=1,2 read(10,*) line enddo do ipg2=1,np_gam2 read(10,*) log_g2(ipg2), gff(ipg2) gff(ipg2) = log(gff(ipg2)) enddo close(10) return 999 print*, "failed to open file ", fnam stop end double precision function gauntff(gam2) implicit none double precision gam2 common/gffdat/ log_g2, gff, lg_gam2_min, step, lg_gam2_max double precision log_g2(NP_GAM2), gff(NP_GAM2), lg_gam2_min, 1 step, lg_gam2_max integer ipg2 double precision lg_gam2, lagrange lg_gam2 = dlog10(gam2) if( lg_gam2.lt.lg_gam2_min .or. lg_gam2.gt.lg_gam2_max ) then stop 'parameter out of range' endif ipg2 = min(max(int((lg_gam2-lg_gam2_min)/step),1)-1, 1 NP_GAM2-NINTERP) gauntff = lagrange(log_g2(ipg2+1), gff(ipg2+1), NINTERP, lg_gam2) return end double precision function lagrange(x, y, n, xval) implicit none integer n double precision x(n), y(n), xval integer i, j double precision l lagrange = 0.d0 do i=1,n l = 1.d0 do j=1,n if( i.ne.j ) then l = l*(xval-x(j))/(x(i)-x(j)) endif enddo lagrange = lagrange + y(i)*l enddo return end