/*
    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

    van Hoof P.A.M., Ferland G.J., Williams R.J.R., Volk K., Chatzikos M., Lykins M., Porter R.L.
    Accurate determination of the free-free Gaunt factor, II -- relativistic Gaunt factors
    2015, MNRAS, 449, 2112
*/

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

static const long gaunt_magic = 20141008L;

static const char* gauntff_file = "gauntff_freqint_Z01.dat";

#define NP_GAM2 161

double log_g2[NP_GAM2], gff[NP_GAM2];

double lg_gam2_min, step, lg_gam2_max;

size_t min(size_t x, size_t y)
{
	return ( x < y ) ? x : y;
}

size_t max(size_t x, size_t y)
{
	return ( x > y ) ? x : y;
}

void read_table(const char* fnam)
{
	char* ptr;
	long magic;
	size_t np_gam2, i, ipg2, len;
	FILE* io = fopen( fnam, "r" );
	if( io == NULL )
	{
		fprintf( stderr, "failed to open file %s\n", fnam );
		exit(1);
	}
	ptr = NULL;
	/* skip copyright statement */
	while( 1 )
	{
		getline( &ptr, &len, io );
		if( ptr[0] != '#' )
			break;
		free( ptr );
		ptr = NULL;
	}
	/* read magic number */
	sscanf( ptr, "%ld", &magic );
	free( ptr );
	ptr = NULL;
	if( magic != gaunt_magic )
	{
		fprintf( stderr, "read_table() found wrong magic number in file %s.\n", fnam );
		exit(1);
	}
	/* read dimensions of the table */
	getline( &ptr, &len, io );
	sscanf( ptr, "%ld",  &np_gam2 );
	free( ptr );
	ptr = NULL;
	assert( np_gam2 == NP_GAM2 );
	/* read start value for log(gamma^2) */
	getline( &ptr, &len, io );
	sscanf( ptr, "%le",  &lg_gam2_min );
	free( ptr );
	ptr = NULL;
	/* read step size in dex */
	getline( &ptr, &len, io );
	sscanf( ptr, "%le",  &step );
	free( ptr );
	ptr = NULL;
	lg_gam2_max = lg_gam2_min + (double)(np_gam2-1)*step;
	/* read atomic number */
	getline( &ptr, &len, io );
	free( ptr );
	ptr = NULL;

	/* next lines are comments */
	for( i=0; i < 2; ++i )
	{
		getline( &ptr, &len, io );
		free( ptr );
		ptr = NULL;
	}

	for( ipg2=0; ipg2 < np_gam2; ++ipg2 ) {
		fscanf( io, "%le %le", &log_g2[ipg2], &gff[ipg2] );
		gff[ipg2] = log(gff[ipg2]);
	}

	fclose(io);
}

double lagrange(const double x[], /* x[n] */
		const double y[], /* y[n] */
		long n,
		double xval)
{
	long i, j;
	double yval = 0.;

	for( i=0; i < n; i++ )
	{
		double l = 1.;
		for( j=0; j < n; j++ )
		{
			if( i != j )
				l *= (xval-x[j])/(x[i]-x[j]);
		}
		yval += y[i]*l;
	}
	return yval;
}

double gauntff(double gam2)
{
	/* use 3rd-order interpolation */
	static const size_t NINTERP = 4;

	long ipg2;
	double lg_gam2;

	assert( gam2 > 0. );

	lg_gam2 = log10(gam2);

	assert( lg_gam2_min <= lg_gam2 && lg_gam2 <= lg_gam2_max );

	ipg2 = min(max((size_t)((lg_gam2-lg_gam2_min)/step),1)-1,NP_GAM2-NINTERP);

	return lagrange(&log_g2[ipg2], &gff[ipg2], NINTERP, lg_gam2);
}

int main()
{
	double gam2;
	read_table( gauntff_file );
	printf( "enter gam2: " );
	scanf( "%le", &gam2 );
	printf( "gam2 = %e, gauntff = %e\n", gam2, exp(gauntff(gam2)) );
	return 0;
}