/*
    Copyright (c) 2014, Peter A.M. van Hoof.
  
    This program is provided 'as-is', without any express or implied warranty. In
    no event will the author be held liable for any damages arising from the use
    of this program.
  
    Permission is granted to anyone to use this program for any purpose, including
    commercial applications, and to alter it and redistribute it freely, subject
    to the following restrictions:
  
    1. The origin of this program must not be misrepresented; you must not claim
       that you created the original program. If you use this program in a product,
       an acknowledgment in the product documentation would be appreciated but
       is not required.
    2. Altered program versions must be plainly marked as such, and must not be
       misrepresented as being the original program.
    3. This notice may not be removed or altered from any further distribution.
  
    Peter A.M. van Hoof
    Royal Observatory of Belgium
    Ringlaan 3
    B-1180 Brussels
    Belgium
    p.vanhoof@oma.be
*/

/*
    If you use any of these data in a scientific publication, please refer to

    van Hoof P.A.M., Williams R.J.R., Volk K., Chatzikos M., Ferland G.J., Lykins M., Porter R.L., Wang Y.
    Accurate determination of the free-free Gaunt factor, I -- non-relativistic Gaunt factors
    2014, MNRAS, 444, 420
*/

#include <stdio.h>
#include <math.h>

double total_gff(double);

int main()
{
	double gam2;
	printf( "enter gamma^2: " );
	scanf( "%le", &gam2 );
	printf( "<g_ff(gamma^2)>: %.6f\n", total_gff(gam2) );
	return 0;
}

double total_gff(double gam2)
{
	static const double al[5] = { 1.43251926625281e+00, 3.50626935257777e-01, 4.36183448595035e-01,
				      6.03536387105599e-02, 3.66626405363100e-02 };
	static const double bl[5] = { 1.00000000000000e+00, 2.92525161994346e-01, 4.05566949766954e-01,
				      5.62573012783879e-02, 3.33019373823972e-02 };
	static const double ah[5] = { 1.45481634667278e+00, -9.55399384620923e-02, 1.46327814151538e-01,
				      -1.41489406498468e-02, 2.76891413242655e-03 };
	static const double bh[5] = { 1.00000000000000e+00, 3.31149751183539e-02, 1.31127367293310e-01,
				      -1.32658217746618e-02, 2.74809263365693e-03 };

	if( gam2 < 1e-6 ) {
		return 1.102635 + 1.186*sqrt(gam2) + 0.86*gam2;
	}
	else if( gam2 > 1e10 ) {
		return 1. + 1./cbrt(gam2);
	}
	else {
		double g = log10(gam2);
		double P, Q;
		if( g <= 0.8 ) {
			P = (((al[4]*g + al[3])*g + al[2])*g + al[1])*g + al[0];
			Q = (((bl[4]*g + bl[3])*g + bl[2])*g + bl[1])*g + bl[0];
		}
		else {
			P = (((ah[4]*g + ah[3])*g + ah[2])*g + ah[1])*g + ah[0];
			Q = (((bh[4]*g + bh[3])*g + bh[2])*g + bh[1])*g + bh[0];
		}
		return P/Q;
	}
}