/*tababun interpolate on table of points to do 'element table' command, */
#include <float.h>
#include "cddefines.h"
#include "abuntablfac.h"
#include "abuntabllg.h"
#include "tababun.h"

double tababun(double r0, 
  double depth, 
  long int iel)
{
	int lgHit;
	long int j;
	double frac, 
	  tababun_v, 
	  x;

#	ifdef DEBUG_FUN
	fputs( "<+>tababun()\n", debug_fp );
#	endif
	/*interpolate on table of points to do 'element table' command,
	 *based on code by K Volk
	 *each line is log radius and abundance. */
	/* atomic number */

	/* interpolate on radius or depth? */
	if( AbunTablLG.lgAbTaDepth[iel-1] )
	{
		/* depth key appeared = we want depth */
		x = log10(depth);
	}
	else
	{
		/* use radius */
		x = log10(r0);
	}

	/* this will be reset below, but is here as a safety check */
	tababun_v = -DBL_MAX;

	if( x < AbunTablFac.AbTabRad[0][iel-1] || x >= AbunTablFac.AbTabRad[AbunTablFac.nAbunTabl-1][iel-1] )
	{
		fprintf( ioQQQ, " requested radius outside range of tababun\n" );
		fprintf( ioQQQ, " radius was%10.2e min, max=%10.2e%10.2e\n", 
		  x, AbunTablFac.AbTabRad[0][iel-1], AbunTablFac.AbTabRad[AbunTablFac.nAbunTabl-1][iel-1] );
		puts( "[Stop in tababun]" );
		exit(1);
	}

	else
	{
		lgHit = FALSE;
		j = 1;

		while( !lgHit && j <= AbunTablFac.nAbunTabl - 1 )
		{
			if( AbunTablFac.AbTabRad[j-1][iel-1] <= x && 
				AbunTablFac.AbTabRad[j][iel-1] > x )
			{
				frac = (x - AbunTablFac.AbTabRad[j-1][iel-1])/(AbunTablFac.AbTabRad[j][iel-1] - 
				  AbunTablFac.AbTabRad[j-1][iel-1]);
				tababun_v = AbunTablFac.AbTabFac[j-1][iel-1] + frac*
				  (AbunTablFac.AbTabFac[j][iel-1] - AbunTablFac.AbTabFac[j-1][iel-1]);
				lgHit = TRUE;
			}
			j += 1;
		}

		if( !lgHit )
		{
			fprintf( ioQQQ, " radius outran dlaw table scale, requested=%6.2f largest=%6.2f\n", 
			  x, AbunTablFac.AbTabRad[AbunTablFac.nAbunTabl-1][iel-1] );
			puts( "[Stop in tababun]" );
			exit(1);
		}
	}

	/* got it, now return value, not log of density */
	tababun_v = pow(10.,tababun_v);


#	ifdef DEBUG_FUN
	fputs( " <->tababun()\n", debug_fp );
#	endif
	return( tababun_v );
}