/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*atmdat_3body derive three-body recombination coefficients */
/*da interpolate on three body recombination by Steve Cota */
#include "cddefines.h"
#include "ionbal.h"
#include "phycon.h"
#include "dense.h"
#include "trace.h"
#include "save.h"
#include "atmdat.h"

#define	MAXZ	28

STATIC void blkdata1(void);
STATIC double da(double z, double temp, double eden);

static double a2[63], 
	  b2[63], 
	  x2[63];

static double a0[83], 
	  x0[83];
static realnum b0[83],
	  b1[83]; 

static double a1[83], 
	  x1[83];

static double tz[83], 
	  zlog7[28], 
	  zlog2[28];

#define RC_INI(rs)	(rs[_r].rc>1 ? DEC_RC_(rs) : (rs[_r].rc==1 ? INC_NDX_(rs) : rs[_r].ini ))
#define DEC_RC_(rs)	(rs[_r].rc--,rs[_r].ini)
#define INC_NDX_(rs)	(_r++,rs[_r-1].ini)

/* this "mapping function" occurs below */
STATIC double xmap(double x[], 
	  double y[], 
	  double x0);

/* inverse routine, also below */
STATIC double xinvrs(double y, 
	  double a, 
	  double b, 
	  double u, 
	  double v, 
	  long int *ifail);

