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