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

	program interp
	implicit none
	double precision gam2, res, total_gff
	print*, "enter gamma^2: "
	read*, gam2
	res = total_gff(gam2)
	print*, "<g_ff(gamma^2)>: ", res
	end

	double precision function total_gff(gam2)
	implicit none
	double precision gam2
	double precision g, P, Q
	double precision al(5) / 1.43251926625281d+00,
	1    3.50626935257777d-01, 4.36183448595035d-01,
	2    6.03536387105599d-02, 3.66626405363100d-02 /
	double precision bl(5) / 1.00000000000000d+00,
	1    2.92525161994346d-01, 4.05566949766954d-01,
	2    5.62573012783879d-02, 3.33019373823972d-02 /
	double precision ah(5) / 1.45481634667278d+00,
	1    -9.55399384620923d-02, 1.46327814151538d-01,
	2    -1.41489406498468d-02, 2.76891413242655d-03 /
	double precision bh(5) / 1.00000000000000d+00,
	1    3.31149751183539d-02, 1.31127367293310d-01,
	2    -1.32658217746618d-02, 2.74809263365693d-03 /

	if( gam2.lt.1d-6 ) then
	   total_gff = 1.102635d0+1.186d0*sqrt(gam2)+0.86d0*gam2
	   return
	else if( gam2.gt.1d10 ) then
	   total_gff = 1.d0+gam2**(-1.d0/3.d0)
	   return
	else
	   g = log10(gam2)
	   if( g.le.0.8d0 ) then
	      P = (((al(5)*g+al(4))*g+al(3))*g+al(2))*g+al(1)
	      Q = (((bl(5)*g+bl(4))*g+bl(3))*g+bl(2))*g+bl(1)
	   else
	      P = (((ah(5)*g+ah(4))*g+ah(3))*g+ah(2))*g+ah(1)
	      Q = (((bh(5)*g+bh(4))*g+bh(3))*g+bh(2))*g+bh(1)
	   endif
	   total_gff = P/Q
	endif
	return
	end