16STATIC double da(
double z,
double temp,
double eden);
34#define RC_INI(rs) (rs[_r].rc>1 ? DEC_RC_(rs) : (rs[_r].rc==1 ? INC_NDX_(rs) : rs[_r].ini ))
35#define DEC_RC_(rs) (rs[_r].rc--,rs[_r].ini)
36#define INC_NDX_(rs) (_r++,rs[_r-1].ini)
62 for( i=0; i <
LIMELM; i++ )
80 for( i=0; i < iup; i++ )
89 fprintf(
ioQQQ,
" 3BOD rate:" );
90 for( i=1; i <= 14; i++ )
94 fprintf(
ioQQQ,
"\n" );
100 fprintf(
save.
ioRecom,
" 3-body rec coef vs charge \n" );
101 for( i=0; i < iup; i++ )
111STATIC double da(
double z,
double temp,
double eden)
143 double alogte , alogne;
176 double eden_limited =
MIN2( eden, 1e5 );
177 alogne = log10(eden_limited);
178 alogte = log10(temp);
182 alogt = alogte - 2.*zlog;
183 alogn = alogne - 7.*zlog;
187 alogt = alogte -
zlog2[nz-1];
188 alogn = alogne -
zlog7[nz-1];
208 if( alogt <= 2.1760913 )
210 if( alogn < (3.5*alogt - 8.) )
214 else if( alogn > (3.5*alogt - 2.) )
217 alogn = 3.5*alogt - 2.;
221 else if( alogt <= 2.4771213 )
231 else if( alogt <= 5.1139434 )
244 double alogt_limited =
MIN2( alogt, 5.0 );
245 double alogt_save = alogt;
246 alogt = alogt_limited;
251 nt = (long)(9.9657843*alogt + 1.);
255 nt = (long)(19.931568*alogt - 19.);
261 if( fabs(alogt-
tz[nt-1]) >= fabs(alogt-
tz[
MIN2(83,nt+1)-1]) )
265 else if( fabs(alogt-
tz[nt-1]) > fabs(alogt-
tz[
MAX2(1,nt-1)-1]) )
271 if( fabs(alogt-
tz[nt-1]) < 0.00001 )
290 xnc =
xinvrs(alogn,c,d,u,v,&jfail);
291 if( xnc <= 0. || jfail != 0 )
314 if( (nt <= 21) && (z == 1.0) )
328 if( alogt <= 2.1760913 )
339 if( alogt <= 2.4771213 )
360 yyb[0] = 6.93410742e03;
366 yyb[0] = 1.36653821e03;
378 c =
xmap(xx,yya,alogt);
381 d =
xmap(xx,yyb,alogt);
384 u =
xmap(xx,yyx,alogt);
393 yyb[0] = 2.25774918e01;
399 yyb[0] = 6.70408707e02;
404 yya[0] =
a2[nt0-(21)];
405 yyb[0] =
b2[nt0-(21)];
406 yyx[0] =
x2[nt0-(21)];
409 yya[1] =
a2[nt-(21)];
410 yya[2] =
a2[nt1-(21)];
411 c =
xmap(xx,yya,alogt);
412 yyb[1] =
b2[nt-(21)];
413 yyb[2] =
b2[nt1-(21)];
414 d =
xmap(xx,yyb,alogt);
415 yyx[1] =
x2[nt-(21)];
416 yyx[2] =
x2[nt1-(21)];
417 u =
xmap(xx,yyx,alogt);
421 xnc =
xinvrs(alogn,c,d,u,v,&jfail);
422 if( xnc <= 0. || jfail != 0 )
433 yya[0] = -8.04963875;
434 yyb[0] = 1.07205127e03;
439 yya[0] = -8.54721069;
440 yyb[0] = 4.70450195e02;
452 a =
xmap(xx,yya,alogt);
455 b =
xmap(xx,yyb,alogt);
459 y =
xmap(xx,yyy,alogt);
462 expp = a - y*alognc + b*pow(xnc,x);
465 da_v = z*pow(10.,expp);
473 da_v *= pow( eden/eden_limited, 0.25 );
474 da_v *= pow( 10., 2.*(alogt_limited-alogt_save) );
498 {
static struct{
long rc;
double ini; } _rs0[] = {
529 for(_i=_r=0L; _i < 28; _i++)
531 {
static struct{
long rc;
double ini; } _rs1[] = {
562 for(_i=_r=0L; _i < 28; _i++)
564 {
static struct{
long rc;
double ini; } _rs2[] = {
650 for(_i=_r=0L; _i < 83; _i++)
652 {
static struct{
long rc;
double ini; } _rs3[] = {
738 for(_i=_r=0L; _i < 83; _i++)
740 {
static struct{
long rc;
double ini; } _rs4[] = {
822 for(_i=_r=0L; _i < 79; _i++)
824 {
static struct{
long rc;
double ini; } _rs5[] = {
831 for(_i=_r=0L; _i < 4; _i++)
833 {
static struct{
long rc;
double ini; } _rs6[] = {
919 for(_i=_r=0L; _i < 83; _i++)
922 {
static struct{
long rc;
double ini; } _rs7[] = {
1008 for(_i=_r=0L; _i < 83; _i++)
1010 {
static struct{
long rc;
double ini; } _rs8[] = {
1092 for(_i=_r=0L; _i < 79; _i++)
1094 {
static struct{
long rc;
double ini; } _rs9[] = {
1101 for(_i=_r=0L; _i < 4; _i++)
1103 {
static struct{
long rc;
double ini; } _rs10[] = {
1189 for(_i=_r=0L; _i < 83; _i++)
1191 {
static struct{
long rc;
double ini; } _rs11[] = {
1257 for(_i=_r=0L; _i < 63; _i++)
1259 {
static struct{
long rc;
double ini; } _rs12[] = {
1325 for(_i=_r=0L; _i < 63; _i++)
1327 {
static struct{
long rc;
double ini; } _rs13[] = {
1393 for(_i=_r=0L; _i < 63; _i++)
1428 x12m = xmapx1 - xmapx2;
1429 y13m = yinit - y[2];
1430 x3 = (xmapx1 + x3)*x13m;
1431 xmapx2 = (xmapx1 + xmapx2)*x12m;
1432 b = ((yinit - y[1])*x3 - y13m*xmapx2)/(x12m*x3 - x13m*xmapx2);
1433 a = (y13m - x13m*b)/x3;
1434 c = yinit - a*xmapx1*xmapx1 - b*xmapx1;
1436 xmap_v = a*xmapx0*xmapx0 + b*xmapx0 + c;
1459 static long itmax = 10;
1472 for( i=0; i < itmax; i++ )
1475 fx = y - a - bxu + v*xlog;
1476 dfx = v*.4342945 - bxu*u;
1480 fxdfx = fabs(fx/dfx);
1481 fxdfx =
MIN2(0.2,fxdfx);
1482 xx = x*(1. -
sign(fxdfx,fx/dfx));
1488 xx = x*(1. -
sign(0.2,fx));
1491 if( (fabs(xx-x)/x) < 1.e-4 )
STATIC double da(double z, double temp, double eden)
STATIC double xinvrs(double y, double a, double b, double u, double v, long int *ifail)
STATIC void blkdata1(void)
STATIC double xmap(double x[], double y[], double x0)
#define DEBUG_ENTRY(funcname)