/home66/gary/public_html/cloudy/c08_branch/source/atmdat_3body.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 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 "punch.h"
#include "atmdat.h"

#define MAXZ    28

STATIC void blkdata1(void);
STATIC double da(double z);

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));
        }

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

/* =================================================================== */
STATIC double da(double z)
{
        /*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);
        alogne = log10(dense.eden);
        alogte = log10(phycon.te);
        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. )
                {
                        ionbal.ilthn += 1;
                        alogn = 9.;
                }

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

        }
        else
        {
                da_v = 0.;
                return( da_v );
        }

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

        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 );
}

---