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