13 static const int MAXZ = 28;
16 STATIC 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++ )
90 for( i=1; i <= 14; i++ )
101 for( i=0; i < iup; i++ )
111 STATIC double da(
double z,
double temp,
double eden)
142 double alogte , alogne;
175 double eden_limited =
MIN2( eden, 1e5 );
176 alogne = log10(eden_limited);
177 alogte = log10(temp);
181 alogt = alogte - 2.*zlog;
182 alogn = alogne - 7.*zlog;
186 alogt = alogte -
zlog2[nz-1];
187 alogn = alogne -
zlog7[nz-1];
207 if( alogt <= 2.1760913 )
209 if( alogn < (3.5*alogt - 8.) )
213 else if( alogn > (3.5*alogt - 2.) )
216 alogn = 3.5*alogt - 2.;
220 else if( alogt <= 2.4771213 )
230 else if( alogt <= 5.1139434 )
243 double alogt_limited =
MIN2( alogt, 5.0 );
244 double alogt_save = alogt;
245 alogt = alogt_limited;
250 nt = (long)(9.9657843*alogt + 1.);
254 nt = (long)(19.931568*alogt - 19.);
260 if( fabs(alogt-
tz[nt-1]) >= fabs(alogt-
tz[
MIN2(83,nt+1)-1]) )
264 else if( fabs(alogt-
tz[nt-1]) > fabs(alogt-
tz[
MAX2(1,nt-1)-1]) )
270 if( fabs(alogt-
tz[nt-1]) < 0.00001 )
289 xnc =
xinvrs(alogn,c,d,u,v,&jfail);
290 if( xnc <= 0. || jfail != 0 )
311 if( (nt <= 21) && (z == 1.0) )
325 if( alogt <= 2.1760913 )
336 if( alogt <= 2.4771213 )
357 yyb[0] = 6.93410742e03;
363 yyb[0] = 1.36653821e03;
375 c =
xmap(xx,yya,alogt);
378 d =
xmap(xx,yyb,alogt);
381 u =
xmap(xx,yyx,alogt);
390 yyb[0] = 2.25774918e01;
396 yyb[0] = 6.70408707e02;
401 yya[0] =
a2[nt0-(21)];
402 yyb[0] =
b2[nt0-(21)];
403 yyx[0] =
x2[nt0-(21)];
406 yya[1] =
a2[nt-(21)];
407 yya[2] =
a2[nt1-(21)];
408 c =
xmap(xx,yya,alogt);
409 yyb[1] =
b2[nt-(21)];
410 yyb[2] =
b2[nt1-(21)];
411 d =
xmap(xx,yyb,alogt);
412 yyx[1] =
x2[nt-(21)];
413 yyx[2] =
x2[nt1-(21)];
414 u =
xmap(xx,yyx,alogt);
418 xnc =
xinvrs(alogn,c,d,u,v,&jfail);
419 if( xnc <= 0. || jfail != 0 )
429 yya[0] = -8.04963875;
430 yyb[0] = 1.07205127e03;
435 yya[0] = -8.54721069;
436 yyb[0] = 4.70450195e02;
448 a =
xmap(xx,yya,alogt);
451 b =
xmap(xx,yyb,alogt);
455 y =
xmap(xx,yyy,alogt);
458 expp = a - y*alognc + b*pow(xnc,x);
461 da_v = z*
exp10(expp);
468 da_v *=
powpq( eden/eden_limited, 1, 4 );
469 da_v *=
exp10( 2.*(alogt_limited-alogt_save) );
493 {
static struct{
long rc;
double ini; } _rs0[] = {
524 for(_i=_r=0L; _i < 28; _i++)
526 {
static struct{
long rc;
double ini; } _rs1[] = {
557 for(_i=_r=0L; _i < 28; _i++)
559 {
static struct{
long rc;
double ini; } _rs2[] = {
645 for(_i=_r=0L; _i < 83; _i++)
647 {
static struct{
long rc;
double ini; } _rs3[] = {
733 for(_i=_r=0L; _i < 83; _i++)
735 {
static struct{
long rc;
double ini; } _rs4[] = {
817 for(_i=_r=0L; _i < 79; _i++)
819 {
static struct{
long rc;
double ini; } _rs5[] = {
826 for(_i=_r=0L; _i < 4; _i++)
828 {
static struct{
long rc;
double ini; } _rs6[] = {
914 for(_i=_r=0L; _i < 83; _i++)
917 {
static struct{
long rc;
double ini; } _rs7[] = {
1003 for(_i=_r=0L; _i < 83; _i++)
1005 {
static struct{
long rc;
double ini; } _rs8[] = {
1087 for(_i=_r=0L; _i < 79; _i++)
1089 {
static struct{
long rc;
double ini; } _rs9[] = {
1096 for(_i=_r=0L; _i < 4; _i++)
1098 {
static struct{
long rc;
double ini; } _rs10[] = {
1184 for(_i=_r=0L; _i < 83; _i++)
1186 {
static struct{
long rc;
double ini; } _rs11[] = {
1252 for(_i=_r=0L; _i < 63; _i++)
1254 {
static struct{
long rc;
double ini; } _rs12[] = {
1320 for(_i=_r=0L; _i < 63; _i++)
1322 {
static struct{
long rc;
double ini; } _rs13[] = {
1388 for(_i=_r=0L; _i < 63; _i++)
1423 x12m = xmapx1 - xmapx2;
1424 y13m = yinit - y[2];
1425 x3 = (xmapx1 + x3)*x13m;
1426 xmapx2 = (xmapx1 + xmapx2)*x12m;
1427 b = ((yinit - y[1])*x3 - y13m*xmapx2)/(x12m*x3 - x13m*xmapx2);
1428 a = (y13m - x13m*b)/x3;
1429 c = yinit - a*xmapx1*xmapx1 - b*xmapx1;
1431 xmap_v = a*xmapx0*xmapx0 + b*xmapx0 + c;
1454 static long itmax = 10;
1467 for( i=0; i < itmax; i++ )
1470 fx = y - a - bxu + v*xlog;
1471 dfx = v*.4342945 - bxu*u;
1475 fxdfx = fabs(fx/dfx);
1476 fxdfx =
MIN2(0.2,fxdfx);
1477 xx = x*(1. -
sign(fxdfx,fx/dfx));
1483 xx = x*(1. -
sign(0.2,fx));
1486 if( (fabs(xx-x)/x) < 1.e-4 )
STATIC double xmap(double x[], double y[], double x0)
STATIC void blkdata1(void)
STATIC double da(double z, double temp, double eden)
#define DEBUG_ENTRY(funcname)
double powpq(double x, int p, int q)
int fprintf(const Output &stream, const char *format,...)
STATIC double xinvrs(double y, double a, double b, double u, double v, long int *ifail)