GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
enormp.c
Go to the documentation of this file.
1 #include <math.h>
2 
3 
4 double Cdhc_enormp(double x)
5 {
6  double x1, x2, x3, x4, ret_val;
7  static double xp[5] = { 7.7105849500132e-5, -0.00133733772997339,
8  0.0323076579225834, 0.0479137145607681, 0.128379167095513
9  };
10  static double xq[3] = { 0.00301048631703895, 0.0538971687740286,
11  0.375795757275549
12  };
13  static double xr[8] = { -1.36864857382717e-7, 0.564195517478974,
14  7.21175825088309, 43.1622272220567, 152.98928504694,
15  339.320816734344, 451.918953711873, 300.459261020162
16  };
17  static double xs[8] = { 1.0, 12.7827273196294, 77.0001529352295,
18  277.585444743988, 638.980264465631, 931.35409485061,
19  790.950925327898, 300.459260956983
20  };
21  static double xt[5] = { 2.10144126479064, 26.2370141675169,
22  21.3688200555087, 4.6580782871847, 0.282094791773523
23  };
24  static double xu[4] = { 94.153775055546, 187.11481179959,
25  99.0191814623914, 18.0124575948747
26  };
27  double yy1, yy2;
28 
29  x3 = (double)0.564189583547756;
30  x1 = fabs(x);
31 
32  if (x1 <= .5) {
33  x4 = x * x;
34  yy1 =
35  (((xp[0] * x4 + xp[1]) * x4 + xp[2]) * x4 + xp[3]) * x4 + xp[4] +
36  1.;
37  yy2 = ((xq[0] * x4 + xq[1]) * x4 + xq[2]) * x4 + 1.;
38  ret_val = x * (yy1 / yy2);
39  }
40  else if (x1 <= 4.) {
41  yy1 = ((((((xr[0] * x1 + xr[1]) * x1 + xr[2]) * x1 + xr[3])
42  * x1 + xr[4]) * x1 + xr[5]) * x1 + xr[6]) * x1 + xr[7];
43  yy2 = ((((((xs[0] * x1 + xs[1]) * x1 + xs[2]) * x1 + xs[3])
44  * x1 + xs[4]) * x1 + xs[5]) * x1 + xs[6]) * x1 + xs[7];
45  ret_val = 1. - exp(-x * x) * yy1 / yy2;
46  if (x < 0.)
47  ret_val *= -1.0;
48  }
49  else {
50  x2 = x * x;
51  x4 *= 1.; /* huh, what? See original FORTRAN */
52  x4 = 0.5;
53  yy1 = (((xt[0] * x4 + xt[1]) * x4 + xt[2]) * x4 + xt[3]) * x4 + xt[4];
54  yy2 = (((xu[0] * x4 + xu[1]) * x4 + xu[2]) * x4 + xu[3]) * x4 + 1.;
55  ret_val = x3 / x1 - yy1 * x1 / (x2 * yy2);
56  ret_val = 1. - exp(-x2) * ret_val;
57  if (x < 0.)
58  ret_val = -ret_val;
59  }
60 
61  return ret_val;
62 }
#define x
double Cdhc_enormp(double x)
Definition: enormp.c:4