/* =================================================================== */
void atmdat_3body(void)
{
	long int i, 
	  iup;

	DEBUG_ENTRY( "atmdat_3body()" );


	if( ionbal.lgNoCota )
	{
		for( i=0; i < LIMELM; i++ )
		{
			ionbal.CotaRate[i] = 0.;
		}
		atmdat.nsbig = 0;
		return;
	}

	if( atmdat.nsbig == 0 )
	{
		/* steve cota only defined things up to 28 */
		iup = MIN2(28,LIMELM);
	}
	else
	{
		iup = MIN3( LIMELM , atmdat.nsbig , 28 ); 
	}

	for( i=0; i < iup; i++ )
	{
		ionbal.CotaRate[i] = (realnum)da((double)(i+1), phycon.te, dense.eden);
	}

	atmdat.nsbig = 0;

	if( trace.lgTrace && trace.lgTrace3Bod )
	{
		fprintf( ioQQQ, "     3BOD rate:" );
		for( i=1; i <= 14; i++ )
		{
			fprintf( ioQQQ, "%8.1e", ionbal.CotaRate[i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

	if( save.lgioRecom )
	{
		/* option to save coefficients */
		fprintf( save.ioRecom, " 3-body rec coef vs charge \n" );
		for( i=0; i < iup; i++ )
		{
			fprintf( save.ioRecom, "%3ld%10.2e\n", i+1, ionbal.CotaRate[i] );
		}
		fprintf( save.ioRecom, "\n");
	}
	return;
}

/* =================================================================== */
STATIC double da(double z, double temp, double eden)
{
	/*lint -e736 loss of precision in assignment in translated data */
	long int jfail, 
	  nt, 
	  nt0, 
	  nt1, 
	  nz;

	static bool lgCalled=false;
	double a, 
	  alogn, 
	  alognc, 
	  alogt, 
	  b, 
	  c, 
	  d, 
	  da_v, 
	  expp, 
	  u, 
	  v, 
	  x, 
	  xnc, 
	  y, 
	  zlog;
	double yya[3], 
	  xx[3], 
	  yyb[3], 
	  yyx[3], 
	  yyy[3];

	/* alogte is base 10 log of temperature and elec density*/
	double alogte , alogne;

	DEBUG_ENTRY( "da()" );

	/* WRITTEN BY S. A. COTA, 2/87
	 * */

	/* MAXZ IS THE MAXIMUM EFFECTIVE NUCLEAR CHARGE ( = IONIC CHARGE + 1 )
	 * WHICH THE DIMENSION STATEMENTS ACCOMODATE.  
	 *
	 * IT IS USED ONLY FOR THE ARRAY ZLOG7 ( = 7 * LOG ( Z ) ) 
	 * AND THE ARRAY ZLOG2 ( = 2 * LOG ( Z ) ) .   THESE ARRAYS
	 * CONTAIN EASILY CALCULATED VALUES, WHICH HAVE BEEN STORED
	 * TO SAVE TIME IN EXECUTION.
	 *
	 * IF MAXZ IS EXCEEDED, THIS PROGRAM SIMPLY CALCULATES THE
	 * LOGS INSTEAD OF LOOKING THEM UP.
	 *
	 * */

	if( !lgCalled )
	{
		lgCalled = true;
		blkdata1();
	}

	/*begin sanity check */
	ASSERT( zlog7[1] > 0. );

	nz = (long)(z + .1);
	// The original routine has major discontinuities
	// The use of eden_limited (here and at the end of the routine) is an attempt 
	// to smooth these problems and still get roughly the same behavior.
	double eden_limited = MIN2( eden, 1e5 );
	alogne = log10(eden_limited);
	alogte = log10(temp);
	if( nz > MAXZ )
	{
		zlog = log10(z);
		alogt = alogte - 2.*zlog;
		alogn = alogne - 7.*zlog;
	}
	else
	{
		alogt = alogte - zlog2[nz-1];
		alogn = alogne - zlog7[nz-1];
	}

	/* CHECK IF PARAMETERS ARE WITHIN BOUNDS.  IF NOT, INCREMENT
	 * APPROPRIATE ERROR COUNTER AND SET TO BOUNDARY IF 
	 * NECESSARY:
	 *
	 * DEFINITION OF ERROR COUNTERS:
	 *
	 * ILT    : LOW T         
	 * ILTLN  : LOW T , LOW  N
	 * ILTHN  : LOW T , HIGH N
	 * IHTHN  : HIGH T , HIGH N 
	 * */
	if( alogt < 0. )
	{
		ionbal.ilt += 1;
		alogt = 0.;
	}

	if( alogt <= 2.1760913 )
	{
		if( alogn < (3.5*alogt - 8.) )
		{
			ionbal.iltln += 1;
		}
		else if( alogn > (3.5*alogt - 2.) )
		{
			ionbal.ilthn += 1;
			alogn = 3.5*alogt - 2.;
		}

	}
	else if( alogt <= 2.4771213 )
	{
		if( alogn > 9. )
		{
			ASSERT( 0 );
			ionbal.ilthn += 1;
			alogn = 9.;
		}

	}
	else if( alogt <= 5.1139434 )
	{
		if( alogn > 13. )
		{
			ASSERT( 0 );
			ionbal.ihthn += 1;
			alogn = 13.;
		}

	}
	
	// limit used temperature to 10^5 K and 
	// extrapolate with T^-2 at end of routine (to avoid discontinuous behavior).
	double alogt_limited = MIN2( alogt, 5.0 );
	double alogt_save = alogt;
	alogt = alogt_limited;

	/* LOCATE POSITION IN ARRAYS */
	if( alogt <= 2. )
	{
		nt = (long)(9.9657843*alogt + 1.);
	}
	else
	{
		nt = (long)(19.931568*alogt - 19.);
	}
	nt = MIN2(83,nt);
	nt = MAX2(1,nt);

	/* CENTER UP SINCE ARRAY VALUES ARE ROUNDED */
	if( fabs(alogt-tz[nt-1]) >= fabs(alogt-tz[MIN2(83,nt+1)-1]) )
	{
		nt = MIN2(83,nt+1);
	}
	else if( fabs(alogt-tz[nt-1]) > fabs(alogt-tz[MAX2(1,nt-1)-1]) )
	{
		nt = MAX2(1,nt-1);
	}

	/* CHECK IF INTERPOLATION IS NEEDED AND PROCEED IF NOT.*/
	if( fabs(alogt-tz[nt-1]) < 0.00001 )
	{
		if( z != 1.0 )
		{
			c = a1[nt-1];
			d = b1[nt-1];
			u = x1[nt-1];
			v = 8.90;
		}
		else
		{
			nt = MAX2(21,nt);
			nt = MIN2(83,nt);
			c = a2[nt-(21)];
			d = b2[nt-(21)];
			u = x2[nt-(21)];
			v = 9.40;
		}

		xnc = xinvrs(alogn,c,d,u,v,&jfail);
		if( xnc <= 0. || jfail != 0 )
		{
			ionbal.ifail = 1;
			jfail = 1;
			da_v = 0.;
			return( da_v );
		}
		alognc = log10(xnc);

		a = a0[nt-1];
		b = b0[nt-1];
		x = -2.45;
		y = x0[nt-1];
		nt0 = nt - 1;

		/* IF INTERPOLATION WAS REQUIRED, 
		 * CHECK THAT NT IS NOT ON THE EDGE OF A DISCONTINUITY,
		 * EITHER AT END OF ARRAYS OR WITHIN THEM,
		 * WHERE VALUES CHANGE ABRUPTLY.
		 * */
	}
	else
	{
		if( (nt <= 21) && (z == 1.0) )
		{
			nt = 22;
		}
		else if( nt <= 1 )
		{
			nt = 2;
		}
		else if( nt >= 83 )
		{
			nt = 82;
		}
		else if( nt == 24 )
		{
			if( alogt <= 2.1760913 )
			{
				nt = 23;
			}
			else
			{
				nt = 25;
			}
		}
		else if( nt == 30 )
		{
			if( alogt <= 2.4771213 )
			{
				nt = 29;
			}
			else
			{
				nt = 31;
			}
		}

		nt0 = nt - 1;
		nt1 = nt + 1;
		xx[0] = tz[nt0-1];
		xx[1] = tz[nt-1];
		xx[2] = tz[nt1-1];

		if( z != 1.0 )
		{
			if( nt0 == 24 )
			{
				yya[0] = 17.2880135;
				yyb[0] = 6.93410742e03;
				yyx[0] = -3.75;
			}
			else if( nt0 == 30 )
			{
				yya[0] = 17.4317989;
				yyb[0] = 1.36653821e03;
				yyx[0] = -3.40;
			}
			else
			{
				yya[0] = a1[nt0-1];
				yyb[0] = b1[nt0-1];
				yyx[0] = x1[nt0-1];
			}

			yya[1] = a1[nt-1];
			yya[2] = a1[nt1-1];
			c = xmap(xx,yya,alogt);
			yyb[1] = b1[nt-1];
			yyb[2] = b1[nt1-1];
			d = xmap(xx,yyb,alogt);
			yyx[1] = x1[nt-1];
			yyx[2] = x1[nt1-1];
			u = xmap(xx,yyx,alogt);
			v = 8.90;

		}
		else
		{
			if( nt0 == 24 )
			{
				yya[0] = 20.1895161;
				yyb[0] = 2.25774918e01;
				yyx[0] = -1.66;
			}
			else if( nt0 == 30 )
			{
				yya[0] = 19.8647804;
				yyb[0] = 6.70408707e02;
				yyx[0] = -2.12;
			}
			else
			{
				yya[0] = a2[nt0-(21)];
				yyb[0] = b2[nt0-(21)];
				yyx[0] = x2[nt0-(21)];
			}

			yya[1] = a2[nt-(21)];
			yya[2] = a2[nt1-(21)];
			c = xmap(xx,yya,alogt);
			yyb[1] = b2[nt-(21)];
			yyb[2] = b2[nt1-(21)];
			d = xmap(xx,yyb,alogt);
			yyx[1] = x2[nt-(21)];
			yyx[2] = x2[nt1-(21)];
			u = xmap(xx,yyx,alogt);
			v = 9.40;
		}

		xnc = xinvrs(alogn,c,d,u,v,&jfail);
		if( xnc <= 0. || jfail != 0 )
		{
			ionbal.ifail = 1;
			jfail = 1;
			da_v = 0.;
			return( da_v );
		}
		alognc = log10(xnc);

		if( nt0 == 24 )
		{
			yya[0] = -8.04963875;
			yyb[0] = 1.07205127e03;
			yyy[0] = 2.05;
		}
		else if( nt0 == 30 )
		{
			yya[0] = -8.54721069;
			yyb[0] = 4.70450195e02;
			yyy[0] = 2.05;
		}
		else
		{
			yya[0] = a0[nt0-1];
			yyb[0] = b0[nt0-1];
			yyy[0] = x0[nt0-1];
		}

		yya[1] = a0[nt-1];
		yya[2] = a0[nt1-1];
		a = xmap(xx,yya,alogt);
		yyb[1] = b0[nt-1];
		yyb[2] = b0[nt1-1];
		b = xmap(xx,yyb,alogt);
		x = -2.45;
		yyy[1] = x0[nt-1];
		yyy[2] = x0[nt1-1];
		y = xmap(xx,yyy,alogt);
	}

	expp = a - y*alognc + b*pow(xnc,x);
	if( expp < 37 )
	{
		da_v = z*pow(10.,expp);
	}
	else
	{
		da_v = 0.;
	}
	ionbal.ifail += jfail;

	da_v *= pow( eden/eden_limited, 0.25 );
	da_v *= pow( 10., 2.*(alogt_limited-alogt_save) );

	return( da_v );
}

/****************************************************************************** */
STATIC void blkdata1(void)
{
	/*block data with Steve Cota's 3-body recombination coefficients */

	/* data for function  da.
	 *
	 * S. A. COTA, 2/1987
	 * */

	long int _i, 
	  _r;
	realnum *const ba0 = (realnum*)b0;
	realnum *const ba1 = (realnum*)b1;
	realnum *const bb0 = (realnum*)((char*)(b0 + 79));
	realnum *const bb1 = (realnum*)((char*)(b1 + 79));

		/* to fix all the conversion errors, change realnum ini to double ini,
		 * but chech that results still ok */
		{ static struct{ long rc; double ini; } _rs0[] = {
			{1,	0.00000e00},
			{1,	2.10721e00},
			{1,	3.33985e00},
			{1,	4.21442e00},
			{1,	4.89279e00},
			{1,	5.44706e00},
			{1,	5.91569e00},
			{1,	6.32163e00},
			{1,	6.67970e00},
			{1,	7.00000e00},
			{1,	7.28975e00},
			{1,	7.55427e00},
			{1,	7.79760e00},
			{1,	8.02290e00},
			{1,	8.23264e00},
			{1,	8.42884e00},
			{1,	8.61314e00},
			{1,	8.78691e00},
			{1,	8.95128e00},
			{1,	9.10721e00},
			{1,	9.25554e00},
			{1,	9.39696e00},
			{1,	9.53209e00},
			{1,	9.66148e00},
			{1,	9.78558e00},
			{1,	9.90481e00},
			{1,	10.01954e00},
			{1,	10.13010e00},
			{0L, 0}
		};
		for(_i=_r=0L; _i < 28; _i++)
			zlog7[_i] = RC_INI(_rs0); }
		{ static struct{ long rc; double ini; } _rs1[] = {
			{1, 	0.00000e00},
			{1, 	6.02060e-01},
			{1, 	9.54243e-01},
			{1, 	1.20412e00},
			{1, 	1.39794e00},
			{1, 	1.55630e00},
			{1, 	1.69020e00},
			{1, 	1.80618e00},
			{1, 	1.90849e00},
			{1, 	2.00000e00},
			{1, 	2.08279e00},
			{1, 	2.15836e00},
			{1, 	2.22789e00},
			{1, 	2.29226e00},
			{1, 	2.35218e00},
			{1, 	2.40824e00},
			{1, 	2.46090e00},
			{1, 	2.51055e00},
			{1, 	2.55751e00},
			{1, 	2.60206e00},
			{1, 	2.64444e00},
			{1, 	2.68485e00},
			{1, 	2.72346e00},
			{1, 	2.76042e00},
			{1, 	2.79588e00},
			{1, 	2.82995e00},
			{1, 	2.86272e00},
			{1, 	2.89431e00},
			{0L, 0}
		};
		for(_i=_r=0L; _i < 28; _i++)
			zlog2[_i] = RC_INI(_rs1); }
		{ static struct{ long rc; double ini; } _rs2[] = {
			{1,	0.},
			{1,	0.09691},
			{1,	0.17609},
			{1,	0.30103},
			{1,	0.39794},
			{1,	0.47712},
			{1,	0.60206},
			{1,	0.69897},
			{1,	0.77815},
			{1,	0.90309},
			{1,	1.00000},
			{1,	1.07918},
			{1,	1.20412},
			{1,	1.30103},
			{1,	1.39794},
			{1,	1.47712},
			{1,	1.60206},
			{1,	1.69897},
			{1,	1.77815},
			{1,	1.90309},
			{1,	2.00000},
			{1,	2.06070},
			{1,	2.09691},
			{1,	2.17609},
			{1,	2.20412},
			{1,	2.24304},
			{1,	2.30103},
			{1,	2.35218},
			{1,	2.39794},
			{1,	2.47712},
			{1,	2.51188},
			{1,	2.54407},
			{1,	2.60206},
			{1,	2.65321},
			{1,	2.69897},
			{1,	2.75967},
			{1,	2.81291},
			{1,	2.86034},
			{1,	2.91645},
			{1,	2.95424},
			{1,	3.00000},
			{1,	3.07918},
			{1,	3.11394},
			{1,	3.17609},
			{1,	3.20412},
			{1,	3.25527},
			{1,	3.30103},
			{1,	3.36173},
			{1,	3.39794},
			{1,	3.46240},
			{1,	3.51188},
			{1,	3.56820},
			{1,	3.60206},
			{1,	3.66276},
			{1,	3.72016},
			{1,	3.76343},
			{1,	3.81291},
			{1,	3.86034},
			{1,	3.90309},
			{1,	3.95424},
			{1,	4.02119},
			{1,	4.06070},
			{1,	4.11394},
			{1,	4.16137},
			{1,	4.20412},
			{1,	4.25527},
			{1,	4.31175},
			{1,	4.36173},
			{1,	4.41497},
			{1,	4.46240},
			{1,	4.51521},
			{1,	4.56526},
			{1,	4.61542},
			{1,	4.66605},
			{1,	4.71600},
			{1,	4.76343},
			{1,	4.81624},
			{1,	4.86629},
			{1,	4.91645},
			{1,	4.96614},
			{1,	5.02119},
			{1,	5.06726},
			{1,	5.11394},
			{0L, 0 }
			};
		for(_i=_r=0L; _i < 83; _i++)
			tz[_i] = RC_INI(_rs2); }
		{ static struct{ long rc; double ini; } _rs3[] = {
			{1,	-4.31396484},
			{1,	-4.56640625},
			{1,	-4.74560547},
			{1,	-4.98535156},
			{1,	-5.15373850},
			{1,	-5.28123093},
			{1,	-5.48215008},
			{1,	-5.63811255},
			{1,	-5.76573515},
			{1,	-5.96755028},
			{1,	-6.12449837},
			{1,	-6.25304174},
			{1,	-6.45615673},
			{1,	-6.61384058},
			{1,	-6.77161551},
			{1,	-6.90069818},
			{1,	-7.10470295},
			{1,	-7.26322412},
			{1,	-7.39289951},
			{1,	-7.59792519},
			{1,	-7.75725508},
			{1,	-7.85722494},
			{1,	-7.91697407},
			{1,	-8.04758644},
			{1,	-8.09447479},
			{1,	-8.15859795},
			{1,	-8.25424385},
			{1,	-8.33880615},
			{1,	-8.41452408},
			{1,	-8.54581165},
			{1,	-8.60400581},
			{1,	-8.65751839},
			{1,	-8.75414848},
			{1,	-8.83946800},
			{1,	-8.91589737},
			{1,	-9.01741695},
			{1,	-9.10663033},
			{1,	-9.18621922},
			{1,	-9.28059292},
			{1,	-9.34430218},
			{1,	-9.42154408},
			{1,	-9.55562973},
			{1,	-9.61459446},
			{1,	-9.72023010},
			{1,	-9.76802444},
			{1,	-9.85540199},
			{1,	-9.93374062},
			{1,	-10.03800774},
			{1,	-10.10044670},
			{1,	-10.21178055},
			{1,	-10.29757786},
			{1,	-10.39561272},
			{1,	-10.45469666},
			{1,	-10.56102180},
			{1,	-10.66205502},
			{1,	-10.73780537},
			{1,	-10.82557774},
			{1,	-10.91007328},
			{1,	-10.98659325},
			{1,	-11.07857418},
			{1,	-11.19975281},
			{1,	-11.27170753},
			{1,	-11.36930943},
			{1,	-11.45675945},
			{1,	-11.53620148},
			{1,	-11.63198853},
			{1,	-11.73875237},
			{1,	-11.83400822},
			{1,	-11.93677044},
			{1,	-12.02933311},
			{1,	-12.13374519},
			{1,	-12.23410702},
			{1,	-12.33664989},
			{1,	-12.44163322},
			{1,	-12.54730415},
			{1,	-12.64975739},
			{1,	-12.76682186},
			{1,	-12.88185978},
			{1,	-13.00052643},
			{1,	-13.12289810},
			{1,	-13.26689529},
			{1,	-13.39390945},
			{1,	-30.00000000},
			{0L, 0 }
			};
		for(_i=_r=0L; _i < 83; _i++)
			a0[_i] = RC_INI(_rs3); }
		{ static struct{ long rc; double ini; } _rs4[] = {
			{1,	4.53776000e05},
			{1,	3.48304000e05},
			{1,	2.80224000e05},
			{1,	1.98128000e05},
			{1,	1.51219797e05},
			{1,	1.21113266e05},
			{1,	8.52812109e04},
			{1,	6.49598125e04},
			{1,	5.20075781e04},
			{1,	3.66190977e04},
			{1,	2.79060723e04},
			{1,	2.23634102e04},
			{1,	1.57683135e04},
			{1,	1.20284307e04},
			{1,	9.17755273e03},
			{1,	7.36044873e03},
			{1,	5.19871680e03},
			{1,	3.97240796e03},
			{1,	3.18934326e03},
			{1,	2.25737622e03},
			{1,	1.72767114e03},
			{1,	1.46202722e03},
			{1,	1.32456628e03},
			{1,	1.06499792e03},
			{1,	9.92735291e02},
			{1,	8.91604858e02},
			{1,	7.59411560e02},
			{1,	6.59120056e02},
			{1,	5.80688965e02},
			{1,	4.66602264e02},
			{1,	4.27612854e02},
			{1,	3.91531494e02},
			{1,	3.34516968e02},
			{1,	2.91021820e02},
			{1,	2.56853912e02},
			{1,	2.17598007e02},
			{1,	1.88145462e02},
			{1,	1.65329865e02},
			{1,	1.41960342e02},
			{1,	1.28181686e02},
			{1,	1.13336761e02},
			{1,	9.17785034e01},
			{1,	8.36242981e01},
			{1,	7.08843536e01},
			{1,	6.58346100e01},
			{1,	5.75790634e01},
			{1,	5.11293755e01},
			{1,	4.37563019e01},
			{1,	3.99226875e01},
			{1,	3.39562836e01},
			{1,	3.00413170e01},
			{1,	2.61871891e01},
			{1,	2.41310368e01},
			{1,	2.08853607e01},
			{1,	1.82632275e01},
			{1,	1.60007000e01},
			{1,	1.42617064e01},
			{1,	1.27951088e01},
			{1,	1.16221066e01},
			{1,	1.03779335e01},
			{1,	8.97864914e00},
			{1,	8.25593281e00},
			{1,	7.39339924e00},
			{1,	6.70784378e00},
			{1,	6.16084862e00},
			{1,	5.57818031e00},
			{1,	5.01341105e00},
			{1,	4.55679178e00},
			{1,	4.13692093e00},
			{1,	3.80004382e00},
			{1,	3.46328306e00},
			{1,	3.17340493e00},
			{1,	2.93525696e00},
			{1,	2.69083858e00},
			{1,	2.46588683e00},
			{1,	2.26083040e00},
			{1,	2.04337358e00},
			{1,	1.89027369e00},
			{1,	1.69208312e00},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 79; _i++)
			ba0[_i] = (realnum)RC_INI(_rs4); }
		{ static struct{ long rc; double ini; } _rs5[] = {
			{1,	1.48992336e00},
			{1,	1.32466662e00},
			{1,	1.10697949e00},
			{1,	9.29813862e-01},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 4; _i++)
			bb0[_i] = (realnum)RC_INI(_rs5); }
		{ static struct{ long rc; double ini; } _rs6[] = {
			{1,	2.12597656},
			{1,	2.08984375},
			{1,	2.06958008},
			{1,	2.05444336},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{1,	2.05},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 83; _i++)
			x0[_i] = RC_INI(_rs6); }

		{ static struct{ long rc; double ini; } _rs7[] = {
			{1,	16.23337936},
			{1,	16.27946854},
			{1,	16.31696320},
			{1,	16.37597656},
			{1,	16.42210960},
			{1,	16.45996284},
			{1,	16.51994896},
			{1,	16.56644440},
			{1,	16.60460854},
			{1,	16.66510773},
			{1,	16.71198654},
			{1,	16.75038719},
			{1,	16.81106949},
			{1,	16.85778809},
			{1,	16.90416527},
			{1,	16.94209099},
			{1,	17.00195694},
			{1,	17.04838943},
			{1,	17.08633804},
			{1,	17.14627838},
			{1,	17.19270515},
			{1,	17.22186279},
			{1,	17.23933601},
			{1,	17.27728271},
			{1,	17.30161858},
			{1,	17.32085800},
			{1,	17.34928894},
			{1,	17.37349129},
			{1,	17.39528084},
			{1,	17.43282318},
			{1,	17.44827652},
			{1,	17.46357536},
			{1,	17.49082375},
			{1,	17.51517677},
			{1,	17.53697205},
			{1,	17.56587219},
			{1,	17.59125519},
			{1,	17.61410332},
			{1,	17.64081383},
			{1,	17.65900803},
			{1,	17.68086433},
			{1,	17.71843529},
			{1,	17.73512840},
			{1,	17.76512146},
			{1,	17.77873421},
			{1,	17.80340767},
			{1,	17.82530022},
			{1,	17.85470963},
			{1,	17.87210464},
			{1,	17.90334511},
			{1,	17.92751503},
			{1,	17.95458603},
			{1,	17.97117233},
			{1,	18.00062943},
			{1,	18.02842712},
			{1,	18.04934502},
			{1,	18.07340050},
			{1,	18.09639168},
			{1,	18.11732864},
			{1,	18.14218903},
			{1,	18.17465591},
			{1,	18.19370079},
			{1,	18.21962166},
			{1,	18.24237251},
			{1,	18.26305962},
			{1,	18.28767967},
			{1,	18.31531525},
			{1,	18.33900452},
			{1,	18.36478043},
			{1,	18.38741112},
			{1,	18.41271973},
			{1,	18.43644333},
			{1,	18.46075630},
			{1,	18.48509216},
			{1,	18.50897980},
			{1,	18.53143501},
			{1,	18.55570030},
			{1,	18.58008003},
			{1,	18.60348320},
			{1,	18.62536430},
			{1,	18.65199852},
			{1,	18.67623520},
			{1,	18.70072174},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 83; _i++)
			a1[_i] = RC_INI(_rs7); }
		{ static struct{ long rc; double ini; } _rs8[] = {
			{1,	1.09462866e10},
			{1,	9.32986675e09},
			{1,	6.15947008e09},
			{1,	1.54486170e09},
			{1,	1.00812454e09},
			{1,	7.00559552e08},
			{1,	6.25999232e08},
			{1,	3.50779968e08},
			{1,	3.11956288e08},
			{1,	3.74866016e08},
			{1,	2.47019424e08},
			{1,	1.73169776e08},
			{1,	1.01753168e08},
			{1,	6.81861920e07},
			{1,	4.61764000e07},
			{1,	3.31671360e07},
			{1,	2.03160540e07},
			{1,	1.40249480e07},
			{1,	1.02577860e07},
			{1,	3.53822650e06},
			{1,	1.32563388e06},
			{1,	9.14284688e05},
			{1,	1.25230388e06},
			{1,	3.17865156e05},
			{1,	4.76750244e03},
			{1,	4.81107031e03},
			{1,	4.88406152e03},
			{1,	4.80611279e03},
			{1,	4.78843652e03},
			{1,	4.65988477e03},
			{1,	1.26723059e03},
			{1,	1.20825342e03},
			{1,	8.66052612e02},
			{1,	7.76661316e02},
			{1,	7.05279358e02},
			{1,	6.21722656e02},
			{1,	5.46207581e02},
			{1,	4.96247742e02},
			{1,	4.26340118e02},
			{1,	3.96090424e02},
			{1,	3.48429657e02},
			{1,	2.37949142e02},
			{1,	2.14678406e02},
			{1,	1.81019180e02},
			{1,	1.68923676e02},
			{1,	1.45979385e02},
			{1,	1.25311127e02},
			{1,	1.05205528e02},
			{1,	9.39378357e01},
			{1,	7.75339966e01},
			{1,	6.68987427e01},
			{1,	5.53580055e01},
			{1,	5.00100212e01},
			{1,	4.14198608e01},
			{1,	3.46289063e01},
			{1,	3.00775223e01},
			{1,	2.60294399e01},
			{1,	2.26602840e01},
			{1,	2.02123032e01},
			{1,	1.76353855e01},
			{1,	1.47198439e01},
			{1,	1.33078461e01},
			{1,	1.17181997e01},
			{1,	1.04125805e01},
			{1,	9.45785904e00},
			{1,	8.42799950e00},
			{1,	7.62769842e00},
			{1,	6.85484743e00},
			{1,	6.25903368e00},
			{1,	5.75135279e00},
			{1,	5.28468180e00},
			{1,	4.87669659e00},
			{1,	4.57353973e00},
			{1,	4.30108690e00},
			{1,	4.05412245e00},
			{1,	3.83283114e00},
			{1,	3.57902861e00},
			{1,	3.43705726e00},
			{1,	3.26563096e00},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 79; _i++)
			ba1[_i] = (realnum)RC_INI(_rs8); }
		{ static struct{ long rc; double ini; } _rs9[] = {
			{1,	3.07498097e00},
			{1,	2.96334076e00},
			{1,	2.92890000e00},
			{1,	2.89550042e00},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 4; _i++)
			bb1[_i] = (realnum)RC_INI(_rs9); }
		{ static struct{ long rc; double ini; } _rs10[] = {
			{1,	-5.46},
			{1,	-5.51},
			{1,	-5.49},
			{1,	-5.30},
			{1,	-5.29},
			{1,	-5.28},
			{1,	-5.37},
			{1,	-5.33},
			{1,	-5.38},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.55},
			{1,	-5.38},
			{1,	-5.19},
			{1,	-5.14},
			{1,	-5.27},
			{1,	-4.93},
			{1,	-3.64},
			{1,	-3.68},
			{1,	-3.74},
			{1,	-3.78},
			{1,	-3.82},
			{1,	-3.88},
			{1,	-3.40},
			{1,	-3.41},
			{1,	-3.32},
			{1,	-3.32},
			{1,	-3.32},
			{1,	-3.32},
			{1,	-3.31},
			{1,	-3.31},
			{1,	-3.29},
			{1,	-3.29},
			{1,	-3.27},
			{1,	-3.16},
			{1,	-3.14},
			{1,	-3.11},
			{1,	-3.10},
			{1,	-3.07},
			{1,	-3.03},
			{1,	-2.99},
			{1,	-2.96},
			{1,	-2.91},
			{1,	-2.87},
			{1,	-2.81},
			{1,	-2.78},
			{1,	-2.72},
			{1,	-2.66},
			{1,	-2.61},
			{1,	-2.56},
			{1,	-2.51},
			{1,	-2.47},
			{1,	-2.42},
			{1,	-2.35},
			{1,	-2.31},
			{1,	-2.26},
			{1,	-2.21},
			{1,	-2.17},
			{1,	-2.12},
			{1,	-2.08},
			{1,	-2.03},
			{1,	-1.99},
			{1,	-1.95},
			{1,	-1.91},
			{1,	-1.87},
			{1,	-1.84},
			{1,	-1.81},
			{1,	-1.78},
			{1,	-1.75},
			{1,	-1.71},
			{1,	-1.69},
			{1,	-1.66},
			{1,	-1.62},
			{1,	-1.60},
			{1,	-1.60},
			{1,	-1.60},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 83; _i++)
			x1[_i] = RC_INI(_rs10); }
		{ static struct{ long rc; double ini; } _rs11[] = {
			{1,	20.30049515},
			{1,	20.28500366},
			{1,	20.25300407},
			{1,	20.16626740},
			{1,	20.15743256},
			{1,	20.11256981},
			{1,	20.04818344},
			{1,	19.99261856},
			{1,	19.94472885},
			{1,	19.86478043},
			{1,	19.83321571},
			{1,	19.80185127},
			{1,	19.74884224},
			{1,	19.70136070},
			{1,	19.65981102},
			{1,	19.60598755},
			{1,	19.56017494},
			{1,	19.52042389},
			{1,	19.47429657},
			{1,	19.44413757},
			{1,	19.40796280},
			{1,	19.34819984},
			{1,	19.32203293},
			{1,	19.27634430},
			{1,	19.25627136},
			{1,	19.22009087},
			{1,	19.18853378},
			{1,	19.14809799},
			{1,	19.12456703},
			{1,	19.08409119},
			{1,	19.05431557},
			{1,	19.02083015},
			{1,	19.00176430},
			{1,	18.96817970},
			{1,	18.93762589},
			{1,	18.91706085},
			{1,	18.89299583},
			{1,	18.87085915},
			{1,	18.85210609},
			{1,	18.83035851},
			{1,	18.80403900},
			{1,	18.78901100},
			{1,	18.77099228},
			{1,	18.75540161},
			{1,	18.74287033},
			{1,	18.72928810},
			{1,	18.71601868},
			{1,	18.70474434},
			{1,	18.69515800},
			{1,	18.68782425},
			{1,	18.68120766},
			{1,	18.67630005},
			{1,	18.67357445},
			{1,	18.67129898},
			{1,	18.67042351},
			{1,	18.67090988},
			{1,	18.67313004},
			{1,	18.67636490},
			{1,	18.68120003},
			{1,	18.68803024},
			{1,	18.69487381},
			{1,	18.70458412},
			{1,	18.71205139},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 63; _i++)
			a2[_i] = RC_INI(_rs11); }
		{ static struct{ long rc; double ini; } _rs12[] = {
			{1,	1.01078403e00},
			{1,	1.97956896e00},
			{1,	3.14605665e00},
			{1,	6.46874905e00},
			{1,	3.16406364e01},
			{1,	3.74927673e01},
			{1,	4.75353088e01},
			{1,	5.27809143e01},
			{1,	5.86515846e01},
			{1,	6.70408707e01},
			{1,	1.14904137e02},
			{1,	1.03133133e02},
			{1,	1.26508728e02},
			{1,	1.03827606e02},
			{1,	8.79508896e01},
			{1,	7.18328934e01},
			{1,	6.19807892e01},
			{1,	5.51255455e01},
			{1,	4.87156143e01},
			{1,	4.58579826e01},
			{1,	4.19952011e01},
			{1,	4.08252220e01},
			{1,	3.78243637e01},
			{1,	3.34573860e01},
			{1,	3.19036102e01},
			{1,	2.92026005e01},
			{1,	2.74482193e01},
			{1,	2.54643936e01},
			{1,	2.46636391e01},
			{1,	2.33054180e01},
			{1,	2.23069897e01},
			{1,	2.12891216e01},
			{1,	2.06667900e01},
			{1,	1.96430798e01},
			{1,	1.87381802e01},
			{1,	1.76523514e01},
			{1,	1.69235287e01},
			{1,	1.62647285e01},
			{1,	1.56806908e01},
			{1,	1.50346069e01},
			{1,	1.42240467e01},
			{1,	1.37954988e01},
			{1,	1.31949224e01},
			{1,	1.27211905e01},
			{1,	1.22885675e01},
			{1,	1.17868662e01},
			{1,	1.12577572e01},
			{1,	1.08565578e01},
			{1,	1.04121590e01},
			{1,	1.00410652e01},
			{1,	9.64534473e00},
			{1,	9.29232121e00},
			{1,	8.92519569e00},
			{1,	8.60898972e00},
			{1,	8.31234550e00},
			{1,	8.04089737e00},
			{1,	7.74343491e00},
			{1,	7.48133039e00},
			{1,	7.21957016e00},
			{1,	6.94726801e00},
			{1,	6.71931219e00},
			{1,	6.45107985e00},
			{1,	6.28593779e00},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 63; _i++)
			b2[_i] = RC_INI(_rs12); }
		{ static struct{ long rc; double ini; } _rs13[] = {
			{1,	-0.43},
			{1,	-0.75},
			{1,	-0.93},
			{1,	-1.20},
			{1,	-1.78},
			{1,	-1.85},
			{1,	-1.95},
			{1,	-2.00},
			{1,	-2.05},
			{1,	-2.12},
			{1,	-2.34},
			{1,	-2.31},
			{1,	-2.42},
			{1,	-2.36},
			{1,	-2.31},
			{1,	-2.25},
			{1,	-2.21},
			{1,	-2.18},
			{1,	-2.15},
			{1,	-2.14},
			{1,	-2.12},
			{1,	-2.14},
			{1,	-2.12},
			{1,	-2.09},
			{1,	-2.08},
			{1,	-2.06},
			{1,	-2.05},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{1,	-2.04},
			{0L, 0 }
		};
		for(_i=_r=0L; _i < 63; _i++)
			x2[_i] = RC_INI(_rs13); }
	/*lint +e736 loss of precision in assignment in translated data */

	DEBUG_ENTRY( "blkdata0()" );
}

/* =================================================================== */
/*xmap mapping function for Cota's 3-body recombination */
STATIC double xmap(double x[], 
	  double y[], 
	  double xmapx0)
{
	double a, 
	  b, 
	  c, 
	  xmapx1, 
	  x12m, 
	  x13m, 
	  xmapx2, 
	  x3, 
	  xmap_v, 
	  yinit, 
	  y13m;

	DEBUG_ENTRY( "xmap()" );

	/* PARABOLIC INTERPOLATION.
	 * */

	yinit = y[0];
	xmapx1 = x[0];
	xmapx2 = x[1];
	x3 = x[2];
	x13m = xmapx1 - x3;
	x12m = xmapx1 - xmapx2;
	y13m = yinit - y[2];
	x3 = (xmapx1 + x3)*x13m;
	xmapx2 = (xmapx1 + xmapx2)*x12m;
	b = ((yinit - y[1])*x3 - y13m*xmapx2)/(x12m*x3 - x13m*xmapx2);
	a = (y13m - x13m*b)/x3;
	c = yinit - a*xmapx1*xmapx1 - b*xmapx1;

	xmap_v = a*xmapx0*xmapx0 + b*xmapx0 + c;

	return( xmap_v );
}

/* =================================================================== */
/*xinvrs do inverse function for Cota's three-body recombination */
STATIC double xinvrs(double y, 
	  double a, 
	  double b, 
	  double u, 
	  double v, 
	  long int *ifail)
{
	long int i;
	double bxu, 
	  dfx, 
	  fx, 
	  fxdfx, 
	  x, 
	  xinvrs_v, 
	  xlog, 
	  xx;
	static long itmax = 10;

	DEBUG_ENTRY( "xinvrs()" );

	/* inverts equation of the form :
	 *
	 *   Y = A + B * X ** U - V * LOG ( X )
	 * */
	*ifail = 0;
	xlog = (a - y)/v;
	x = pow(10.,xlog);
	xx = 0.;

	for( i=0; i < itmax; i++ )
	{
		bxu = b*pow(x,u);
		fx = y - a - bxu + v*xlog;
		dfx = v*.4342945 - bxu*u;

		if( dfx != 0. )
		{
			fxdfx = fabs(fx/dfx);
			fxdfx = MIN2(0.2,fxdfx);
			xx = x*(1. - sign(fxdfx,fx/dfx));
		}
		else
		{
			/* >>chng 96 feb 02 this added in case dfx ever 0
			 *  suggested by Peter van Hoof */
			xx = x*(1. - sign(0.2,fx));
		}

		if( (fabs(xx-x)/x) < 1.e-4 )
		{
			xinvrs_v = xx;
			return( xinvrs_v );
		}
		else
		{
			x = xx;
			if( x < 1e-30 )
			{
				xinvrs_v = 100.;
				*ifail = 1;
				return( xinvrs_v );
			}
			xlog = log10(x);
		}
	}
	xinvrs_v = xx;
	*ifail = 1;
	return( xinvrs_v );
}