GRASS Programmer's Manual  6.4.2(2012)
enormp.c
Go to the documentation of this file.
00001 #include <math.h>
00002 
00003 
00004 double enormp(double x)
00005 {
00006     double x1, x2, x3, x4, ret_val;
00007     static double xp[5] = { 7.7105849500132e-5, -0.00133733772997339,
00008         0.0323076579225834, 0.0479137145607681, 0.128379167095513
00009     };
00010     static double xq[3] = { 0.00301048631703895, 0.0538971687740286,
00011         0.375795757275549
00012     };
00013     static double xr[8] = { -1.36864857382717e-7, 0.564195517478974,
00014         7.21175825088309, 43.1622272220567, 152.98928504694,
00015         339.320816734344, 451.918953711873, 300.459261020162
00016     };
00017     static double xs[8] = { 1.0, 12.7827273196294, 77.0001529352295,
00018         277.585444743988, 638.980264465631, 931.35409485061,
00019         790.950925327898, 300.459260956983
00020     };
00021     static double xt[5] = { 2.10144126479064, 26.2370141675169,
00022         21.3688200555087, 4.6580782871847, 0.282094791773523
00023     };
00024     static double xu[4] = { 94.153775055546, 187.11481179959,
00025         99.0191814623914, 18.0124575948747
00026     };
00027     double yy1, yy2;
00028 
00029     x3 = (double)0.564189583547756;
00030     x1 = fabs(x);
00031 
00032     if (x1 <= .5) {
00033         x4 = x * x;
00034         yy1 =
00035             (((xp[0] * x4 + xp[1]) * x4 + xp[2]) * x4 + xp[3]) * x4 + xp[4] +
00036             1.;
00037         yy2 = ((xq[0] * x4 + xq[1]) * x4 + xq[2]) * x4 + 1.;
00038         ret_val = x * (yy1 / yy2);
00039     }
00040     else if (x1 <= 4.) {
00041         yy1 = ((((((xr[0] * x1 + xr[1]) * x1 + xr[2]) * x1 + xr[3])
00042                  * x1 + xr[4]) * x1 + xr[5]) * x1 + xr[6]) * x1 + xr[7];
00043         yy2 = ((((((xs[0] * x1 + xs[1]) * x1 + xs[2]) * x1 + xs[3])
00044                  * x1 + xs[4]) * x1 + xs[5]) * x1 + xs[6]) * x1 + xs[7];
00045         ret_val = 1. - exp(-x * x) * yy1 / yy2;
00046         if (x < 0.)
00047             ret_val *= -1.0;
00048     }
00049     else {
00050         x2 = x * x;
00051         x4 *= 1.;               /* huh, what? See original FORTRAN */
00052         x4 = 0.5;
00053         yy1 = (((xt[0] * x4 + xt[1]) * x4 + xt[2]) * x4 + xt[3]) * x4 + xt[4];
00054         yy2 = (((xu[0] * x4 + xu[1]) * x4 + xu[2]) * x4 + xu[3]) * x4 + 1.;
00055         ret_val = x3 / x1 - yy1 * x1 / (x2 * yy2);
00056         ret_val = 1. - exp(-x2) * ret_val;
00057         if (x < 0.)
00058             ret_val = -ret_val;
00059     }
00060 
00061     return ret_val;
00062 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines