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 #define gaunt_magic 20140210 #define gaunt_magic2 20140510 c use 3rd-order interpolation #define NINTERP 4 c uncomment the following line to interpolate on the non-thermally averaged Gaunt factors c #define USE_NONAV_GAUNT 1 #ifdef USE_NONAV_GAUNT #define gauntff_file "gauntff_nonav.dat" #define NP_GAM2 151 #define NP_U 276 #else #define gauntff_file "gauntff.dat" #define NP_GAM2 81 #define NP_U 146 #endif program interp implicit none double precision gam2, u, gauntff, res call read_table( gauntff_file ) #ifdef USE_NONAV_GAUNT print*, "enter eps_i, w: " #else print*, "enter gam2, u: " #endif read*, gam2, u res = gauntff( gam2, u ) print*, "gam2 = ", gam2, ", u = ", u, ", gauntff = ", 1 res end subroutine read_table(fnam) implicit none character fnam*(*) common/gffdat/ gff, lg_gam2_min, lg_u_min, step, lg_gam2_max, 1 lg_u_max double precision gff(NP_GAM2,NP_U), lg_gam2_min, lg_u_min, step, 1 lg_gam2_max, lg_u_max integer magic, np_gam2, np_u, i, ipu, ipg2 character line*2800 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 .and. magic.ne.gaunt_magic2 ) then print*, "read_table found wrong magic number in file", fnam stop endif c read dimensions of the table read(10,*) np_gam2, np_u c read start value for log(gamma^2) read(10,*) lg_gam2_min c read start value for log(u) read(10,*) lg_u_min c read step size in dex read(10,*) step lg_gam2_max = lg_gam2_min + dble(np_gam2-1)*step lg_u_max = lg_u_min + dble(np_u-1)*step c read atomic number when present if( magic.eq.gaunt_magic2 ) then read(10,*) line endif c next lines are comments do i=1,9 read(10,*) line enddo do ipu=1,np_u read(10,*) (gff(ipg2,ipu), ipg2=1,np_gam2) enddo c the remainder of the file contains the uncertainties of the gff values c we will not read those here... close(10) return 999 print*, "failed to open file ", fnam stop end double precision function gauntff(gam2, u) implicit none double precision gam2, u common/gffdat/ gff, lg_gam2_min, lg_u_min, step, lg_gam2_max, 1 lg_u_max double precision gff(NP_GAM2,NP_U), lg_gam2_min, lg_u_min, step, 1 lg_gam2_max, lg_u_max integer i, ipg2, ipu double precision lg_gam2, lg_u, x(NINTERP), interp(NINTERP) double precision lagrange lg_gam2 = dlog10(gam2) lg_u = dlog10(u) if( lg_gam2.lt.lg_gam2_min .or. lg_gam2.gt.lg_gam2_max .or. 1 lg_u.lt.lg_u_min .or. lg_u.gt.lg_u_max ) then stop 'parameter out of range' endif ipg2 = min(max(int((lg_gam2-lg_gam2_min)/step),1)-1, 1 NP_GAM2-NINTERP) ipu = min(max(int((lg_u-lg_u_min)/step),1)-1,NP_U-NINTERP) do i=1,NINTERP x(i) = lg_gam2_min + dble(ipg2+i-1)*step enddo do i=1,NINTERP interp(i) = lagrange(x, gff(ipg2+1,ipu+i), NINTERP, lg_gam2) enddo do i=1,NINTERP x(i) = lg_u_min + dble(ipu+i-1)*step enddo gauntff = lagrange(x, interp, NINTERP, lg_u) 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