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,...)
double pow(double x, int i)
STATIC double xinvrs(double y, double a, double b, double u, double v, long int *ifail)